This is a straightforward paper evaluating the inclusion of P cycle processes in JULES. The work is obviously important for the importance of P in estimating the global C sink and the important role of JULES in the modelling community. I have several comments to hopefully help further strengthen the manuscript.

1. Given that the site is a nutrient fertilization experiment, it seems to be a missed opportunity to not evaluate the model performance in response to nutrient fertilization. The evaluation of the CO2 response is obviously still important, but I don’t see enough rationale to do it given that there is limited data available to actually evaluate the simulated CO2 responses. Can we learn more by simulating the nutrient fertilization response where there is actual data?
2. Some details of the model description is not available. For example, how P interacts with N to affect allocation, plant growth, stoichiometry and nutrient uptake?
3. Can the authors spend some efforts highlighting the novelty of this work? I understand that adding P cycle into JULES is a significant work, but some readers could take it as “yet another model with P cycle”. Considering the spread of the simulated CO2 responses among CNP models in Fleischer et al. 2019, having yet another model doesn’t necessarily reduce the model spread and our knowledge gap. Is there anything specific about having a P cycle in JULES would help to resolve? In other words, what unique features JULES-CNP can provide relative to other models?

Specific comments:

L31: Unclear – what does this 60% mean? Is it the 60^th percentile of P availability (doesn’t sound like a particularly P limited site)? Or does it mean the total P/available P at the site represents 60% of soil across the Amazon? Be good to make it clearer.

L32: What is the eCO2 treatment? Show a number is useful.

L32: The model is able to reproduce observed plant and soil P pools under ambient conditions – to what extent are these values provided as parameters and/or targeting values to tune the model?

Abstract: Given that you evaluated the model at a nutrient fertilization site, does it make sense to apply the nutrient treatment to the model and try to evaluate the model against observation? It’s quite unclear what to learn from evaluating the model against a CO2 treatment where you don’t have any data. The rationale about the CO2 effect isn’t very clear in the abstract.

L42: net primary production.

L43: C, not carbon.

L47: It’s not either or. There are other nutrients involved.

L99: Hou et al. 2019 – are you sure this is a DGVM?

Figure 1. Do you have P leaching flux and biochemical P mineralization flux? What about atmospheric P deposition?

From the equation it seems that you would have a continuous build-up of P in the system, because there is no way out from the occluded P pool.

L239: This exudate term is confusing. Can I consider it as part of autotrophic respiration? It’s a missed opportunity where this exudate term isn’t pumped into soil to facilitate plant-soil interaction, which could potentially be used to alleviate plant nutrient stress. See for example (Jiang et al., 2020).

L284: What timestep does retranslocation occur?

Table 2: Would be good to show P retranslocation coefficient.

Method: Description on N-P interaction is missing. I suppose P is not entirely independent from N processes? I think it should be useful to describe how P affects N, and vice versa in the model.

L404: What about N only model? What about Jmax?

L428: I may have missed this, but does the model assume fixed or variable CP ratios?

L450: What does the symbol -//- mean?

L469: What does a step increase mean?

L664: What is SMCL?

L767: Jiang et al., not Jing :)

This study implements phosphorus cycle processes into the JULES land surface model. The authors calibrate and test their model using observation from the Amazon Fertilization Experiment site and perform a theoretical elevated OC2 experiment. The inclusion of the P cycle in models is critical to our ability to predict ecosystem responses to elevated CO2 and climate change in low-P systems such as the tropics and subtropics, and the majority of land surface models are now adding these processes. 

My main concern about this paper is not with the model development per se, which is largely robust, but with the very extensive calibration. The authors calibrate their model very intensively using data from a very well instrumented site, which leaves me wondering if JULES CNP can be applied to any other site. In particular, the leaf C:P ratio, which the authors themselves show is a very sensitive parameter (Fig. 4) can be highly spatially variable and I do not see how the model can be run at sites where this data is not available, not to say anything about globally. Ideally, I would like to see the model validated at a site at which it has not been calibrated, but I understand this can be very difficult as it would involve obtaining more data. Alternatively, the authors could perform a more comprehensive parameter sensitivity analysis and include a discussion of the generality of their model.

Additionally, I think the paper would benefit from a more extensive discussion of the implementation of exudates and the knock on effect on CUE. The question of what to do with excess carbon under nutrient limitation is one that all models face and there have been a variety of solutions: respire it (what JULES CNP does too), down-regulate photosynthesis, decrease tissue nutrient content etc. I am not saying that the choice made here is necessarily wrong , but it does have implications for the model results. In particular, changes in carbon use efficiency are a direct result of this modeling choice and do not necessarily have an interpretable meaning.

Detailed comments:

