Summary

This manuscript presents an overview over the performance of ERA5, ERA5-LAND and dynamical downscaling with the WRF model (to dx~3km and dx~1km) for calculating the surface energy balance (SEB) of mountain glaciers over Western Canada. Four glaciers are chosen for evaluation, where observations of the relevant variables (e.g., turbulent heat fluxes, temperature, radiation, wind speeds, etc.) are available during the summer season. The authors derive the simulated variables for the SEB, after some corrections, from the model output and directly compare the results with the observed values. Furthermore, they run the WRF model in multiple configurations for paramterizations to find the "optimal" setup for a satisfactory calculation of the SEB. Results suggest that dynamical downscaling with WRF does not automatically outperform ERA5, except for the wind speed and direction - mostly due to the higher horizontal resolution. Generally speaking, both ERA5 and WRF are useful for calculating the SEB, while a corect simulation of the meteorological fields over the glaciers would require even higher horizontal resolution at the hectometric range.

The manuscript is extensive and has a valuable purpose in discussing the challenges of dynamical downscaling over glaciated environments and suggesting an "optimal" setup for future applications. However, in some sections, the authors need to argue in more detail on why they apply a new method; some reasonings are given in the discussion, while they would be already required in the methods section. The interpretation the the WRF  results is sometimes lacking an important factor - namely terrain resolution. Comments and suggestions are given in the list below.

Major comments

Calculation of the surface fluxes from model output via the bulk method. I agree that the modelled albedo values strongly differ from the observations; however, while only reading the methods it is difficult to follow the argumentation why the authors decide to calulate the turbulent fluxes with the observed albedo via the bulk method instead of directly using the values for sensible & latent heat fluxes from model output. Is this a common method to utilize output from an atmospheric model for glacier SEB modelling- was this approach also used in previous studies?

The authors mention in the discussion the unsatifactory performance from the turbulent fluxes from the direct model output (lines 497--508), but for the general understanding of the manuscript, it would make sense to add these paragraphs directly after they introduce the new method (ca. line 385).

Furthermore, changing one parameter to derive a quantitity from the rest of the modelled output might lead to physical inconsistencies, beause all the other variables used for the bulk method still indirectly depend on the "wrong" albedo. Did the authors calculate the SEB with directly modelled turbulent fluxes?

Interpretation - terrain resolution. The authors argue that the poor performance of wind speed and direction simulation yields from the inability to simualte the katabatic galcier wind. The authors could check whether the "bad" model performance only happens during the wind directions coresponding to the down-glacier wind - the model seems to perform better during synoptically-forced conditions. However, glacier winds are not the only meteorological phenomenon present over mountain glaciers; such as thermally-induced circulations, downslope windstorms, etc, which are all mostly governed by the topography (Goger et al, 2022). Therefore, well-resolved topography is essential for the correct simulation of the wind field - tis also explains the general bias reduction of wind speed & direction for small horizontal grid spacings (dx=1.1km and dx=370m). This is an important point which should be mentioned in the discussion and interpretation of the results. Publications from idealized simulations argue that at least 10 points across a valley are necessary to simulate the relevant processes well, and that the correct representation of topography is likely more important than the choice of parameterization schemes (Wagner et al, 2014).

TOPSIS and minRMSE configurations. Maybe I missed it, but do the authors somewhere list the final WRF model setup of TOPSIS and minRMSE, like Table 2 for the REF run? This might be of use for furutre dynamical downscaling studies.

Minor comments

line 50: which simplified assumptions?

line 57: make a new paragraph

line 83: An extensive analysis of real-case, high-resolution large-eddy simulations over a glacier is provided by Goger et al (2022), and Sauter & Galos (2016) performed semi-idlealized LES over a glacier and evaluated the calculation of turbulent fluxes. 

line 85: "Downscaling to several kilometers": Several kilometers might not be the optimal target for mountain glaciers embedded in highly complex terrain, which requires likely horzintal grid spacings of less than 1km.

lines 134-203: I understand that it is important to mention the most commonly used parameterization schemes in WRF, but this is too lengthy for an introduction - perhaps it's enough to mention this configuration in the methods and finally say how it performs within the ensemble.

line 213: You can place the optimal configuration of parameterizations from the introduction here.

line 220: What do you mean exactly by "reflect different time windows during melt season"? 

line 398: " none of these altered WRF configurations yield a strong impact on the calculated Q_M from the SEB model": Did you reset the albedo for calculating the turbulent fluxes here as well?  Becuase then this relative agreement is not very surprising.

line 478: ...."do not distinguish between ice and snow categories": It's true that the land use category does not distinguish between snow and ice. However, after intialization, WRF indeed initalizes snow cover on glacierzed surfaces. The authors mention observed snow cover at one of the glaciers during the time window of interest - is this snow cover present in WRF as well? If yes, the snow cover indeed has an influence on the SEB in the model. 

Figures 8 and 10: Please add a background grid to the figure, this improves their readability.

References

Goger, B., Stiperski, I., Nicholson, L., and Sauter, T. (2022): Large-eddy simulations of the atmospheric boundary layer over an Alpine glacier: Impact of synoptic flow direction and governing processes, Q. J. R. Meteorol. Soc, 148, 1319–1343, https://doi.org/10.1002/qj.4263

Sauter, T. and Galos, S. P. (2016): Effects of local advection on the spatial sensible heat flux variation on a mountain glacier, The Cryosphere, 10, 2887–2905, https://doi.org/10.5194/tc-10-2887-2016

