Evaluating the overall quality ("general comments"), 

The authors conduct a twin model experiment to test the ability of two competing parameter estimation techniques (eta-4DVAR vs. 4DEnVar) to constrain 54 parameters within the land model ORCHIDEE.  The authors use synthetic observations to generate atmospheric CO2 station data in which to retrieve the ‘true’ parameters from prior parameter distributions.   The authors find that the 4DEnVar approach performs the best in terms of the RMSE statistic of global NBP and in terms of the posterior parameter values relative to the true values.   The authors claim this demonstrates strong potential for 4DEnVar to be used with real data, and should be widely applicable to other land surface models.

This reviewer found the topic relevant to the current state of earth system science which requires a wide range of observations to calibrate model performance to improve forecasts/projections.   The potential for such a problem to be ill-constrained and suffer from equifinality stood out to this reviewer given the use of a single data set (atmospheric CO2) to constrain a multi-dimensional problem.  This author recommends a more nuanced discussion of equifinality for this application both in the twin model experiment and for potential applications of using real data (see scientific questions below).  This reviewer also would have appreciated a better description of how the parameter values were sampled (perturbed) from their respective distributions – and also to what extent the authors could have presented their results for the 4DEnVar using parameter distributions rather than point estimates.  It was somewhat surprising how well the posterior parameter values improved the global distribution of NBP, given the limited range of the CO2 station footprint.  More discussion related to the distribution of land PFTs relative to the CO2 station footprint may have been helpful to aid this discussion.

Individual scientific questions/issues ("specific comments") 

Lines 45-50:   If you calibrate against CO2 observations only, and do not adjust for model biases in land carbon pools – how accurate can your projections/forecasts of the carbon cycle be?

Line 95:  “We demonstrate the potential of 4DEnVar using synthetic observational data and compare its performance with that of 4DVar with finite differences.”

Your criteria for testing the differences between the methods is vaguely stated here.  Are you judging success based on which method best identifies the ‘true parameters’.  

Lines 135:140:   Given the use of pre-calculated transport fields that relate atmospheric concentrations to surface fluxes,  you do not make use of a dynamic atmospheric ensemble to generate this relationship based upon actual atmospheric forcing.  This seems a bit like using a background climatology to get the concentration/surface flux relation.    Also the land seems to be decoupled from your atmosphere – in the sense the dynamics that drive the CO2 concentration/flux relationship is not the same as the actual met forcing driving the land model.   This is perhaps not as important given the authors are generating ‘perfect’ obs, but during implementation for real parameter estimation should this not have an important impact?

Line 150-155, Figure 1:  One would expect that the CO2 sites were also chosen such that land surface areas sensitive to atmospheric CO2 also coincide with the range of land surface PFT types in this analysis.  Any consideration of this?

Line 280-85 – Same question as before, these sites were chosen to be sensitive to land surface fluxes, however, do they provide good sampling of the most important PFTs?  Sampling of North America and South America look poor.   Sampling of Africa does not include the tropics at all…..

Line 286:  How was your prior parameter distributions chosen?  From the uniform distribution shown in the figures where you only show upper and lower bounds?   Or from a normal distribution as described in Equation 14 and 15?

Figure 3:  Showing the ensemble mean parameter behavior  doesn’t give any information on the ensemble distribution.  Maybe I am misinterpreting the implementation of the 4DEnVar method, but can’t you show this in terms of the true, prior and posterior *distributions* instead of the ensemble mean behavior?  

Figure 6:  Same question as before – can you convey this information in terms of distributions (histograms)?

Figure 6: (SLA panel) Why do most of the prior values for the parameters all start at the same value?  Were they not being perturbed independently?

Table A1:   What does proportion mean in this context?

Line 430:   As far as I can tell, it is still not defined what the ‘proportion’ of the parameters is.  Is it based on global land area coverage, or land coverage that coincides with spatial footprints from your chosen network?  These could be two completely different things.  Did you do any comparison of the MAD statistic based on % of land area covered by the CO2 network spatial footprint?  Are they strongly related?  Would be nice to see a land surface map with PFT distribution.

Table A2:    Is the partial derivative  averaged over space and time?  Therefore the 4dEnvar itself doesn’t account for any seasonality (time-variation) in the relationship?  

Figure 7:  I was surprised at how well the  True – Posterior 4dEnvar net carbon flux (top right panel) performed given the limitation of the spatial footprint influencing the station CO2, thus informing the biogenic contribution to CO2.    This is promising, but be aware,  that the ability to match the net carbon flux gives no guarantee that the component fluxes are well simulated.   An interesting complement to this plot would be to compare the true component fluxes of GPP and ecosystem respiration against the posterior component fluxes of GPP and ecosystem respiration for the complex case.

