The authors of the study examine the surace energy balance of the largest Icelandic glaciers from 2000 to 2021, taking into account the ablation period from April to September. They apply an energy balance estimation using MODIS-derived albedos. High-resolution WRF data with an hourly resolution of 2 km are used as meteorological forcing. The results show the large spatial and temporal variability of the melting energy. The energy balance terms are presented in detail, resolving the different glacier areas, elevations and seasonal and annual patterns. The special feature of the study, in my view, is the coverage of almost the entire ice-covered area (97%) of Iceland and especially the focus on the influence of light-absorbing particles on the energy balance. These particles come from sand deserts and volcanic eruptions.

The paper is excellently written, and the text is easy to understand. It was enjoyable to read the manuscript.

In my opinion, the publication can make a valuable contribution to the knowledge of glaciers in Iceland. However, I have some points that should be (better) addressed before publication. I hope that the following major, specific and technical comments can help improve the manuscript.

Major comments

1. Novelty and differentiation from Gunnarsson et al. (2021): I like the idea of building on a previous study using the derived gap-filled and post-processed MODIS dataset. My main issue is that some of the results and conclusions are similar to those of the 2021 TC paper. You mentioned several times that SWnet is modulated by the albedo and that the melt patterns are mainly driven by SWnet. Therefore, it is obvious that e.g. the elevation gradient of the albedo in Gunnarsson et al. (2021) is consistent with the elevation gradients of the melt energy presented here. Therefore, at the end of the introduction, I would recommend clearly elaborating what the similarities and differences are to Gunnarsson et al. (2021),  how this study extends Gunnarsson et al. (2021) and what makes this study unique.
2. Energy balance: wind speed? surface temperature measured or iteratively solved? humidity at the surface?: The explanations regarding the energy balance estimates are not detailed enough in my opinion. According to Section 3.3, air temperature, surface temperature, incoming long and short wave radiation, barometric pressure and specific humidity are used. The wind speed which is of major importance for the turbulent fluxes, is missing. You could use the output wind fields from WRF. According to equation (5) wind speed is needed? Furthermore, what is used for the surface humidity at height z0? The WRF output or is the surface humidity assumed to be 100 %? Another question I have is the surface temperature. Is the surface temperature used from WRF? According to the text the SEB was solved iteratively for surface temperature (Line 205), but according to Section 3.3, the surface temperature is used from WRF.
3. Uncertainties, simplifications and limitations: I think you are aware of the uncertainties, simplifications and limitations of the study. Nevertheless, in my opinion, these are too little discussed in the present study. One idea would be to add a subsection after the validation and collect and discuss the different issues. If you add a section before presenting and discussing the results you directly show that you are aware of these issues. In the following I will highlight only some of the things I thought of. Sometimes you have mentioned the points I thought of when reading the manuscript, but only at a later stage.
   • If you have a tephra layer above the snow and ice surface you can have surface temperatures > 273.15 K in reality.
   • If I understand correctly, you calculate  the albedo from an 11 day average. So if there is a thin snow cover from a summer snowfall lasting only 2 days, for example, this will be underrated in your approach. It is not essential for the results but I think you still can mention such limitations.
   • There are limitations within the WRF data. Line 162: I would recommend at least discussing the uncertainties when you combine WRF datasets with a totally different forcing (ERA-Interim versus NCEP). In my opinion your combined forcing dataset is not consistent anymore. It is ok to use the combined dataset but this issue has to be mentioned and discussed.
   • Your SEB estimation has limitations and simplifications. Can you name more effects and discuss them including references? E.g. the bulk approach, LWout  parameterisations.
   • In the 2021 TC paper you write: “[...] Vatnajökull, and boundaries in 2007 and 2008 were used for Langjökull and Hofsjökull, respectively. This was selected as a midpoint representing an average glacier area during the period 2000–2019. This needs to be considered when interpreting rapid changes at the glacier terminus, as some areas in 2000 were part of an active glacier but might in 2019 be dead ice or land.” and for example “It is important to consider how representative point-based in situ observations are (observing ∼ 120–180 m2; Kipp and Zonen, 2019), compared with the spatial footprint of the MODIS data (0.25 km2), especially in glaciated areas with high spatial albedo variability and MODIS sub-pixel variability as is observed in the bare-ice areas of the Icelandic glaciers.“ I think such considerations should be made here as well, adapted to the used data.
   • Line 257 a temperature bias up to 1.15 K. I suggest removing slightly, and discuss the uncertainty resulting from the bias. 1 K makes a difference.