Wagner, J. S., A. Gohm, and M. W. Rotach (2014): The impact of horizontal model grid resolution on the boundary layer structure over an idealized valley. Mon. Wea. Rev., 142, 3446–3465, https://doi.org/10.1175/MWR-D-14-00002.1

This paper investigates the use of ERA5 reanalysis data set, with and without dynamical downscaling to evaluate surface energy balance modeling of glaciers in western Canada. The authors look specifically at the variables that are used to calculate the energy available for melt when evaluating the different forcing sets. For downscaling they use the Weather Research Forecasting (WRF) model, and include several tests with various parametrization combinations in WRF and try to determine the best combination for the western region in Canada.  I believe this paper is a nice contribution to anyone interested in glacier mass balance studies in this and similar regions and should be published after addressing the comments and suggestions.

Major Comment

Page 12, line 257. Does the Fitzpatrick 2019 paper actually support neglection of heat flux from rain (or ground?) The only thing I could find was that “while readings from periods affected by precipitation on the analyser windows were removed” for eddy covariance data. That does not warrant negligible contribution from rain. I suggest removing reference to the 2019 paper.   However, the Fitzpatrick 2017 paper is a nice citation. And I found this paragraph interesting: “QR provided a negligible contribution (<1%) to the total melt energy over the recorded period. However, over daily and sub-daily timescales, QR was observed to have a considerable influence on SEB and ablation during heavy rainfall.”

Precipitation is one of the outputs you want to get correct in WRF. Though it might not have been important for the study window in this manuscript, I do believe a WRF parameterization combination that gets the precipitation correct (along with other metrics) is desired (and as stated in the Fitzpatrick 2017 paper, “QR was observed to have a considerable influence on SEB and ablation during heavy rainfall”). I would have liked to see which parameter combinations score the best when precipitation also is included 

Minor Comments:

Page 1, line 22: I suggest rephrasing to “…and are increasingly losing a considerable amount of mass …..”

Page 2, line 39: Why does SEB models not require precipitation?

Page2, line 43: I suggest rephrasing to “ …fewer than 100 sites worldwide, and only a handful in Western Canada….”

Page 3, section starting on line 62. Here I suggest including a citation of the recent work by Eidhammer et al 2021 (https://hess.copernicus.org/articles/25/4275/2021/), where they use the detailed snow model Crocus within the WRF-Hydro model to estimate glacier melt (and streamflow). They used 1 km downscaled WRF simulations over a glacier in Norway for four seasons. 

Page 3: In the discussion of using dynamical downscaling, you might want to add a comment related to the paper by Lundquist et al. 2019 with the title: “Our Skill in Modeling Mountain Rain and Snow is Bypassing the Skill of Our Observational Networks”. https://doi.org/10.1175/BAMS-D-19-0001.1. I think that this can add to the argument in this manuscript to use downscaling for SEB modeling

Page 3, line 22. I suggest adding a citation to the paper by Liu et al 2011 “High-Resolution Simulations of Wintertime Precipitation in the Colorado Headwaters Region: Sensitivity to Physics Parameterizations “ (https://doi.org/10.1175/MWR-D-11-00009.1). They tested several different WRF physics parameterizations over the Colorado headwaters region

Page 3, line 94. The Thompson-Eidhammer scheme (https://doi.org/10.1175/JAS-D-13- 0305.1) has also been used for Glacier studies

Page 5, Table 5, caption: What is meant by “full days”?

Page 5, line 147: I suggest rephrasing: “…the accumulation zone of the Conrad glacier in 2016. “

Page 6, line 163: The way I read the sentence, the reference to Table 1 indicates that there is some information in regards to the melting surface with intermittent fresh snowfall in the Table 1. I do not think the reference to Table 1 is necessary here. 

Page 8, lines 2 and 3: I suggest to clarify that both 3D and 2D ERA 5 fields (I assume some 2D fields are used at initialization) are used as forcing data for the WRF model.

Page 8, line 193. The way I read this line, the d1 domain for all the 4 glaciers are the same. However, Figure 3 shows that d1 is different between Kaskawulsh and the other glaciers. Can you please clarify?

Page 8, line 198. I assume that you mean that many physics variables are updated every 2.2 s, not outputted? And most likely, the radiation and land surface variables are probably not calculated every 2.2 seconds, but perhaps somewhere between every 5 or 30 minutes? Also make sure if the hourly outputs indeed are hourly averages. Typically, most of the WRF outputs are instantaneously outputs, with some of them being accumulated.

 Page 8, line 199: Table 2 does not describe any of the output saved. I suggest remove the reference to Table 2

Page 9, line 208. I wonder if it would be helpful to add a delta elevation from AWS in table S1 as well. In this way, it would be easier to see the actual elevation difference.

Page 9, line 214. On page 3, it is stated that the Microphysics by Morrison is the most commonly used in glacier studies, but in this work, the Thompson microphysics is used. Please clarify this discrepancy. Also see page 27, line 597

Page 13, line 277: It is stated that the goal is to evaluate daily timeseries of simulated energy available for melt. As shown in Fitzpatrick 2017, the QR was shown to have considerable influence on SEB when considering daily and sub-daily timscales. I am wondering if ignoring QR in this study is then valid?

Page 13, line 286. The study by Eidhammer et al. 2021 also shows that the albedo, when using Noah-MP does not perform well over glaciers (especially in the ablation region).

Page 18, line 371: Please specify that this is the WRF REF data.

Page 19, section 3.3 I think that this section should come before section 3.2, since you are using the results from section 3.3 in describing results in section 3.2

Page 19, line 404: Level 3, not Lavel 3

Page 23, line 490: Did you consider employing different roughness lengths over snow versus ice?