Line 431:  I think a more nuanced discussion of equifinality is required here and/or in the Discussion.  In addition to the pure number of parameters attempted to calibrate simultaneously – equifinality can arise for a number of different reasons – 1) a single parameter type being compensated within the large list of PFTs, 2) the station CO2 concentration is influenced through the NBP, which is a confluence of both photosynthetic and respiration processes, which can easily compensate for each other to provide a net biogenic carbon flux consistent with station CO2 data.   I understand that this is a twin model experiment, so the following do not necessarily contribute here, but if this setup were to be applied to real data additional equifinality challenges present themselves including 1) the model state itself (carbon, water nutrient pools) have not been constrained by any data, thus parameters will compensate for biases due to model state initialization problems 2) The biogenic fluxes (controlled by parameters) would seem to contribute just a portion of the land-atmosphere carbon exchange which includes other large fluxes from fossil fuel, fires and ocean fluxes which would have to be measured accurately--  3) atmospheric model transport errors, influencing the relationship between land carbon flux and station CO2 data.  

Line 512:  Given the significant challenges related to equifinality mentioned above, I am not sure this setup shows “great potential” to constrain parameters.   I might be more realistic and state that it demonstrates that 4DEnvar shows more potential than eta-4DVar.

Purely technical corrections 

Abstract:

Awkward:  “These models rely on parameterisations that necessitate to be carefully calibrated”.    These models rely on parameterizations that require careful calibration.

Introduction:

Can you describe in terms more accessible to general community what isotropic means in this context?

“corrections to CO2 surface fluxes are isotropic in time and space.”

Line 115:   I wouldn’t use the terminology ‘assimilation’ routine to describe photosynthesis or carbon uptake routine.  Assimilation is often used within data assimilation context, a component of this analysis, which is not what this is describing.

Section 2.1.4.  I think it’s also worth mentioning that you didn’t optimize the prior biogenic fluxes by constraining them with carbon pool observations (LAI, biomass, soil carbon etc).

Shouldn’t Figure A1 include the land biogenic fluxes for the truth simulation, just for relative perspective?  After all the parameter optimization is based on the influence of biogenic fluxes on the atmospheric CO2.

Line 305:  LAImax:   The absolute max value that LAI can be? Can you clarify?  Does this mean carbon cannot be allocated to leaf carbon once achieving this level?

Line 444:  LAI or LAImax ?

Summary of manuscript

The authors explore the assimilation of atmospheric CO₂ concentration data for parameter calibration of the ORCHIDEE model using the 4DEnVar data assimilation method. Through carefully designed data assimilation experiments, they demonstrate the capability of 4DEnVar in assimilating atmospheric CO₂ concentration for parameter calibration, and highlighted its superiority over the ε-4DVar method in terms of computational efficiency, parameter recovery, and fitting to CO₂ concentration observations.

General comments

Due to the continuous evolution of land surface models (LSMs), the use of 4DVar and gradient descent method for calibrating LSM parameters faces significant challenges in maintaining the tangent linear and adjoint models. Thus, it is necessary to explore adjoint-free variational methods. However, as the authors mentioned, the results obtained here for the ε-4DVar are not equivalent to a standard 4DVar, and no conclusions can be drawn regarding the comparison between the 4DEnVar and standard 4DVar methods. In light of this, to what extent can this study provide insights and practical guidance for the application of 4DEnVar to other LSMs and the assimilation of real, multi-source observations? I believe the manuscript would benefit from a clearer articulation of its research significance.

The manuscript focuses on the introduction, application, and evaluation of the 4DEnVar and ε-4DVar methods throughout the methodology, results, and discussion sections. However, this focus is not well reflected in the title. Perhaps the authors could consider revising the title in light of related works, such as Yaremchuk et al. (2016).

In the comparison of the assimilation results between the 4DEnVar and ε-4DVar methods, the authors repeatedly attribute the poorer performance of the ε-4DVar method to the fact that it falls into a local minimum. However, for the 4DVar method, whether the parameter iteration converges to a local minimum undoubtedly depends on factors such as the a priori parameter vector. This study employed only one a priori parameter vector, and its generation process was not clarified. This raises concerns about the reproducibility and generalizability of the findings. In other words, would different conclusions be reached if a different a priori parameter vector was used?

A more detailed description and presentation of the methods and results are needed. The manuscript currently lacks an explanation of the parameter set's value range and sampling approach. It would be beneficial to include formulas that demonstrate how the selected parameters influence ecosystem processes such as photosynthesis, respiration, and other carbon cycle components. Personally, I would appreciate seeing the distribution of the parameter ensemble and the spread of the ensemble simulations, as presented in Pinnington et al. (2020).

