Padrón et al report an analysis of drivers of the terrestrial carbon sink in the CMIP6 ensemble scenario SSP126 where warming is limited to 2 oC. This is a useful study of the latest CMIP model results in a policy relevant scenario, showing that terrestrial carbon sink projections by 2100 (cumulative NBP) in the ensemble vary from 56 to 207 Pg C, mean 144 and standard deviation 47 Pg C. Using linear regression Padrón et al partition this variability among sensitivity to CO2, temperature (T), soil moisture (SM), and differences in baseline temperature and soil moisture. Their methods show that the greatest proportion of this variance is explained by sensitivity to T and SM combined, with sensitivity to CO2 as the second most important driver of variability. Based on these results, they conclude that the gamma feedback (climate) is greater than the beta feedback (physiological) under this policy-relevant scenario and thus climate sensitivities require the greatest attention. They also show compensating drivers of cumulative NBP variability such that reduction of uncertainty in response to one driver would not greatly reduce overall NBP variability.

Overall this is a well written and executed study. The analysis of the relatively low-warming SSP126 scenario is timely and to my knowledge has not been done before. I have several comments and criticisms that I hope will help to make the analysis and conclusions more robust and impactful. First, I think that for a number of reasons the method has low-biased the estimation of the impact of CO2 sensitivity on NBP variability. Second, I encourage a little more quantification and thought into exactly what is quantified and communicated.

First: underestimation of the impact of CO2 sensitivity on NBP variability. CO2 sensitivity is estimated as the sensitivity of GPP to CO2 in the 1 % per year increases in CO2 simulations (1pctCO2-bgc) in which CO2 ranges from 350 to 800 ppm. This method assumes 1) the CO2 sensitivity of NBP is the same as that for GPP, 2) that CO2 sensitivity is linear across the range 350 to 800 ppm, and 3) that there are no interactions between CO2 sensitivity and either T or SM. It is likely that all three of these assumptions will low-bias the estimate of the impact of model CO2 sensitivity on cross-model NBP variability.

• While GPP sensitivity to CO2 is likely the main driver of NBP sensitivity to CO2, as asserted in the current ms, the assumption ignores potential changes in turnover rates that can also occur in response to CO2, which can be substantial. Using cross-model GPP sensitivity to CO2 will result in a lower correlation with NBP variability than using NBP sensitivity to CO2. Further, for T and SM sensitivity, NBP is used, biasing results in favor of T and SM sensitivity. Comparing the sensitivities of GPP to CO2 to NBP to T and SM is not a like-for-like comparison. Sensitivity of NBP to CO2 should be estimated and used in the regression analysis.
• The CO2 response over 350 to 800 ppm is likely not linear in these models, it almost certainly is not at the leaf scale which drives model CO2 responses. The SPP126 simulations max out at 446 ppm. There is likely saturation in the CO2 response for many models somewhere between 450 and 800 ppm. CO2 sensitivities should be estimated over the range of CO2 concentrations that preserve linearity over the range 350 – 446 (i.e. concentrations can be higher but responses must be linear over the range). A supplemental figure showing NBP against CO2 for the 1pctCO2-bgc simulations would be useful.
• There are interactions between CO2 and T and SM. Interactions with T are likely the most important for this discussion. At high T it is well known that CO2 can alleviate some of the reductions in photosynthesis due to interactive effects on photo-respiration. This could alleviate GPP reductions in high T years that I’m not sure would be removed by detrending NBP. I’m not sure there is an easy way to account for this, and that is OK. But some acknowledgment of this effect and some attempt to quantify it would help make results more robust.

Second: encourage more thought into exactly what is quantified and communicated. I suggest quantifying statement in the abstract, cumulative NBP variability etc. Also, as well as putting these numbers in the context of current annual emissions, I think it might also be useful to present them as a proportion of the assumed emissions in the SSP126 scenario (if someone has calculated those). Why is the proportion of variance in NBP variability to CO2 sensitivity not quantified on ln 354? I encourage the authors to think about what is best to present given this is a study of the global carbon cycle. Most figures are presented in the units per meter squared. When aggregating to broad zonal regions I suggest it is more informative to present results as the absolute sum across the whole area – this would make it easier to relate the regions and sensitivities directly to the global aggregate numbers. Finally, have differences in grid-square area been taken into account when presenting the global aggregate drivers of NBP variability?  

Technical Comments

Title: suggest switching “controls” for “drivers” as control suggests some degree of intention.

While I see some of the benefits of the narrative style with the methods spread throughout the results (e.g. lns 157-164, 179-189, 263-283, etc), I think it is more practical to have the methods all in one place where they can be found easily and assessed side by side.

Fig 6: can probably go into the supplement.

Fig 7: A little hard to read, I think the sensitivities could be presented more clearly if they were presented as in Fig 8. I recognise that would necessitate removal of uncertainties from the figure and I appreciate the effort made to quantify uncertainty but is clarity of communication is the trade off. Fig 8 as is could go to the supplement.

Fig 8: Suggest adding a dot for the actual cumulative NBP. Also I really think this would be better off expressed in global sums rather than per meter squared. The wite dot could be a little larger.

Ln 30: There are several commentaries explaining why Wang et al 2020 is not a reliable analysis.

Ln 37: Are there other disturbances that release C directly back to the atmosphere?

