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?