4. Validation: Some of the meteorological forcing variables are validated. What about barometric pressure, specific humidity and probably wind speed. How did you downscale the barometric pressure from the 2 km WRF grid to the 463 m MODIS grid? The MODIS data are validated in the 2021 TC paper. That’s good. But I could not find a validation of the calculated energy balance terms (SWin, LWin, LWout, SHF, LHF) and the resulting potential melt energy? In my opinion the validation of the results could be done within the discussion of the results or in a specific subsection of the results and discussion section. I understand that there may not be directly comparable data. But you could use other studies in Iceland on single glaciers, or use studies on the Greenland ice sheet, or in Svalbard, or in Scandinavia to at least assess the range of the calculated values. Furthermore, you could also compare relative values with Björnsson (1971, 1972), and more recent studies. In the abstract you write: “Validation was performed using observations from various glaciers spanning distinct locations and elevations with good visual and statistical agreement.” after the sentence: “ The SEB was reconstructed from April through September for 2000–2021 at a daily timestep with a 500 m spatial resolution.” So I expected a statistical validation of the SEB.
5. Estimated SW radiative forcing from LAPs: Line 443–445: If I understand correctly, you use the same climate forcing e.g. for the year 2010, first with the mean albedo (2000–2021), then with the observed albedo in 2010. Besides LAPs the observed albedo in 2010 could also be influenced by climate, or? In Line 481 you state that for example snowfall has an impact on SWfLAP,. Please explain your setup in detail and discuss this issue. I like the approach and the investigation, it would just be good if you could show that you are aware of the limitations and possible influence of e.g. snowfall and temperatures on the observed albedo in certain years.

Specific comments

Titel: The title 'Modeling of surface energy balance' is very general. Perhaps sth. about LAPs, volcanic impacts, .. could be added.

Abstract: There are no numbers from your results in the abstract. Maybe the mean melt enhancement (in %) from LAPs could be added to the abstract.

Line 7: What is the difference between annual variability and inter-annual variability? By annual variability, do you mean intra-annual variability or seasonal variability? For me, annual variability is the same as inter-annual variability. You wrote seasonal and inter-annual in the heading of section 5.1 and in the conclusion Line 493: “[...] melt-season and inter-annual variability [...]”.

Line 32: Could you explain a little more in depth what “ high precipitation sustaining a seasonal snow pack and glaciers” means. In which month/season is the precipitation peak? Which months are the driest?

Introduction: The introduction is rather long with 1650 words. Please revise this section and check which sentences are really needed for the motivation of the study. The historical background in Line 73–102 is very interesting. Nevertheless, I think these paragraphs can be shortened.

Methods: In contrast to the introduction the Methods section is rather short, especially the presentation of the surface energy balance and the parameterisation of the different terms (cf. Major comment 2). Is there a storage (snow/ice temperature) which is not mentioned or how is the cold content from winter (Line 378-379) resolved by the estimation. The sub-surface heat flux which could transport cold content to the surface is assumed to be zero (Line 209).

Line 132–132: Did you derive the albedo again from the MODIS product? Or did you use the dataset from Gunnarsson et al. (2021). If you used the dataset, I recommend rephrasing the sentence accordingly. 

Line 140–150: The final used albedo product has a daily resolution or? I would name the final used temporal resolution in this paragraph. I was first confused with the 11 days buffer.