Ln 61: Note the editor’s note for Keenan 2021

Ln 69: can delete “consider it important to instead”

This manuscript takes 11 ESMs from the CMIP6 archive and attempts to unravel the root causes of the uncertainty in the NBP simulated by each model over the (roughly) 2 deg warming scenario. The models' sensitivity to CO2 (through the 1 percent runs) and to temperature and soil moisture are investigated for both short and long timescales. I found the paper figures to be generally well designed (though see my comment about Fig 3) but the text could be confusing at times. There is a lot of rather convoluted steps/arguments in producing the T/SM/sT/sSM metrics and it sometimes was hard to understand exactly what they were telling me about the models. I think the paper is publishable, but needs revisions for clarity. An obvious target for clarity/context would be to discuss the results of this work in the context of previous efforts as discussed in the introduction (principally Arora et al. 2020). I found myself comparing this work to that paper and not understanding why they differed strongly in some cases.

Main comments:

Fig 3: I found this to be a strange figure. So if a model has a positive correlation it counts towards the blue end of the colour scheme, whereby if it is negative it counts towards the red. However this seems to have no consideration of how positive or negative a model was. I think it would treat a model that is 0.9 the same as one that is 0.0009, which seems to be a bit too ambiguous. Also what if all 11 models are +0.0001 vs. all models are >0.9, as is they would appear the same in this figure but arguably the situation where all 11 are >0.9 is more interesting than the situation where the models are rather ambiguous (close to 0). I would suggest reconsidering this figure.

Deserts - how did you mask the deserts? I think the grid cell sizes of the ESM precludes removing many deserts, e.g. the Atacama. Instead it seems like only a few were removed (Sahara, around Middle East, and Gobi) but I am not sure why those made the cut but not, for example the deserts of western Australia or the US SW. What impact does it have keeping them in? Greenland makes sense since there is no vegetation at all but some of the world's dry regions have been fingered as influential in the global C cycle (e.g. Ahlström et al. 2015), so exactly where masking applies could have impact I would assume.

Fire - Fire is mentioned on line 335 but ignored otherwise, why? I see you mention which models do fire in Table S1.

Smaller comments:

Line 44: 'with  drought-related observed decreasing trends in leaf area' consider rewording, confusing.

Line 103: CanESM has an implict N cycle (empirical downregulation scheme see Arora and Scinocca 2016)

L 160: This explanation is confusing. Perhaps spell it out in a bit more detail.

L 182: Why 22.5 degrees and not 25 or 30 or some other number? It just seems awfully precise for a seemingly arbitrary limit.

L 235: remind reader that both use CLM?

L 290: 'underestimation  of  the  land  carbon  sink  modelled  by  NorESM2-LM and  CanESM5,' where is this shown? I can't seem to see any figure where CanESM5 sticks out with an underestimation of the land C sink but it is mentioned here and line 357, indeed in Figure 1 it seems to have one of the highest cumulative NBP. What am I missing?

CanESM5: Other papers (Arora et al. 2020) have suggested that CanESM5 has the largest land C uptake (at least for the 4X CO2 simulations) so it is surprising that it is suggested to be underestimated for the land C sink. Can you clarify how the same model appears to be on the low/high end depending on the analysis? I realize these are different scenarios but I would have assumed high CO2 sensitivity would follow in both (but be exaggerated in the 4XCO2 run), but I don't see high CO2 sens in Fig 7. I assume I missed something here as you mention the Arora et al. paper in the intro but don't return to place your results in context of those other works.

L 352: And assumedly many of them use Nemo for their ocean so model commonalities are not just atm/land. Never mind all who use Farquar photosynthesis etc. 

L 378: 'Outperfom' seems out of place, consider swapping it out with something like 'be more important than'

Lit cited:

Ahlström, A., Raupach, M. R., Schurgers, G., Smith, B., Arneth, A., Jung, M., Reichstein, M., Canadell, J. G., Friedlingstein, P., Jain, A. K., Kato, E., Poulter, B., Sitch, S., Stocker, B. D., Viovy, N., Wang, Y. P., Wiltshire, A., Zaehle, S., and Zeng, N.: The dominant role of semi-arid ecosystems in the trend and variability of the land CO2 sink, Science, 348, 895–899, 2015.

Arora, V. K., Katavouta, A., Williams, R. G., Jones, C. D., Brovkin, V., Friedlingstein, P., Schwinger, J., Bopp, L., Boucher, O., Cadule, P., Chamberlain, M. A., Christian, J. R., Delire, C., Fisher, R. A., Hajima, T., Ilyina, T., Joetzjer, E., Kawamiya, M., Koven, C. D., Krasting, J. P., Law, R. M., Lawrence, D. M., Lenton, A., Lindsay, K., Pongratz, J., Raddatz, T., Séférian, R., Tachiiri, K., Tjiputra, J. F., Wiltshire, A., Wu, T., and Ziehn, T.: Carbon–concentration and carbon–climate feedbacks in CMIP6 models and their comparison to CMIP5 models, Biogeosciences, 17, 4173–4222, 2020.

Arora, V. K. and Scinocca, J. F.: On constraining the strength of the terrestrial CO2 fertilization effect in an Earth system model, https://doi.org/10.5194/gmd-2015-252, 2016.