L 241 As far as I understand from the description here, what the authors term ‘exudates’ is just excess C that is respired by the plants. However, the term normally refers to carbohydrates released by plants into the soil for a potential benefit in additional nutrients. This is a complex process and extremely difficult to include in models, so I am not suggesting the authors include it here, but maybe a different term for this flux can be used here. I am also not clear what ‘spread’ refers to here.

L 367 Table 2 Check the notations here - doe eta_CP refer to litter or soil CP?

L 460 by 1000 times, do you mean 1000 years?

L 564 Figure 3 Would it be possible to label the panels more clearly? I have to assume that the one on the left is fluxes and the one on the right, pools?

L 581 Figure 4  Could you discuss why a change in biomass production does not result in a change in plant C?

The authors present an attempt at describing the implementation of a phosphorus cycle into a major land surface model and its evaluation. They compare simulated carbon variables with few observations and provide another model-based quantification of P effects on NPP, C stocks under 2 years of elevated CO2 which is compared to existing model predictions. The inclusion of P cycles in ESM  is certainly a timely and important endeavor given the growing evidence of the importance of phosphorus cycling on land surface conditions and greenhouse gas fluxes. However, my main concern is that this study adds little to existing studies: (1)  the evaluation of the model falls short thereby no new insights could be gained; (2) the eCO2 experiment is a repetition of intermodel comparison of Fleischer et al. (2019) and provides no new insights. 

