My comments are given in the attached .pdf file.

This paper describes new modeling of Greenland snow and ice in the ORCHIDEE surface model. It is well written and pleasant to read. The work is substantial and certainly deserves to be published. The idea of using MAR as a forcing and comparing it to its outputs is very pertinent. Nevertheless, it seems to me that certain modeling choices need to be discussed and model validation a little improved.

Major comments :

1 - For Figures 1 to 4, you give the raw distributions of MAR and each simulation. But to show differences would be also interesting. In addition, for the others (Figure 5 to 8), you give the differences, but not the raw distributions. As it is possible for you to have additional information, it would be nice to provide all the details, i.e. both the raw distributions and the differences for each variable analyzed (and thus for each figure 1 to 8).

2 - By the way, regarding this validation, why not provide the PDF of each quantity, i.e. compare the MAR (or MODIS) PDF with the ORCHIDEE PDFs ? When you look at Figures 1 to 8, it's hard to see whether one version is much better than another. An objective method of comparison is missing. Comparing PDFs could be a solution but there is perhaps another ways.

3 - It seems that spatial statistics (correlation, etc.) are missing for each field analyzed. In other words, two or three objective criteria to determine whether in Figure 2 (for example) the OPT-12L spatial distribution (e) is better than the others compared to MAR (a). That is, each panel (b, c , d, e) should have its spatial correlation (spatial rmse, etc.) with MAR. The fact that, for example, the spatial distribution of OPT-3L refreezing (Figure 3e) is closer to MAR (Figure 3a) is not trivial to see with the naked eye. Anyway, I hope this comment is understandable. You lack objective statistical criteria in your assessment of all the figures showing comparisons of spatial distribution. The simple Table 2 obtained via a spatial average is not enough.

4 - The modeling choices made could have been discussed. In particular, the parametrization of snow albedo. Moreover, I don't understand why this new parametrization compared to Wang et al. (2013) is in section 2.1 (existing parm) and not rather in section 2.2 (new param). This new parametrization is a bit outdated today when there are more robust parametrizations in land surface models accounting for spectral albedo and solar absorption calculation as in CLM with SNICAR (Flanner and Zender 2006) or ISBA-ES (Decharme et al. 2016). What's more, this more robust representation already exists for ES (Decharme et al. 2016). Why not use it ? This choice is debatable in view of the importance of albedo on the SMB. Please discuss about that.

Other comments:

1 - Independent MAR observations of the Greenland SMB based on GRACE data could have been used in section 5.3 (Schlegel et al. 2016, Wang et al. 2024).

2 - From what I understand, some parameterizations that are in section 2.1 are new compared to Wang et al. (2013), and should therefore be in section 2.2. : snow fraction, albedo

3 - On the implementation of the ice layer (section 2.2.2), why didn't you use ES directly to model this. On line 293, you say that snow density is limited to 750kg/m³. But if you had raised this limit to 900 or 950kg/m³, it seems to me that all the "snow" equations converge to "ice" equations, at least that's what comes out when we compare equation 26 to 20, 27 to 21, etc. In theory, if ES is done right, snow that has reached a certain density should be able to become ice. It would then be sufficient to initialize the height and density of the snowpack accordingly (e.g. the last 6 layers with an ice density and a total height of 100m for example). I don’t know if it’s possible but this could be discussed.

4 - Line 629 - 632: I understand here that the improved runoff modeling in OPT-12L would not be due to bias compensation. Well, I'm really not sure. What I understand from looking at your results is that to improve runoff compared to MAR, you need to set an albedo lower than MAR (Figure 6e), which inevitably induces a surface temperature (and surely an internal temperature of the snowpack) that is too high (Figure 7e) compared to MAR. I have the impression that this is also what Figure (8f) reveals. So, to obtain the same runoff than MAR, ORCHIDEE-ICE have to simulate a lower albedo than MAR to capture more energy. If it is true, it is perhaps due to the non representation of solar absorption by snow or a poor simulation of snowpack density.

5 - This last remark also underlines the fact that other important variables concerning the internal properties of the snowpack could be shown/analyzed, such as the temperature and density of the simulated snowpack compared with MAR, etc. This would enable a better understanding of the processes involved in the improvements related to one or another process claimed by the authors.

References :

Flanner, M. G., and C. S. Zender (2006), Linking snowpack microphysics and albedo evolution, J. Geophys. Res., 111, D12208, doi:10.1029/2005JD006834.

Schlegel, N.-J., Wiese, D. N., Larour, E. Y., Watkins, M. M., Box, J. E., Fettweis, X., and van den Broeke, M. R.: Application of GRACE to the assessment of model-based estimates of monthly Greenland Ice Sheet mass balance (2003–2012), The Cryosphere, 10, 1965–1989, https://doi.org/10.5194/tc-10-1965-2016, 2016.

Wang W., Yunzhong Shen, Qiujie Chen, Fengwei Wang, High-resolution mascon solutions reveal glacier-scale mass changes over the Greenland Ice Sheet from 2002 to 2022, Geophysical Journal International, Volume 236, Issue 1, January 2024, Pages 494–515, https://doi.org/10.1093/gji/ggad439