Natel and coauthors are interested in developing emulators to allow easier integration of ecosystem models (like LPJ-GUESS) in broader frameworks that couple multiple models together. LPJ-GUESS is sufficiently computationally expensive that an emulator could be valuable for the multimodel frameworks (in particular LandSyMM). They use two machine learning based approaches: random forests (RF) and neural network (NN). Both emulators were trained using LPJ-GUESS outputs for some historical and future simulations. The emulators differed in both their performance and the main variables they were sensitive to but both were much faster than LPJ-GUESS itself.

The paper is generally well written and easy to follow. The work falls well within GMD's area of interest. I think the work is suitable for publication but have several questions that I would like to see answered beforehand.

Main comments:

1. I liked that the authors used two different ML-based approaches in their emulators and then attempted to understand/interpret what each emulator was sensitive to. This is valuable information but it felt like only half the story. What was missing was what LPJ-GUESS is sensitive to. If the point of the emulators was to allow cheaper approximation of LPJ-GUESS (the 'model') then the most important thing is that the emulator is responding in the same manner and to the same variables as LPJ-GUESS. There are many plots showing trajectories of pools and fluxes for LPJ-GUESS and the emulators (e.g. Fig 2 and 4) but no similar plots showing sensitivity of LPJ-GUESS as there is of the emulators (e.g. Fig 6). I realize this is more challenging with LPJ-GUESS since it is inherently a completely different kind of model, but I struggle to understand how one can trust either emulator without knowing if it is actually mimicking the model's sensitivities (which under this circumstance has to be assumed to be perfect).

2. I am quite skeptical of the claim that (L 342) 'NNs excel at modeling continuous relationships, making them more capable of generalizing to unseen data, particularly when extrapolation is required'. My read of Muckley et al. (2023) does not support this contention. Muckley et al. test out the performance of linear regressions and black box RF and NN models. For the interpolation tasks, the linear model was poor but when it came to extrapolation it could out perform the black box models in some of the tests (~40%). The authors state the '{linear regressions}... may be desirable over complex algorithms in many extrapolation problems because of their superior interpretability...'. Lakshminarayanan et al. (2017) nicely demonstrate this for a toy example using a NN whereby the NN extrapolates poorly (their Fig 1 - left panel is bounds of 5 NNs). They show that the uncertainty bound via an ensemble technique can be created that encompasses the true function. So, getting to my main concern, given that NNs do not extrapolate well (same with decision tree-based methods), how can we trust the NN/RF models when they are forced to extrapolate? The approach here doesn't have any way to uncertainty bound the emulator results so it can extrapolate (poorly) blindly to the user. I don't expect the authors to fix this problem right now, but I would like to see more discussion about this difficulty and how it could be addressed for emulators as their use if becoming more common.

Minor Comments:

Supplement - Can you explain more about the disturbance interval of 100 years and how that is applied? Also with fire off, the disturbance is then what? I see land use change is also not used.

L 84 - This sentence is a bit confusing as I found it less clear what the features were selected for.

L 112 - I wonder about the influence of this post-processing step for non-negative C stocks. How often did this come up? Were the instances where this came up for regions with very low stocks such that the inaccuracy would be small? (e.g. true value is 0.1 kg C/m2, so going negative is fairly reasonable but if it was really supposed to be 10 kg C/m2 then that is a big problem). This also demonstrates a problem with using an off-the-shelf NN whereby it has no knowledge of boundaries that one like a physics-informed ML model could.

L 137 - sampled by grid cell, time, or ?

L 153 - with the meteorology/climate of MPI..., not the actual climate model itself.

L 155 - Does LPJ-GUESS not have a dependence for computational cost on the number of PFTs present or the soil permeable depth?

L 158 - If it takes 5000 sims to train the emulator but actually running the model that uses the emulator only happens 1000 times then you may end up with no net benefit. Also, sorry if I missed it, how many simulations did you need to train the emulator? 

