This is a nicely written and well-executed work implementing phosphorus cycle and data assimilation framework into a process-based ecosystem model (TECO). The novelty of this work lie in the model development and its coupling with data assimilation. By comparing CNP model against observations and their respective C-only and CN coupled models, the authors show a superior performance of the newly developed model.

Overall, I enjoy reading the manuscript, and support its publication. There are several occasions where I think some justifications/modifications would further improve the quality of the manuscript. Below I list my main questions/concerns.

There are certain processes where N and P interact. For example, some models consider N to have an effect on soil P biochemical mineralization rate (e.g. ORCHIDEE, ELM, etc.). It does not seem that the authors have adopted these NP interacting processes in their model. Furthermore, others consider N and P to have joint effect on C processes, such as photosynthesis. In this work, it seems that nutrient effect on photosynthesis is realized via downregulation of leaf surface area. There have been some empirical relationships derived on how N and P affects photosynthetic traits (e.g. Vcmax, Jmax; Ellsworth et al., 2022; Walker et al,, 2014), and these relationships have been incorporated into models. What is the authors’ consideration on not following these conventional approaches?  

Solution P is part of labile P, and some work suggests the need to explicitly model solution P in addition to labile P (e.g. Reed et al., 2015). In this work, how does the author consider this suggestion and what is the rationale for only simulating labile P?

It seems that the data assimilation framework was only applied to the CNP model, and then the authors reported CN and C models to overestimate observations. I find this logic to be a bit problematic. Without data assimilation, does CNP model still achieve better match with observations? Alternatively, how does C-only model coupled with data assimilation perform relative to observations? If it can achieve similar performance as compared to CNP model, what benefits of having a CNP model?

This study developed a coupled carbon-nitrogen-phosphorus model, TECO-CNP Sv1.0, based on the Terrestrial ECOsystem (TECO) model. The developed model was used to simulate C, N, and P pools and fluxes in a phosphorus-limited subtropical forest site in East China. In addition, a parameter optimization algorithm was also incorporated into the model framework to improve the model's performance. Overall, the manuscript provides detailed information on the model structure, parameters, and performance. However, I still have some questions on the soil pool structure and calibration processes of the model.

Four inorganic P pools, including labile P, sorbed P, secondary P, and occluded P, are set in TECO-CNP. This structure is different from other CNP models. For example, the labile P pool in ORCHIDEE-CNP includes both dissolved and sorbed P. In CLM-CNP, inorganic P pools include labile P (including solution P), secondary P, and occluded P. In a global P dataset developed by He et al. (2023) Biogeosciences, the soil inorganic P is divided into labile P, moderate P, and occluded P. Can you explain the differences in inorganic P pool structure among these models? I am confused about the definition of labile P pool. In addition, how did you initialize these inorganic P pools？

What is the advantage of TECO-CNP compared with other CNP models?

How were the soil P pools at the Tiantong site measured? Did you compare the measured soil P pool with other studies? They seem lower than other studies (Fig. 5c).

The simulation results of C-only, CN, and CNP versions were compared in this study to prove the good performance of the CNP model. Did you calibrate these three versions individually? Was the same parameter optimization algorithm applied to all three versions? Or you just calibrate the only CNP model, and apply the same parameters to other versions. Did you simulate the C cycle in this site by using C-only or CN versions before the development of CNP? How well did these two models perform? Many parameters are constant values, such as Vre. I guess these parameters were not calibrated but derived from the literature. Were these parameters suitable for subtropical forest ecosystems?

The vegetation and soil pools simulated by the three versions were compared in section 3.1. What about C and nutrient fluxes? In Fig.7, I cannot identify the simulated NEE by the three model configurations.

Minor

L491-493. Please list the equations of P loss from SOM pools.

Equation 54. What is the meaning of Pl

L607. Please correct the reference of Xu et al.

Fig 8. What is the meaning of the posterior distribution of parameters? Do they change with time?

Table. How did you identify these target parameters? Did you conduct a sensitivity analysis?

Reference: He, Xianjin, et al. "Global patterns and drivers of phosphorus fractions in natural soils." Biogeosciences 20.19 (2023): 4147-4163.