General comments

Bakketun et al. present the results of a new land surface reanalysis refining the snow scheme compared to the forcing CARRA reanalysis and focusing on the assimilation of in-situ snow depths with a new ensemble assimilation framework derivated from the well-known Ensemble Kalman Filter. The introduction presents clearly the context and objectives of this work and the production of an Arctic land surface reanalysis with improved snow modelling and snow data assimilation has certainly a great potential for various related scientific fields.

However, in general, the methods are not described with a sufficient accuracy level to allow their reproductibility (parameters used for perturbation, interpolation, observation errors, inflation, etc.). I believe the manuscript shoud be revised in that spirit. As a result, it is sometimes difficult to evaluate the relevance of the proposed methodology, see in particular my detailed comments relative to L149-151, L153-163 ; L214-217 ; L250. My understanding after reading the Methods section was that the main advantage of the EnKF was lost by the local assumption for application while the estimation of the snow model states are much more challenging to obtain than in the context of particle filters or particle batch smoothers. More explanations are necessary to understand this key methodological choices, because although I may have misunderstood something, at this stage, my feeling was that the chosen methodology might not be the most appropriate for this application. When reaching the results Section, I understand that LETKF is actually more complex than described, probably including neighbour pixels in the assimilation process but the readers should have fully understood that after reading the Methods even for readers not aware of all the variants of ensemble data assimilation algorithms.

Another critical aspect of this manuscript is the fact that the independence of observations is not considered as a mandatory property for evaluating the added value of data assimilation. Although some datasets are more or less independent, the results do not focus on this part. More critically, when considering only this dataset, I am afraid that the results are not really convincing to demonstrate the added value of their land surface reanalysis compared to the CARRA reference. As a result, I stay with a mitigated opinion where I see a really nice attempt to produce snow reanalyses with state-of-the-art methods, but leading to not really convincing results, questioning whether the chosen assimilation methodology was really appropriate, and with an unclear understanding of the behaviour of the assimilation algorithm due to unsufficient methodology description and maybe lack of examples on specific cases.

Therefore, I would recommend a major review for this manuscript to give the opportunity to authors to improve their description of methods, refine their evaluation process focusing on independent dataset, and maybe find new results emphasing better the added value of their methodology compared to their benchmark.

This being said, I am not fully convinced that assimilating only in situ snow depth observations has really a significant potential in mountainous areas with sparse observations with low spatial representativity, and I am really wondering why the authors do not prefer to start by assimilating snow cover fractions already assimilated in CARRA.

Detailed comments

L52 Note that this reference presents twin experiments (with synthetic observations). To the best of my knowledge, real observations of surface temperature have not yet been assimilated in a snow cover modelling system. I suggest to be more accurate on that point.

L54-56 The authors classify the DA algorithms according to their management of uncertainties but the algorithms also differ in their ability to preserve the multivariate consistency of model variables. It would also be useful to emphasize this point which is also a strength for instance of particle filters.

L59-60 Variable transformation is also a common practice to fit with the Gaussian assumption.

L61 Unclear at this step of the introduction which model refers to « the multi-layer snow model ».

L85-86 How was estimated this lookup table and is it variable in space ?

L87 and elsewhere : The display of units should be improved (space between kg and m^-2)

L88 I guess there is a typo in the unit (kg m^-2 and not kg m^-3). Can you explain the philosophy behind excluding observations « if the model exceeds 25 kg m^-2  » ?

L90-93 Does this mean that pseudo-observations and in-situ observations are considered to have the same observation error ?

L93-94 Are there other variables to describe the snow cover state in the model ? Are they kept at the same value during this process ?

L94-95 Can you explain the reason behind that ? I guess this is because surface temperature can not be updated ?

L96 « For modelling the land surface » → I understand that the idea is to produce a complementary offline simulation with an improved surface model compared to the CARRA coupled system, but this could be more explicitely explained in this transition.

L101 « the reduced heat flux between soil and snow » As a coauthor of Monteiro et al, 2024, I know what this is about, but I can imagine than most other readers would need more explanations.

L112 To encourage the reproductibility of the results, could you be more accurate on the spatial and temporal interpolation methods used for the different forcing variables ?

L129-130 Snow albedo and snow age do also need to be updated