The authors may need to consider citing and discussing some recent studies, such as Douglas et al. (2025).

Specific comments

Line 34: “i.e.” to “i.e.,”.

Lines 50-51: Pay attention to the spacing before or after the paragraph.

Line 54: “4DVar” to “four-dimensional variational (4DVar)”. Please check the use of abbreviations in the manuscript to ensure they are correct.

Lines 56-58: It is necessary to add references here, such as Talagrand and Courtier (1987).

Line 83: The citation format is incorrect and needs to be changed from “Pinnington et al. (2020)” to “(Pinnington et al., 2020)”.

Line 92: The space between 'approaches' and ',' is extra.

Line 92-94: The sentence is not concise and clear. It is recommended to revise it as follows: “Although tangent linear or adjoint models are not required for methods such as GA, MCMC, or emulator-based approaches, these methods necessitate defining a large ensemble, making them unfeasible for use in this study due to the time-consuming nature of model simulations.”

Line 123: The period currently at the beginning of the line should be placed at the end of the previous line.

Line 151 and Figure 1: You mentioned the stations are selected according to their 6-month averaged sensitivity. Which six months were chosen? Given seasonal variations, it would seem more reasonable to select a full year or multiple years. Additionally, you may need to clarify whether any climate pattern, such as ENSO or IOD, occurred during the sensitivity analysis period and the simulation period. Please provide a more detailed description.

Line 160: The version of the Global Fire Emission Database used in the study is outdated, or why used this one?

Line 171: Please verify that the equations are correctly written. For example, vectors should be in italics, while matrices should not.

Line 278: It is suggested to consider organizing the default parameter values in a table and placing them in the supplement.

Line 284: It is necessary to specify how the a priori parameter vector was obtained.

Lines 289-290: “Vcmax” and “℃” should not be italicized.

Line 312: It is recommended to provide some references regarding this setup.

Lines 313-314: A more detailed explanation of the parameter range settings and the rationale behind them is needed.

Lines 364-365: Although RMSD and MAD are common statistical metrics, I still recommend that the authors provide their calculation formulas and explanations here. Since the observations have already been synthesized, are the simulation results involved in the calculation also synthesized? Furthermore, both RMSD and MAD, in terms of their form, resemble the cost function, as they include the critical term representing the difference between observations and simulations. In assimilation experiments, reductions in these metrics are expected. It would be valuable to explore additional metrics with distinct physical interpretations (e.g., coefficient of determination, R²) to comprehensively assess method performance.

Line 367: There should be a space between the number and the unit.

Line 370-372: I don't fully understand why configurations with more ensemble members (e.g., 350, 400) result in a smaller RMSD reduction. Could the authors provide an explanation?

Line 395: The use of 'seem to' here makes the experiment appear insufficiently rigorous.

Line 412: Use exponential notation and change “GtC/year” to “Gt C year⁻¹”.

Line 523: The line break in the link seems to be problematic.

Figure 1. The website should include the date of the last access.

Figure 3. “triangle” to “triangles”.

Figure 4: It is recommended to consistently retain two decimal places.

Figure 7: The presentation of the last subplot can be improved, for example, the current color scheme does not match well with that of the other subplots.

Figure A1: Revise unnecessary capitalization and add a comma at the end of the sentence.

References

Douglas, N., Quaife, T., and Bannister, R.: Exploring a hybrid ensemble–variational data assimilation technique (4DEnVar) with a simple ecosystem carbon model, Environmental Modelling & Software, 186, 106361, https://doi.org/10.1016/j.envsoft.2025.106361, 2025.

Pinnington, E., Quaife, T., Lawless, A., Williams, K., Arkebauer, T., and Scoby, D.: The Land Variational Ensemble Data Assimilation Framework: LAVENDAR v1.0.0, Geosci. Model Dev., 13, 55-69, 10.5194/gmd-13-55-2020, 2020.

Talagrand, O. and Courtier, P.: Variational Assimilation of Meteorological Observations With the Adjoint Vorticity Equation. I: Theory, Quarterly Journal of the Royal Meteorological Society, 113, 1311-1328, https://doi.org/10.1002/qj.49711347812, 1987.

Yaremchuk, M., Martin, P., Koch, A., and Beattie, C.: Comparison of the adjoint and adjoint-free 4dVar assimilation of the hydrographic and velocity observations in the Adriatic Sea, Ocean Modelling, 97, 129-140, https://doi.org/10.1016/j.ocemod.2015.10.010, 2016.