Line 176–177: Did you adjust the WRF output to the original IslandsDEM or to the MODIS grid (463 m)?

Line 187: What does “environmental” lapse rate mean in this context?

Line 201: Due to the usage of 5 days backward/forward in case of the MODIS data and the original hourly WRF forcing and three different spatial resolutions it was not directly clear for me what is the final temporal and spatial resolution. Maybe you can add to the sentence in Line 201: “[...] using estimations of daily SEB with a resolution of 463 m.” if I am correct. Or sth. similar indicating the final spatio-temporal resolution.

Section 4.1: Is there a reason why you present R2 for T2 and LW but not for SW?

Line 249: Maybe there is a misunderstanding from my part, but for me SW, LW, .. are the energy balance components. So maybe you mean: “The downscaled meteorological forcing [...]” instead of “The downscaled energy balance components [...]”.

Results and discussion: Sometimes it is difficult to recognise which are results of the study and which are results evaluated with the help of other studies. Separating the results and discussion into two different sections would help here. With this, the discussion could also be conducted more independently of the order of the graphs. Furthermore, the discussion could be expanded. Especially the comparison with other studies with numbers would be helpful. This comment adds to the validation of the energy balance terms (major comment 4). The comparison in Line 316–320 is very general and all studies are cited at the end of the paragraph. Readers will be interested in a more in depth comparison of what is similar and what is different. Besides the calculated energy balance terms and the available melt energy the gradients could be compared (Line 310–315). Furthermore, you can discuss that you found positive albedo trends over the study period in northern Vatnajökull in the TC 2021 paper, but no significant trends were found in this study.

Line 285: In my opinion, you cannot see the inter-annual variability with Figure 3. You can see the seasonal and spatial variability. But extreme positive or negative years are not visible.

Line 291: I understand between 10 and 15 % of the mean annual (2000-2021) melt energy was observed. If so, think of adding ‘mean annual’.

Line 321–322: Please add a reference to “other Northern Hemisphere glaciers and ice sheets”.

Line 324–325: Can you add a short statement how a negative correlation between LWnet and SWnet increases the contribution of the sensible heat fluxes? You mean the relative contribution?

Line 328: How do you know that the albedo was mainly driven by climatology? From the applied method or from another source. I recommend adding a short explanatory sentence or a reference. In the MODIS data you just see the evolution of the albedo, but in the first place you do not see the reason, for example, for a sudden decrease.

Line 347–348: Can you add a reference to “warm southerly winds and precipitation” and I guess you mean liquid precipitation or? So maybe add “liquid” to precipitation or change precipitation to “rain”.

Line 344: I think you can partly restructure the discussion. Here you are already talking about the impacts of volcanic eruptions before the subsection “5.2 Impacts of volcanic eruptions and other LAP events” starts. 

Line 350: How do you know that the LAP deposits are from the near pro-glacial areas and not from somewhere else?

Line 352–354: You probably got the information about “clear skies” and “cold temperatures” from the weather stations. But where does the information about the winds come from? I could not find the information in the manuscript. Please add somewhere a sentence with reference, maybe in the methods.

Line 372: Do you have an assumption or can you discuss why cloud cover and LWnet were not significantly correlated?

Line 385: When you cite explicit numbers, a direct reference would be good. Please add a reference to 0.06 km3. The same applies to Line 390.

Line 442: “The impacts … 2004, 2010, 2011 and 2019 were assessed..” and in Line 444 “observed albedo in 2010, 2011 and 2019.” How was 2004 assessed? Using the observed albedo in 2004 or 2005? Because the event was in fall 2004 if I understood correctly. Please specify. 

Line 443 and Figure 9: Inconsistent. In the text you write: “mean albedo for the study period (2000-2021)” in the caption you write: “average albedo (2000–2021 mean excluding 2010, 2011 and 2019 in the mean)”.