L128-132 Direct insertion of snow depth has been applied in some studies with multilayer snow models by preserving the ratio between layers and other snow properties. The assumptions are strong in that case, but it could be mentioned to emphasize the added value of more advanced ensemble approaches.

L149-151 If EnKF is applied independently at each grid point, how can the increments be propagated in space ? This is rather counter-intuitive in the context of assimilation of in-situ snow depths in a spatialized system. Please provide explanations.

L153-163 Could you mention if vector x includes the 36 prognostic variables described line 129 ? If yes, how do you manage snow age and snow albedo ? What is the implication of computing error covariances with such a large vector ? What guarantees that the resulting analysed vector fills the discretization rules of the snow scheme ? Considering that for a given snow depth, the vertical discretization of ISBA-ES is fixed, I am not sure what is the interest to consider the 12 layer snow depths as independent variables in the state vector. I imagine it helps constraining unobserved density and heat contents. But could you justify this choice and provide more details to help the understanding beyond the mathematical formalism ?

Then, does vector x only include the prognostic variables of one single point ? If yes, how do you propagate information in space ? I see further from the results, that increments are spatialized but I do not see from this formalism how this happens.

L165 This assumption would mean that observations are as representative as the open patch than as the forest patch. However, in a large majority, observations are operated on open areas. Could you comment on that ?

L182-183 This is true but monovariate perturbations may produce unrealistic meteorological conditions with physical inconsistencies between meteorological variables, potentially leading to unrealistic model states.

L185 As mentioned before, the length of the control vector might be reduced thanks to the bijectivity between total snow depth and layer snow depths.

L187-189 The description of perturbations refers to several references which is not accurate enough for reproductibility. Please provide a table with perturbation methods and parameters (variance, auto-correlation time constants, etc.) for each variable.

L190-203 This is an interesting idea but even in Appendix A, the perturbation parameters to prescribe are unclear and not provided (I guess at least some horizontal distance defining the statistical properties of the random advection ?). Could you think again the description of the method as a basis allowing reproductibility ?

L212 « we initialized zero snow members with the ensemble-mean of members with snow. » It is unclear how is computed the mean. The different members have different layer thicknesses. Do you mean that the properties (density, heat) are averaged for each numerical layer regardless its depth in the snowpack ? Is it physically realistic ? Would not it be sufficient to extract the values of one single member with low snow depth rather than potentially mixing the properties of thin snowpacks with properties of thick snowpacks ?

L214-217 If the algorithm provides unrealistic states and crashes happen, I am wondering if applying thresholds is really the appropriate option compared to improvements in the assimilation formalism. To my mind, the main added value of the EnKF compared to particle filters is lost if it is applied independently at each pixel point. At the local scale and to assimilate only snow depths, particle filters would be best suited algorithms to retrieve all variable states without any possibility to obtain unrealistic model states, and they have already been largely applied within SURFEX even with more complex snow schemes. The added value of the proposed approach is unclear for me.

Table 1 : I do not think a threshold on temperature is necessary as snow heat, snow density and snow mass are sufficient to diagnose snow temperature. Or do you mean the threshold on temperature is applied to constrain the heat content ?

Section 2.6 Could you provide the elevation distribution of the different observation datasets ? Is a threshold applied on elevation as commonly done in NWP ?

I also understand that optical snow cover maps were not used here while they were assimilated in CARRA. This should be more explicit and justified as it could be questionable that an offline land surface reanalysis assimilates less observations than the original coupled reanalysis.

L231 How are initialized the simulations ? Do you apply any spinup to initialize the soil temperature profiles ? What is the assumption for snow depth ? I guess some pixels are affected by permanent snow ?

L242 Do you use independent observation datasets (not used in the assimilation process) to evaluate the system ?

L250 The fact that Figure 3 shows spatial patterns of increments means that something is definitely missing in Section 2.3 to produce the spatial propragation of increments when applying EnKF at the pixel scale. I guess this is associated with « the Gaussian localization function » mentioned L256 but this should be explained in Section 2.3 (providing methodological details as well as control parameters) as after my first reading I thought only one pixel was considered in vector x. The paper should be self-sufficient to understand everything without reading other papers describing LETKF. This raises numerous questions by the way about the localization radius. Is it constant or spatially / temporally variable ? Does it depend on the density of the observation network ? Etc.

L256 « The patterns of these increments are similar to the Gaussian localization function ». Without any knowledge about this localization function, I do not understand this comment.