L 207 - Could you give the real values in addition to the percent. It would be nice to see how much in clock time these cost (acknowledging it is system dependent).

L 215 - I think the lack of decline in performance is simply due to training with the most extreme ends of the scenarios. This ensured that you were interpolating as much as possible. This is likely the only reasonable approach given the these techniques do not extrapolate well (see one of my main comments). But it means that the emulator always requires retraining for new scenarios and the scenarios always need be more extreme than what the actual system should realistically experience. I think some aspects of this bear mentioning.

L 218 - 'greatest accuracy' - by NN? Unclear as written.

Fig 2 and 4 - What about adding new plots presenting these and the fluxes as cumulative plots so the impact of over/under predicting over time are visible? The fluxes is important as it has impact on how much C the land surface takes up/releases. The stocks as it changes how much C is emitted during disturbance or land use change. A cumulative plot can show the effect across the simulated period.

Lakshminarayanan, B., Pritzel, A., and Blundell, C.: Simple and scalable predictive uncertainty estimation using deep ensembles, arXiv [stat.ML],31st Conference on Neural Information Processing Systems (NIPS 2017), Long Beach, CA, USA. https://arxiv.org/pdf/1612.01474

Carolina Natel et al. developed both a Random Forest model and a neural network model to emulate the dynamics of ecosystem carbon fluxes and carbon pool changes. These machine learning models have been widely applied to land models or components in the past and have consistently demonstrated effectiveness. Similarly, this study shows reasonable performance in emulating the target variables. The paper is well-written, with well-documented data and code. Overall, this is a solid modeling paper. Below, I have a few specific comments:

1. Interpretability vs. Physical Consistency

My primary concern is the balance between model interpretability and physical consistency. While SHAP has been used to interpret the ML models, it does not ensure that the emulators capture established physical knowledge embedded within the land model (in this case, LPJ-GUESS). I encourage the authors to further explore this aspect, as it is fundamentally important to understand the functional relationship emerging from ML emulators.

For example, in land models, atmospheric CO₂ concentration is a key driver of vegetation productivity, while temperature (T) strongly influences soil carbon stocks. Ideally, such first-order relationships should also be reflected in the trained ML emulators. One way to test this would be to leverage factorial LPJ-GUESS simulations with:

(a) Future SSP climate scenarios +  historical atmospheric CO2 levels.

(b) Future SSP CO2 levels + historical climate conditions (e.g., repeated climate from 2010–2020).

If the trained ML emulators can reproduce the results of these factorial runs, it would provide strong evidence that the emulators have captured critical relationships between environmental drivers and carbon dynamics.

2. Justification for Annual Time Step

Further justification is needed regarding the choice of an annual time step. Land models typically operate at much finer temporal resolutions (e.g., daily, hourly, or at least monthly). It would be helpful to explain why the annual scale was selected and how potential loss of information at shorter timescales may affect the emulator's performance.

3. Capturing Inter-Annual Variability

Given the focus on annual time steps, evaluating the emulator’s ability to capture inter-annual variability in carbon fluxes (in addition to long-term trends) would be an important validation metric. Although the training is performed at the grid cell level (random sampling), it may also be valuable to include spatially aggregated fluxes (e.g., global/regional totals) as part of the loss function. This could improve the model's ability to represent inter-annual variability at regional or global scales.

4. Treatment of Disturbance Intensity

The results show disturbance as one of the most important features. However, it remains unclear how disturbance intensity (e.g., fractional area burned by wildfire, or land-use/land-cover change) is handled. How does the ML model represent partially disturbed grid cells? Additional clarification on this point would be needed.

5. Recommendations for future work: Since land models simulate continuous, time-dependent changes in carbon fluxes and pools, it may be worthwhile to explore time-series ML models (e.g., RNNs, LSTMs, or Transformers) in future work. Such models could potentially outperform static models like Random Forests and ANNs by better capturing temporal dependencies and dynamics.