Conclusion: One of your conclusions is the influence by high climate variability. To support this statement you could create some monthly (2000–2021) and annual plots of the different forcing variables placed in the appendix. These plots could also support the discussion of the climatic influences on the melt patterns in the “Results and discussion”’ section.

Technical, minor comments

Line 30: I think the dot within 103.000 is wrong. I assume you wanted to use a thousands separator. If so, you should rather use a comma. The same applies to Line 44 and Line 49. Please also check the other sections. If you want to use a thousands separator you should also use it everywhere. For example in Line 44 “3400” and Line 221 “1005”.

Line 174: I recommend adding the url to the dataset here: “(https://www.lmi.is/, last access: June 1, 2020)”. And maybe you have a link that points specifically to the dataset and not to the main page.

Line 179 and all further units of Kelvin: “6–7° K” should be “6–7 K” without the degree sign. Furthermore here you do not use a space before the degree sign, while placing a space in Line 181.

Line 179: Here you write “6–7 K”.. in Line 181 you write “4.5 K km-1 to 8 K km-1”. I would try to be consistent throughout the manuscript and either write “6–7” or “6 to 7”.

Line 253: "reported by Schmidt et al. (2017)”. Only the year in parentheses.

Equation (1) and LIne 208: Maybe use the common abbreviations HS and HL, or QH and QL (sometimes QE), or SH and LH for the sensible and latent heat flux here and elsewhere in the manuscripts and plots.

Line 435: Please remove space after “respectively” and before the comma.

Line 513 (Data availability): The reference NLSI (2019) is missing in the bibliography.

Figures and tables:

Most people know what T2, SWin and LWin mean. Nevertheless all tables and figures should be completely readable without the main manuscript. Therefore, it would be good to explain all abbreviations in all captions if they are not explained in a legend within the plot. E.g. T2 SWnet, LWnet SHF, LHF.

Figure 1: Can you add the latitude and longitude to the axes? The Vatnajökull map is missing the full glacier name. Some of the weather station names are not readable and can only be assigned to Table A1 by excluding the others. A scale in all maps would also be handy.

Figure 2: Please add to the caption a short description of what can be seen in the different rows (RAV2, ICEB and FCST) of the tables. And perhaps a reference to the table with the additional statistics can be added.

Figure 3: These sentences are redundant: “Note that the colour scale varies between months. Note that the scale varies between panels.” The second one would be enough. Is it impossible to see anything in e.g. April if the same scaling is used for all panels? I understand the problem, but it would be extremely helpful using the same scale for all panels visualising the seasonal evolution and to support the statements made in section 5.1.

Figure 4: The vertical scale varies between the panels as well. The label intervals of the x-axis are random. The LWnet panel for Hofsjökull has 5 intervals in a range of 8 W in 2 W steps. The others, for example Drangajökull, have only 2 intervals in the range of 5 W: minimum and maximum for LWnet. I think the different scales of the x-axis are chosen for a reason. If it is possible otherwise, I think this would be preferable. The LWnet gradient for Hofsjökull looks steeper than that of Vatnajökull. This is only due to the scaling of the axis.

Figure 5: Have the grid points that are in one bin been weighted in any way? Or do some bins only consist of, for example, two grid points and others bins of 300?

Figure 6: Please add somewhere in the caption “mean monthly” and an explanation of SWnet, LWnet, SHF, and LHF. Furthermore you can try to use thicker lines for melt energy, albedo and cloud cover to increase their visibility.

Figure 7: I like the box plots, but an explanation of what we can see would be good. There are different variants of box plots. Which percentiles, mean/median, ….. Furthermore, the caption needs a short statement explaining SWnet, LWnet, SHF, and LHF. Here the information that the scale varies is missing.

Table 2: Are all results statically significant? In the 2021 TC paper you explained statically significance using the p value. Maybe you add the p value presentation to the manuscript as well.

Table A1: I recommend writing the full names of the column headings in the caption. Elevation (Ele.), Number of air temperature in 2 m measurements (N. T2 obs.),...