L259-260 « The stations with mean negative innovations are situated in mountainous areas and have small impact on the analysis. » Why? Do you reduce some localization radius as a function of elevation ? Methodological details are missing to understand.

L261 « domain-averaged ». There are two different domains between Figure 2 and Figure 3. To which domain is applied this average ?

Figure 5 presents ensemble correlations between between layer variables and total snow depth. The results are useful to understand the assimilation process. However, it is not explained how this information is represented at a given depth in the snowpack as the different ensemble members have different snow depths and layer thicknesses. Do you choose a member ? Or compute median thicknesses ? Please provide the details in the text.

L267-270 The correlations between total snow depth and mass of each layer could be commented at the light of the vertical discretization rules of the snow scheme, explaining easily the weak correlation between the total snow depth and the mass of the surface and bottom layers.

L284-291 I would have expected that the effect of assimilation on density would have been described here (comparing densities obtained with CARRA-Land-Pv1 and CTRL), as done before with snow depth.

L292-296 The observations used to compute errors are not described. Are they independent from the dataset used for assimilation ? If not, I think it is unsufficient to evaluate a data assimilation algorithm with the assimilated data themselves. Then, the fact that assimilation deteriorates snow water equivalent is really problematic as the main advantage the authors emphasize in their method is its ability to attribute snow depth errors in both mass and density. Unfortunately, the results suggest that this idea fails. A much more simple algorithm assuming that the simulated snow density is correct may have been a better assumption ?

L300-301 If I understand well, OBS-ONLY-CARRA is the only relatively independent dataset to evaluate the added value of assimilation. So this is of course more challenging, but the description of results in the whole section should be mainly focused on this dataset, as imrproving snow depth at the assimilated observations is definitely not a proof that the system is valid at large scale. The word ‘validation’ is inappropriate in the title of the subsection (prefer evaluation). Unfortunately, the fact that the whole method deteriorates scores compared to CARRA (Fig 8c, value on the right of each column in Tab 3) is rather problematic while the goal of the study was to improve the simulation of snow cover with a dedicated offline reanalysis.

L307 It should be explained that CRPS is identical to MAE in the case of a deterministic forecast.

L333-334 I would like to be optimistic but compared to independent observations (OBS-ONLY-CARRA) not assimilated in the offline reanalysis, the results do not exhibit improvements of snow depths (Fig 8c, Table 3) and exhibit significant a deteriotation of density (Fig 7). I think a more qualified discussion of results is needed allowing a better questioning of the choices for data assimilation.

L343-357 I think the discussion about the link between density and mass could be improved.

« We emphasize that these stations are not collocated with assimilated snow depth observations ». In a spatialized assimilation system, I think the goal is to improve snow cover at large scale, not only at the assimilated observations.

L358 « By exploring the different observation sets used in CARRA and CARRA-Land-Pv1, we evaluate the analysis performance where observations are not available for assimilation ». Actually, only one dataset is independent and the results do not really focus on this one. I think the evaluation dataset used in this paper should be questioned and probably extended. A leave-one-out approach as in Cluzet et al., 2022 could be a valuable approach to consider all the snow depth datasets. Independent observations of satellite snow cover fractions could also be considered for a spatialized evaluation of the added value of data assimilation, especially in a context where such observations were assimilated by CARRA and are no longer assimilated in CARRA-Land-PV1.

L361 It is difficult to have a critical reading of this discussion due to the lacks in the methodology description. For instance, authors say that « the distance based localization used in CARRA-Land-Pv1 cause several stations in OBS-ONLY-CARRA to be unreachable in the analysis. C ». But the localization distance has not been defined neither provided, and the typical observation density is not known. I think a discussion should more rely on quantitative considerations provided to the reader.

L384-391 Authors explain again the theoretical added value of their framework compared to the simple assimilation process of CARRA. However, my feeling is that they have not found the results illustrating how considering variable background uncertainties is able to improve assimilation compared to their reference benchmark. If authors have illustrations on specific cases, as suggested, I would encourage to present the material in the paper to better emphasize the potential added value, if it can not be demonstrated with finalized scores.

L393 implemented

L396-397 Again, this conclusion seems rather optimistic considering only independent evaluation datasets.

L399-401 Again, my feeling is that this theorical advantage is not completely supported by results.

My comments are in the attached pdf file