Potential new contributions could have been  (1) resolving root exudates. But (guessing from the incomplete model description) it seems a simply a relabelling of the 'excess NPP' (i .e NPP which cannot be allocated to new biomass growth (Thornton et al 2007, Goll et al 2012) as now 'root exudates'. On top of that, no attempt has been made at its evaluation and the authors seem to confuse observed BP as NPP (BP + root exudates). Thus, I am not sure anything new has been learnt here about. (2) the use of response of biota to nutrient addition from the AFEX experiment. 

Overall, several shortcomings have been identified which are listed in the following.

1 The model description is incomplete 

• it is not clear how the different configurations of JULES differ (C,CN,CNP). In reality the cycles are closely intertwined, so you must have made some simplification to switch them on and off. As interactions between NP cycles are not explained one can only guess how NP affects C fluxes. This prevents the reader from understanding the implications of the model results as major underlying model assumptions are not given.
• it is not clear which modelling approaches are novel and which are based on concepts from previous studies/models. The majority of process representation seem to be taken from earlier work (like early P work in JSBACH, CLM, ORCHIDEE, CABLE).  The authors fail repeatedly to credit earlier works (most references given are related to JULES) and to justify their modelling choice. 
• The presentation of model equations is poor making and many inaccuracies make it hard to follow (see minor points related to eq listed below).
• Model input parameteres are not given

2 Some of the assumptions / choices are in contrast to current understanding and consensus while no explanation was given.

• The assumption that CNP cycles are in steady-state with present day conditions (1999-2019) is not appropriate. There are multiple lines of arguments, that the historic increase in CO2 has led to a progressive limitation of nutrients (e,g, Luo et al 2004, Goll et al 2012, Penuelas et al 2013) and that present day land carbon cycle is not in equilibrium.  The non-steady state of the present day CNP cycles is accounted for in the majority of  modelling exercises (including Fleischer et al 2019. Trendy modelling protocol). The historic increase in CO2 is likely the more dominant factor affecting the present day state of C cycle compared to (progressive) NP limitation. Thus this omission is a major shortcoming, in particular as model predictions which account for this exist (i.e. Fleischer et al 2019). You should at least test what the implications of omitting this on the results are / better redo the whole experiment.
• There are several highly uncertain parameters in your model. It is not clear why you varied in the sensitivity test, only (a few) parameters which happen to be among the few observed ones and not choose more uncertain parameters? Besides, the impact of varying stoichiometry has been investigated in earlier models with a comparable plant P cycle (Goll et al 2012). 
• Some processes usually included in model have been omitted without giving any rationale why. E.g. Why is phosphatase mediated mineralization being omitted? It seems all major P models account for it. I don’t imply the  author must account for it, but if they choose not to, an explanation should be given. However, Fig2 indicates modelled soil organic P is quite high compared to observation, which would point towards missing biochemical mineralisation (which reduces organic P by enhancing its turnover) is problematic; Why is plant internal nutrient and carbon storage being omitted? Previous modelling studies showed the importance of accounting for such storage pools, see e.g.  Yang, Xiaojuan, et al. "Global evaluation of ELM v1 and the role of the phosphorus cycle and non-structural carbon in the historical terrestrial carbon balance." AGU Fall Meeting Abstracts. Vol. 2020. 2020. 
• What about N losses like leaching or erosion, inputs from atmospheric deposition? Are they omitted? why?

3 The evaluation is insufficient and conducted poorly

• The observed NPP is based on biomass increments (Fleischer et al 2019; SI table 1). It should be thus referred to BP and not NPP. Also the observed CUE is BP/GPP and NPP/GPP. There is little support that CNP performs better than the C version.
• The short period of 2 years considered in the evaluation is inappropriate. Fluxes like NPP and GPP show large interannual variations. The fit of modelled fluxes with a long-term mean of fluxes (e.g GPP from Fleischer et l 2019) could be merely by chance. The justification for using only two years ( line 470 ) is not plausible as you don’t use only the soil P measurement for evaluations (see first point). You should evaluate modelled variables over a longer time period.
• The soil P measurements were used to calibrate the model not to evaluate (misleading labelling in table 3)
• There are several datasets available to evaluate nutrient cycles in ESM. Some of the remote sensing products have a spatial resolution sufficiently high to compare to site simulations.  e.g. Sun et al 2020, Hou et al 2020. Sun, Y., Goll, D. S., Chang, J., Ciais, P., Guenet, B., Helfenstein, J., Huang, Y., Lauerwald, R., Maignan, F., Naipal, V., Wang, Y., Yang, H., and Zhang, H.: Global evaluation of the nutrient-enabled version of the land surface model ORCHIDEE-CNP v1.2 (r5986), Geosci. Model Dev., 14, 1987–2010, https://doi.org/10.5194/gmd-14-1987-2021, 2021.
• There seems to be eddy covariance tower nearby. Why hasn’t been any data other than long-term everage GPP being used to evaluate the model?
• It is misleading to refer to a control plot of AFEX as a fertilizer experiment. I assume there is no data available from the fertilizer experiment, otherwise you should take advantage of this data to evaluate the model (see main point). I would suggest dropping AFEX and referring to the plots as nearby plots of K34 tower.

more specific points:

Model description

• Which processes / C Fluxes are affected by N and P limitation?  How are interactions among N P limitations accounted for? Does P affect N fluxes and vice versa?
• Eq6,7,9,10: the flux decPi,n - you should specify the i here - otherwise it seems the same fluxes are substrate from two different pools / i.e double accounting of fluxes.
• Eq9,10: What are the factors 0.46 and 0.54 - how are they derived? They should be parameters listed in Table2.
• Eq11: It is not clear how the delta P_O_S equals the sum of all P subpools. Do you mean the change of all P subpools?
• 205-207: can you indicate the equation describing the fluxes  (as done before e.g. line 195)
• The calculation is missing for 
• total exudates,  and the subcomponents related to growth and spread (eq16)
• Plant P:C ratio (eq20)
• What is dC/dt in eq 19,20
• Eq19 how does this NPP differ from the one in  eq17 - or is it the same?
• Eq28: how does R_in related to R_n?
• Eq29: what is DEM_DPMn and DEM_RPMn? What is the rationale behind this formulation?
• Eq30: why is the plant demand a function of soil P availability and not plant P demand?
• Eq31-36: give the rationale behind the choice of equation.
• Table 1: misses variables e.g. desorption fluxes
• Table 2: misses value for parameters which are PFT or depth dependent
• Eq39,40 is nleaf the same as N?
• Line 248: : Zaehle & Friend has no P cycle, the references seems inappropriate
• Line 438: why do you need to cap these fluxes? Why can these max be assigned to observed stocks? You should explain.
• Line 458: why do you prescribe LAI? Isn’t that computed prognostically by JULES?

The choice of various different ways to label fluxes make it hard to read the equations. I would suggest adopting a more homogenous way. Here are some examples how inhomogeneous labelling:

• compare fluxes in eq9 
• Eq10 and eq12 use two different spellings to refer to desorption. Same goes for occlusion in eq 12 and 14. 
• For carbon fluxes (eq16) you again use also a mix of labels (e.g. greek letters for exudate, NPP. etc). 
• See labelling of CP ratio in eq1-3 vs eq19ff
• You use P subscript to refer to potential and phosphorus.
• Eq 42-44 vs equations before: concerning leaf, root, wood c mass
• I didn’t list all the inconsistencies here, there are more  but it shouldn't be the job of the referee ...

others

Line 23: N -> nitrogen

Line 24: not only tropical systems (Huo et al 2020 and references therein): Hou, E., Luo, Y., Kuang, Y. et al. Global meta-analysis shows pervasive phosphorus limitation of aboveground plant production in natural terrestrial ecosystems. Nat Commun 11, 637 (2020). https://doi.org/10.1038/s41467-020-14492-w

Line 30: you don’t use any data from the fertilizer experiment, only from the control plot. This is misleading.

Line 31-35: Given this study introduces a new model, I would expect more information on the model evaluation and performance, less on model prediction. 

Line 54: nutrient -> nitrogen

Line 58: you need to add references which indicates that much progress regarding has been made

Line 64/65: Reference is missing. E.g. Huo et al 2020 and references therein

Line 68-73: This formulation is unclear and the statement is not very nuanced.

Line 73:  Most model studies don’t show that P is stronger than N with respect to limiting the ecosystem response to eCO2 (e.g Goll et al 2012, Yang et al 2014)  Goll, D. S., Brovkin, V., Parida, B. R., Reick, C. H., Kattge, J., Reich, P. B., van Bodegom, P. M., and Niinemets, Ü.: Nutrient limitation reduces land carbon uptake in simulations with a model of combined carbon, nitrogen and phosphorus cycling, Biogeosciences, 9, 3547–3569, https://doi.org/10.5194/bg-9-3547-2012, 2012.

Line 73: what is separate knowledge?

Line 86-88: reference does not support this claim

Line 88: 2016 not so recent 

Line 133: this is not clear. I would assume light competition affects plant cover via photosynthesis and growth, and not directly.  You should be more clear how canopy processes are resolved.

Line140: it is not clear how N limitation on SOC can applied to vegetation. What thats’ the evidence/rationale for that ?

Figure 1: This figure caption is too short and does not reflect the text. E.g where is the depth information here? It would ne nice to indicate in the figure which equation explains which flux.

Line 159: you list only 7 pools, not 8.

Line 165: are they rigid ratios?

L226: this section also give C fluxes, should be relabelled

Line 228: I guess this should read: ‘In JULES-CNP, NPP represent the potential amount’ 

Line 228-235: this is hard to follow. It seems you list  definitions of potential NPP, actual NPP and BP as a justification for your modelling choice. I struggle to make the connection.

Line 244-246: maybe first follow the BP idea before. Also I think it needs some more explanation to model litter production as a function of NPP and not plant pools. 

Line 284: Zaehle & Friend has no P cycle. 

Line 300-306,  I can’t follow the equation because of the confusing way of labelling to C:P ratios (see point above).

Line 308: remove ‘However’

Line 381-382: this repeats info given before (line 376 ..)

Line 391-394: Please give information how data from plots were aggregated. You also indicate vegetation information from AFEX in Tab3.. How were vegetation  C stock derived.

Line 398: you need to explain how photosynthesis is implemented in JULES. Is it a big leaf approach? 

Line 428-436: you should explain which observable pools were assigned to which modelled pools.

Line: 459 : what’s the equilibrium criteria, e.g. do you apply the trendy protocol criteria? From the equations occluded and parental material pools cannot equilibrate. I guess you considered them outside the boundary of your model.

Line 460: 1000 times what? Years , timesteps, …

Line 470: why did you limit it to two years? Your evaluation is based on C variables, so your argumentation based on soil P measurement seems invalid. 

Line 483: Be clear about the timestep use here.

Line 500-503: how was this done exactly? E.g. Did you vary the ratio one by one? The relative differences between C:P ratios is critical for the availability of P. Did you test for this effect (as has been done in earlier studies (e.g. Goll et al 2012)).  The N:P ratios are usually much more constrained than the C:P ratios ? Did you take this information into account when varying C:P ratios?

Line 514: which pools exactly? How did you calibrate them?

Table 3: P pools were used to calibrate the model ( stated on line 514) not to evaluate them.

Table 5: The ‘observed NPP’ from Fleischer et al 2019 is based on biomass increments, thus you should compare it to modelled BP. This also affects CUE. 

Figure 3: you should add some measure of uncertainties of the obs, e.g. variation among plots

Line 555: how do you compute P limitation? I couldn’t find an equation.

Line 655: the units of GPP are misleading, i.e. monthly GPP as yr-1.

Line 674: why is it the higher soil water content which enhances uptake and not something? This needs to be demonstrated or explained in more detail.

L690: this is grassland study. You need to explain why this can be compared to a tropical forest.

L709: you should make the link between high competition for P and unclear role of P for plant CO2 response.

L711: it’s odd to refer to the site as well document. From your work it seems there is hardly any data available for model evaluation available.

Line 712: you included all major processes in detail ? Then the role of P should be clear now or?

Line 720-725: unclear formulations. Weren’t the soil P pools optimized in JULES?

Lin 7 27: inappropriate reference

Line 829: ‘a cornerstone’  - a  more humble forumation could be use here