TECO-CNP Sv1.0: a coupled carbon-nitrogen-phosphorus model with data assimilation for subtropical forests

Subtropical forests play a crucial role in the global carbon cycle, yet their carbon sink capacity is significantly constrained by phosphorus availability. Models that omit phosphorus dynamics risk overestimating carbon sinks, potentially undermining the scientific basis for carbon neutrality strategies. In this study, we developed TECO-CNP Sv1.0, a coupled carbon-nitrogen-phosphorus model based on the Terrestrial ECOsystem (TECO) model, which explicitly captures key biogeochemical interactions and nutrient-regulated carbon cycling. The model simulates how plant growth and carbon partitioning respond to both external soil nutrient availability and internal physiological constraints, enabling plant acclimation to varying nutrient conditions. Using observations from a phosphorus-limited subtropical forest in East China, we first evaluated the model's performance in estimating state variables with empirically calibrated parameters. Compared to the C-only and coupled C-N configurations, the CNP model more accurately reproduced the observed pools of plant and soil C, N, and P. To systematically optimize model parameters and reduce uncertainties in predictions, we further incorporated a built-in data assimilation framework for parameter optimization. The CNP model with optimized parameters significantly improved carbon flux estimates, reducing root mean square errors and enhancing concordance correlation coefficients for gross primary productivity, ecosystem respiration, and net ecosystem exchange. By explicitly incorporating phosphorus dynamics and data assimilation, this study provides a more accurate and robust framework for predicting carbon sequestration in phosphorus-limited subtropical forests.

1 Introduction

Accurately representing phosphorus (P) cycling in land surface models (LSMs) is crucial for projecting terrestrial carbon sink dynamics under climate change (Wieder et al., 2015). As an essential element, P availability regulates plant growth and ecosystem productivity (Walton et al., 2023; Vitousek et al., 2010). For instance, nutrient addition experiments in an old-growth Amazon rainforest demonstrated that net primary productivity increased exclusively with P addition (Cunha et al., 2022). Likewise, in subtropical mature forests, soil P availability was found to exert dominant control over plant functional traits at both species and community levels (Cui et al., 2022). Recent global syntheses have revealed a more widespread distribution of terrestrial P limitation than previously recognized (Hou et al., 2020; Du et al., 2020; Xia and Wan, 2008; Elser et al., 2007). More concerning is that P limitation is expected to intensify (Wang et al., 2023a; Luo et al., 2022) due to factors such as N deposition-induced N : P stoichiometric imbalance (Peng et al., 2019; Lu and Tian, 2017; Du et al., 2016; Peñuelas et al., 2013) and reduced P availability under elevated CO₂ concentration (Wang et al., 2023b). Consequently, incorporating P limitation into LSMs has become a pressing challenge for improving carbon cycle projections (Fisher and Koven, 2020; Achat et al., 2016; Reed et al., 2015).

To address this challenge, several modeling groups have incorporated a prognostic P cycle into their existing frameworks over the past decade, including CASACNP (Carnegie-Ames-Stanford Approach; Wang et al., 2010), JSBACH (Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg; Goll et al., 2012), CLM-CNP (Community Land Model; Yang et al., 2014), among others. These pioneering efforts in coupled carbon-nitrogen-phosphorus (C-N-P) modeling have laid a solid foundation for increasing incorporation of P cycling in LSMs (e.g., Goll et al., 2017; Nakhavali et al., 2022) and demographic vegetation models (Knox et al., 2024), shedding light on how P limitation constrains ecosystem productivity under elevated atmospheric CO₂ (Wang et al., 2024b; Fleischer et al., 2019; Medlyn et al., 2016). However, current C-N-P models often yield “right answers for wrong reasons” (Jiang et al., 2024a), largely due to two key limitations: (1) calibration and validation data are predominantly derived from a narrow range of ecosystems, with most coupled C-N-P models relying on in-situ data from tropical regions, particularly Hawaii and the Amazon (e.g., Nakhavali et al., 2022; Yang et al., 2014; Goll et al., 2012, 2017; Zhu et al., 2016b), and (2) oversimplified representations of P cycling processes (Achat et al., 2016; Reed et al., 2015), such as the absence of physiological mechanisms governing vegetation P uptake (Jiang et al., 2019). Addressing these gaps requires advancing the coupled C-N-P model with improved mechanistic process-based representations and broader ecosystem applicability (Jiang et al., 2024a).

Subtropical forest ecosystems are recognized as important carbon sinks in the global carbon cycle (Pan et al., 2024; Keenan and Williams, 2018; Yu et al., 2014). In particular, East Asian monsoon subtropical forests exhibit high carbon sink capacity, with an average net ecosystem productivity of about 400 g C m⁻² yr⁻¹ (Yu et al., 2014). These ecosystems are likely subject to substantial phosphorus limitation, as evidenced by a meta-analysis of nutrient addition experiments showing that forest productivity exhibits the strongest standardized response to P addition in the subtropical regions (25–40 latitude; Hou et al., 2021). Moreover, intensive nitrogen deposition may further exacerbate P limitation (Zhu et al., 2016a; Yu et al., 2014). Accurately projecting the future carbon sink capacity of subtropical forests is crucial for assessing their role in climate change mitigation (Friedlingstein et al., 2023; Requena Suarez et al., 2019; Grassi et al., 2017). However, substantial uncertainties remain in current model projections of subtropical carbon dynamics (Wei and Xia, 2024), highlighting the urgent need for improved carbon cycle predictions through better representation of coupled C-N-P interactions in these regions.

In this study, we develop TECO-CNP Sv1.0, an advanced version of the Terrestrial ECOsystem (TECO) model (Weng and Luo, 2008, 2011), incorporating detailed mechanistic representations of coupled C-N-P cycling processes, such as dynamic plant growth response to soil available nutrient through modified growth rates and allocation patterns, and the combined physical and physiological controls on phosphorus uptake. Additionally, we integrated a data assimilation module based on a Bayesian probabilistic inversion approach (Xu et al., 2006; Ma et al., 2017; Shi et al., 2016, 2018; Zhou et al., 2021), providing an efficient framework for model reparameterization and broader applications. Based on comprehensive observations from a P-limited subtropical evergreen broadleaf forest in eastern China, we further test two key hypotheses: (1) the CNP model can reproduce ecosystem state variables through traditional spin-up and manual parameter tuning, and (2) the built-in data assimilation system can substantially improve carbon flux predictions.

2 Materials and Methods

2.1 TECO developments

The TECO-CNP model has evolved from its precursor, the Terrestrial Ecosystem model (TECO, Weng and Luo, 2008). The TECO model is a process-based ecosystem model encompassing eight organic carbon pools and a plant non-structural carbohydrate (NSC) pool (Weng and Luo, 2008). The representation of the NSC pool in TECO is advantageous for capturing the seasonal decoupling of growth and nutrient acquisition within plants (Zavišić and Polle, 2018; Jones et al., 2020) and for managing carbon that is not utilized for plant growth under nutrient-limited conditions (Nakhavali et al., 2022; Haverd et al., 2018). The TECO model has been part of model intercomparison ensembles (Zaehle et al., 2014; De Kauwe et al., 2014) and has been applied across diverse ecosystem types, such as grassland (Weng and Luo, 2008; Zhou et al., 2021), temperate coniferous forests (Luo et al., 2003; Weng and Luo, 2011; Jiang et al., 2017) and deciduous broadleaf forests (Jiang et al., 2017) and northern peatland (Ma et al., 2017, 2022; Huang et al., 2017).

Simplified N and P cycling were incorporated in the TECO successively (Shi et al., 2016; Du et al., 2018, 2021), where the structure of the carbon processes was expressed as a matrix form (Luo et al., 2003; Xu et al., 2006; Weng and Luo, 2011). Thus, the photosynthesis was simulated aided by an external model; for instance, Shi et al. (2016) utilized MAESTRA to generate the gross primary productivity. The processes related to the N and P cycle were only represented in a parsimonious way in the matrix versions. For example, the nutrient uptake process was simplified at a constant rate, and the interactions of carbon, nitrogen, and phosphorus were treated implicitly (Shi et al., 2016; Du et al., 2021).

In this study, we developed TECO-CNP, a coupled C-N-P model based on the full version of TECO, which fundamentally differs from previous matrix-based approaches. This new model explicitly represents the mechanistic processes of nutrient cycling (Sect. 2.2), with a focus on the regulation of carbon cycling by nutrients. Specifically, the model incorporates four key nutrient-carbon interactions: (1) growth rate limitations controlled by internal plant nutrient concentrations and nutrient supply-demand relationships; (2) allocation patterns dependent on nitrogen and phosphorus availability; (3) decomposition processes constrained by microbially-mediated nutrient availability; and (4) carbon costs associated with nutrient uptake and fixation. These process-based implementations, which aim to provide a more realistic representation of terrestrial biogeochemical cycles, are described in detail in the following sections.

2.2 Model description

We introduce a comprehensive biogeochemical N and P cycle into the full TECO, named TECO-CNP Sv1.0. Key processes of N and P cycling and their interactions with the carbon cycle have been represented using reliable mechanistic assumptions based on our experimental measurements or validated by state-of-the-art LSMs. In the following sections, we first document an overview of the carbon cycle and highlight the effects of nutrient limitation on the carbon cycle in Sect. 2.2.1. We then describe the shared and specific N and P cycling processes in Sect. 2.2.2 and 2.2.3, respectively.

2.2.1 Nutrient-limited carbon cycle

The carbon cycle in the new model builds upon the TECO model, incorporating processes such as photosynthesis, plant growth controlled by allocation and phenology, autotrophic and heterotrophic respiration, litter production, and carbon transfer (Fig. 1). See Luo et al. (2003) and Weng and Luo (2008) for detailed descriptions. These processes regulate the dynamics of plant, litter, and soil pools (Fig. 2). Nutrients directly or indirectly constrain them. For instance, plant growth rates and carbon allocation strategies are directly influenced by internal nutrient availability within pools and the availability of soil-accessible nitrogen and phosphorus. Additionally, resource limitations adhere to Liebig's law of the minimum, where the nutrient-constrained process is hindered only by the most limiting resources (Rastetter, 2011).

Figure 1 The schematic diagram of the biogeochemical processes of the carbon, nitrogen, and phosphorus cycles and associated interactions in TECO-CNP. Representation of carbon cycling processes controlled by nitrogen and phosphorus in TECO-CNP. Solid lines indicate carbon cycling processes (labelled 1–7) comprise (1) photosynthesis, (2) carbon allocation, (3) plant growth, (4) autotrophic respiration, (5) litter production, (6) carbon transfer, and (7) heterotrophic respiration. These processes are controlled directly by nitrogen and phosphorus (black control characters) or indirectly (colorless control characters). Dashed lines indicate the common processes that control the dynamics of soil-available nitrogen and phosphorus, simplified as plant uptake, mineralization, immobilization, biogeochemical mineralization, and external input and loss. Irregular pink shapes represent competition for soil available nitrogen and phosphorus between plants and microorganisms. Min., mineralization; BMin., biochemical mineralization; Imm., Immobilization.

Figure 2 Model structure of TECO-CNP. The model represents the nine organic carbon, nitrogen, and phosphorus stocks within the plant (denoted as Q1–4), litter (Q5–7), and soil (Q7–9). Fluxes among these organic pools are depicted by black arrows. Specific N and P fluxes are indicated by dark red arrows, with associated processes labeled accordingly. Min denotes mineralization, and Imm denotes immobilization. The circled numbers (1–7) correspond to the carbon cycling processes in Fig. 1.

The canopy-level photosynthesis is simulated using a two-leaf model, which consists of a radiation sub-model and a coupled sub-model of stomata, photosynthesis, and transpiration for both sunlight and shaded leaves (Wang and Leuning, 1998). Leaf photosynthesis is estimated by the equations derived from the Farquhar model (Farquhar et al., 1980) and a stomatal conductance model (Ball et al., 1987; Leuning, 1995). The photosynthesis of a single leaf is then scaled up to the canopy level (Wang and Leuning, 1998). We hypothesize that plant photosynthesis is downregulated as photosynthetic surface area decreases when nutrient limits plant growth. Plant growth is adjusted based on the nutrient limitation factor calculated at each time step, meaning that plants tend to reduce growth under low nutrient conditions to avoid nutrient deficiency within the organism (Veneklaas et al., 2012). Accordingly, the nutrient-constrained growth rate (GP_a) is dependent on the potential growth rate (GP_p) and nutrient limitation scalar for plant growth (L_GP) as the following equation:

$$\begin{array}{}\text{(1)}& {\mathrm{GP}}_{\mathrm{a},i}={\mathrm{GP}}_{\mathrm{p},i}\cdot L_{\mathrm{GP}},\end{array}$$

where subscript i indicates leaf (i=1), wood (i=2), root (i=3) or reproduction (i=4) (Table 1). The difference between actual and potential plant growth is referred to as excess carbon, which implicitly represents the carbon lost from the NSC pool through various pathways to cope with nutrient limitations.

Table 1 Variables for carbon cycling processes in TECO-CNP.

^* i indicates leaf (i=1), wood (i=2), root (i=3) or reproduction (i=4), j indicates metabolic litter (j=5) or structure litter (j=6), and m indicates fast SOM (m=7), slow SOM (m=8) and passive SOM (m=9).

The nutrient limitation scalar for plant growth incorporates both the nutrient status of plant tissues and soil nutrient supply (Fig. 1b), which can be expressed as:

$$\begin{array}{}\text{(2)}& L_{\mathrm{GP}}=L_{\mathrm{in},\mathrm{leaf}}\cdot L_{\mathrm{sp}},\end{array}$$

where L_in,leaf and L_sp represent the nutrient limitation factors derived from leaf nutrient concentration (Eqs. 3–5) and the nutrient demand-supply process (Eqs. 6–8), respectively. Shifts in leaf nutrient concentrations act as a potential limiting factor for plant growth, implying the mechanism by which changes in leaf nutrient concentration can impact photosynthesis (Ellsworth et al., 2022; Sterner and Elser, 2002). Description of limitation factors that account for plant tissue's nutrient concentration can be given by:

$$\begin{array}{}\text{(3)}& {\displaystyle }{L}_{\mathrm{in},i}=min\left(L_{\mathrm{in},\mathrm{N},i},L_{\mathrm{in},\mathrm{P},i}\right),\text{(4)}& {\displaystyle }{L}_{\mathrm{in},\mathrm{N},i}={\displaystyle \frac{R_{\mathrm{N},i}}{R_{\mathrm{N},i}+k_{\mathrm{CN}}}},\text{(5)}& {\displaystyle }{L}_{\mathrm{in},\mathrm{P},i}={\displaystyle \frac{R_{\mathrm{P},i}}{R_{\mathrm{P},i}+k_{\mathrm{CP}}}},\end{array}$$

where R_N and R_P represent the C : N ratios and C : P ratios, respectively. k_CN and k_CP are empirical parameters. A study by Cui et al. (2022) reveals that the Tiantong site is identified as a P-limited ecosystem, as indicated by the leaf N : P thresholds from Koerselman and Mueleman (1996). Thus, we adopted the values of k_CP (0.0006 g P g C⁻¹) in Wang et al. (2010) to achieve N limitation when N : P < 16 (g N g P⁻¹), and otherwise, plant growth is limited by P. k_CN (0.01 g N g C⁻¹) is given based on the results of Linder and Rook (1984).

The nutrient demand-supply limitation factor is calculated as a function of plant nutrient uptake and demand. When nutrient demand is not satisfied, the value of the limitation factor falls below one, thereby impacting plant growth. This assumption aligns with field findings that reveal an increase in plant productivity following nutrient addition (Cunha et al., 2022; Liang et al., 2021). Description of nutrient demand-supply limitation factor (L_sp) can be given by:

$$\begin{array}{}\text{(6)}& {\displaystyle }{L}_{\mathrm{sp}}=min\left(L_{\mathrm{sp},\mathrm{N}},L_{\mathrm{sp},\mathrm{P}}\right),\text{(7)}& {\displaystyle }{L}_{\mathrm{sp},\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp (-\mathrm{12}\cdot \frac{F_{\mathrm{up},\mathrm{N}}}{F_{\mathrm{dm},\mathrm{N}}}+\mathrm{6})}},\text{(8)}& {\displaystyle }{L}_{\mathrm{sp},\mathrm{P}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\exp (-\mathrm{12}\cdot \frac{F_{\mathrm{up},\mathrm{P}}}{F_{\mathrm{dm},\mathrm{P}}}+\mathrm{6})}},\end{array}$$

where F_up,N and F_up,P represent plant nutrient uptake for N and P, respectively, which is determined by both supply and demand (Eq. 23). F_dm,N and F_dm,P represent the plant required N and P to sustain a given NPP (Eq. 24). We implemented a logistic function to represent the phosphorus limitation factor, which provides a more mechanistically sound representation of nutrient limitation compared to the simple linear ratio. This formulation ensures a smooth transition between phosphorus-limited and phosphorus-sufficient conditions, with values bounded between zero and one. The coefficients were carefully selected to maintain appropriate sensitivity in the transition zone while avoiding unrealistic sharp thresholds. This sigmoidal response more accurately reflects the gradual physiological adjustments of plants to varying nutrient availability. It is consistent with a theoretical understanding of the effects of nutrient limitation on plant growth. The method of determining whether plants are nutrient-limited based on the supply-demand method is widely employed in many models, for example, CASACNP (Wang et al., 2010), CLM-CNP (Yang et al., 2014), and ORCHIDEE (revision 4520; Goll et al., 2017).

The carbohydrates available for plant growth will be redistributed among the plant pools based on their actual growth rates. A prescribed proportion of those allocated to reproductive processes (Sitch et al., 2003; Smith et al., 2001), such as flower formation, fruit development, and seed production, is stored in the reproductive pool. Vegetation growth is assumed to take priority over reproduction (Zust et al., 2015; Tang et al., 2021). Thus, the plant's reproductive allocation is zero when the leaf area index (LAI) is below the minimum threshold. When LAI exceeds the minimum threshold, 12 % of the available carbon is allocated to the reproduction pool. The remaining carbon is subsequently distributed among leaf, wood, and root pools based on a resource limitation allocation scheme.

The dynamic allocation for leaf, wood, and root is regulated by light availability, soil water supply, canopy phenological status (Luo et al., 1995; Denison and Loomis, 1989; Arora and Boer, 2005), and plant's internal nutrient status (Fig. 1b). This allocation strategy permits a reduction in photosynthetic surface area and enhanced root growth under nutrient limitation, exemplifying a structural adjustment in line with the observations (Keith et al., 1997; Thomas et al., 2015; Yan et al., 2016). The allocation fractions for leaf, wood, and root are given by:

$$\begin{array}{}\text{(9)}& {\displaystyle }{b}_{\mathrm{C},\mathrm{leaf}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{L}}\cdot L_{\mathrm{in},\mathrm{leaf}}}{\mathrm{1}+\mathit{\omega }\left(\mathrm{2}-\mathit{\kappa }-W\right)}},\text{(10)}& {\displaystyle }{b}_{\mathrm{C},\mathrm{wood}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{w}}\cdot L_{\mathrm{in},\mathrm{wood}}+\mathit{\omega }\left(\mathrm{1}-\mathit{\kappa }\right)}{\mathrm{1}+\mathit{\omega }\left(\mathrm{2}-\mathit{\kappa }-W\right)}},\text{(11)}& {\displaystyle }\begin{array}{rl}{b}_{\mathrm{C},\mathrm{root}}& ={\displaystyle \frac{(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{L}}\cdot L_{\mathrm{in},\mathrm{leaf}}-{\mathit{\epsilon }}_{\mathrm{w}}\cdot L_{\mathrm{in},\mathrm{wood}})+\mathit{\omega }\left(\mathrm{1}-\mathit{\kappa }\right)}{\mathrm{1}+\mathit{\omega }\left(\mathrm{2}-\mathit{\kappa }-W\right)}}\\ & =\mathrm{1}-b_{\mathrm{C},\mathrm{leaf}}-b_{\mathrm{C},\mathrm{wood}},\end{array}\end{array}$$

where b_C, leaf, b_C, wood and b_C, root represent the carbon fractions available for growth allocated to leaf, wood, and root, respectively. W is the root zone soil water availability stress factor (Arora and Boer, 2005). The soil water availability is weighted by the existing fraction of roots in each soil layer (Weng and Luo, 2008; Arora and Boer, 2005). κ represents the availability of light (Arora and Boer, 2005). Parameters ε_w, ε_L, and ω are calibrated based on the broadleaf evergreen PFT parameters given in Arora and Boer (2005). L_in, wood and L_in, leaf represent the limitation factor determined by the nutrient status of tissues (Eqs. 3–5), designed to capture the reduction of carbon allocated to leaf and wood as an adaptation to nutrient limitation (Binkley et al., 1995; Yan et al., 2016) and the negative correlation between fine root biomass and soil fertility (Fortier et al., 2019).

Canopy phenology is represented by annual variation in LAI. The beginning of a growing season is determined by growing degree days. Leaf senescence is triggered by low air temperatures and soil moisture (Arora and Boer, 2005), resulting in a reduction of LAI. The litter production rates of wood and roots are prescribed. The phenological parameters are adjusted according to the vegetation characteristics in the studied evergreen forest (Table S1 in the Supplement).

Carbon transfer between litter pools and soil organic pools through microbial decomposition (Luo and Reynold, 1999; Weng and Luo, 2008). The decomposition of litter and soil organic matter (SOM) is diminished when the amount of available inorganic N and P restricts nutrient immobilization during decomposition:

$$\begin{array}{}\text{(12)}& {\displaystyle }{D}_{\mathrm{a},j}=D_j\cdot L_{\mathrm{de}},\text{(13)}& {\displaystyle }{D}_{\mathrm{a},m}=D_m\cdot L_{\mathrm{de}},\end{array}$$

where j indicates metabolic litter (j=5) or structure litter (j=6), and m indicates fast SOM (m=7), slow SOM (m=8) and passive SOM (m=9). D_a is the nutrient-constrained decomposition rate, and D is the default decomposition rate controlled by the soil temperature and moisture (Weng and Luo, 2008). L_de is the limiting factor of decomposition, and the calculation involves dividing the un-limited net mineralization rate by the size of the inorganic nutrient pool, which can be addressed in the following equations:

$$\begin{array}{}\text{(14)}& {\displaystyle }{L}_{\mathrm{de},\mathrm{N}}=max\left(\mathrm{0},\mathrm{1}+{\displaystyle \frac{F_{\mathrm{N},\mathrm{net}}^{\prime }}{N_{\min }}}\right),\text{(15)}& {\displaystyle }{L}_{\mathrm{de},\mathrm{P}}=max\left(\mathrm{0},\mathrm{1}+{\displaystyle \frac{F_{\mathrm{P},\mathrm{net}}^{\prime }}{P_{\mathrm{lab}}}}\right),\text{(16)}& {\displaystyle }{L}_{\mathrm{de}}=min\left(L_{\mathrm{de},\mathrm{N}},L_{\mathrm{de},\mathrm{P}}\right),\end{array}$$

where $F_{\mathrm{N},\mathrm{net}}^{\prime }$ and $F_{\mathrm{P},\mathrm{net}}^{\prime }$ represent the net mineralization rate for nitrogen and phosphorus, respectively, assuming no nutrient limitation on mineralization (Wang et al., 2010).

2.2.2 Shared processes in the N and P cycle

The shared processes of N and P cycling include plant uptake, resorption, allocation, transfer from plant to the soil through litterfall, and transfer between organic litter and soil pools via biological mineralization and N, P biological immobilization (Fig. 2). Underlying these processes, TECO-CNP incorporates two key N-P interaction mechanisms: P uptake regulated by a nutrient balance scalar and a cost-benefit approach-based regulation of phosphatase production. To avoid duplication, the shared processes were described collectively.

The organic N (Q_N) and P pools (Q_P) are coupled with carbon pools through flexible stoichiometry within plant, litter, and soil pools. Inorganic nutrient components consist of one inorganic soil N pool (N_min) and four inorganic soil P pools, including labile P (P_lab), sorbed P (P_S), secondary P (P_SS), and occluded P (P_O). Labile P represents readily bioavailable inorganic phosphate for biotic uptake and soil leaching. Sorbed P is weakly bound to soil surfaces in dynamic equilibrium with labile P. Through petrochemical processes, sorbed P transforms into secondary mineral P, which eventually becomes occluded P with minimal bioavailability. The key variables of N and P cycling are listed in Tables 2 and 3, respectively. Key parameter values were derived from site-specific field observations of plant functional traits and biogeochemical properties, as well as from validated studies chosen based on careful consideration of the ecosystem characteristics of the study site (Tables S1–S3 in the Supplement and Table 4).

Table 2 Common variables for N and P cycle modeling.

* χ indicates N or P. Subscripts i, m, and j refer to the values in Table 1.

Table 3 Specific variables in N and P cycle modeling.

Table 4 Parameters for nitrogen and phosphorus cycling in TECO-CNP.

^a For reference codes, see Table S4. ^b Ref 11 – Mollier et al. (2008).

The initial size of the organic nutrient pool is determined by the size of the carbon pool and the carbon-to-nutrient ratio. The dynamics of organic nitrogen and phosphorus transfer from donor to recipient pools within plants, litter, and soil are coupled with carbon cycling through flexible stoichiometry. The dynamics of plant nutrient pools can be expressed as:

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{Q}_{\mathit{\chi },i}\left(t\right)=F_{\mathrm{new},\mathit{\chi },i}-Q_{\mathrm{C},i}\cdot {\mathit{\tau }}_i\cdot R_{\mathit{\chi },i}^{-\mathrm{1}},\text{(18)}& {\displaystyle }{F}_{\mathrm{new},\mathit{\chi },i}=F_{\mathrm{new},\mathrm{C},i}\cdot R_{\mathit{\chi },i}^{-\mathrm{1}}+(Q_{\mathrm{C},i}\cdot R_{\mathit{\chi },i,\mathrm{0}}^{-\mathrm{1}}-Q_{\mathrm{C},i}\cdot R_{\mathit{\chi },i}^{-\mathrm{1}})\end{array}$$

where subscript χ = N, P, $F_{\mathrm{new},\mathit{\chi },i}$ represents the newly input nutrients from non-structural nutrient pool to sustain plant growth (Table 2), $F_{\mathrm{new},\mathrm{C},i}$ is determined by the newly input carbon from the NSC pool to plant pool i and stoichiometric ratios (Eq. 25). $R_{\mathit{\chi },i,\mathrm{0}}$ and R_χ,i denote the initial and updated C : N (or C : P) ratios of plant pool i. Q_C,i and τ_i represent the carbon pool size and turnover rate of plant pool i. The dynamically constrained nutrient redistribution process in plants (Eq. 18) follows the principles of stoichiometric homeostasis theory (Sterner and Elser, 2002) and helps avoid excessive flexibility in stoichiometry during model simulations (Meyerholt and Zaehle, 2015; Goll et al., 2017).

Nutrients newly acquired from root uptake (F_up,χ) and tissue resorption (F_res,x) enter the labile nutrient pool, which buffers the nutrient dynamics and mitigates imbalances between supply and demand (Weng et al., 2017). Thus, the dynamics of plant labile nutrient pools are modeled as:

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathrm{NS}}_{\mathit{\chi }}\left(t\right)=F_{\mathrm{up},\mathit{\chi }}+F_{\mathrm{res},\mathit{\chi }}-\sum _i{F}_{\mathrm{new},\mathit{\chi },i}.\end{array}$$

Since the reproduction pool is designed as a long-term pool supporting a series of reproductive events, from flower bud formation to fruiting, no resorption is prescribed in this pool. The relocation of nutrients from senesced plant tissues (F_res,χ) is modeled as:

$$\begin{array}{}\text{(20)}& F_{\mathrm{res},\mathit{\chi }}=\sum _i{\mathit{\alpha }}_{\mathit{\chi }}\times Q_{\mathrm{C},i}\cdot {\mathit{\tau }}_i\cdot R_{\mathit{\chi },i}^{-\mathrm{1}}(i\ne \text{reproduction}),\end{array}$$

where α_χ is the resorption rate and the second term represents the loss of carbon from plant pool i (Table 4). We assume that the different plant organs have the same and fixed resorption rate to simplify this process. Additionally, we prescribe a higher resorption rate for P at 0.4 compared to N at 0.2, considering the higher phosphorus use efficiency in the P-limit habitat (Xu et al., 2020).

Litter nutrient dynamics is given by:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{Q}_{\mathit{\chi },j}\left(t\right)=F_{\mathrm{P}\mathrm{2}\mathrm{L},ij}-Q_{\mathrm{C},j}\cdot {\mathit{\tau }}_j\cdot R_{\mathit{\chi },j}^{-\mathrm{1}},\end{array}$$

where F_P2L,ij represent the nutrient flux from plant pool i to metabolic litter (j=5) and structure litter (j=6):

$$\begin{array}{}\text{(22)}& F_{\mathrm{P}\mathrm{2}\mathrm{L},ij}=\left\{\begin{array}{l}\left(\mathrm{1}-{\mathit{\alpha }}_{\mathit{\chi }}\right)Q_{\mathrm{C},i}\cdot R_{\mathit{\chi },i}^{-\mathrm{1}}\cdot {\mathit{\tau }}_i\cdot r_{i,j},\\ i=\mathrm{1},\mathrm{2},\mathrm{3}\\ \\ Q_{\mathrm{C},i}\cdot R_{\mathit{\chi },i}^{-\mathrm{1}}\cdot {\mathit{\tau }}_i\cdot r_{i,j},\\ i=\mathrm{4}\end{array},\right.\end{array}$$

where r_i,j represents the fraction of plant carbon to different litter pools.

The TECO-CNP model exclusively considers the active uptake of inorganic P through specialized transporters on the root surface (Schachtman et al., 1998), as inorganic P is the form most readily absorbed by plants (Bieleski, 1973). Plants possess specific transporters and mechanisms dedicated to transmembrane transport, ensuring they can acquire P even from soil solutions with low P concentrations, where the P concentration can be as low as one-thousandth of the intracellular concentration (Schachtman et al., 1998). Therefore, we assume that plants absorb only inorganic P from the soil. Similarly, we also only consider the plant uptake of inorganic N. Soil labile nutrients taken up by plants are generally contingent upon both nutrient demand for growth (Wang et al., 2010) and root uptake capacity (Grant et al., 1999, 2001; Goll et al., 2017) that are related to root morphology and soil nutrient concentrations. The nutrient demand-supply scheme has been widely employed in most coupled C-nutrient models (Achat et al., 2016). We assume plants will not consume nutrients beyond their luxury consumption demand for assimilating nutrients (Van Wijk and Williams, 2003; Chapin, 1980). Therefore, the F_up,x is determined by either the nutrient demand (F_dm,χ) or the nutrients supplied by soil (F_sp,χ), whichever is lower:

$$\begin{array}{}\text{(23)}& F_{\mathrm{up},\mathit{\chi }}=\left\{\begin{array}{ll}{F}_{\mathrm{dm},\mathit{\chi }}& \left(F_{\mathrm{dm},\mathit{\chi }}<{\mathrm{UC}}_{\mathrm{root},\mathit{\chi }}\right)\\ F_{\mathrm{sp},\mathit{\chi }}& \left(F_{\mathrm{dm},\mathit{\chi }}>{\mathrm{UC}}_{\mathrm{root},\mathit{\chi }}\right)\end{array}\right..\end{array}$$

The F_dm,χ is determined by the invested carbon for newly formatted tissues (NPP_i) and C : nutrient ratios. The actual demand is considered as the difference between the demand for growth and resorption capacity:

$$\begin{array}{}\text{(24)}& {\displaystyle }{F}_{\mathrm{dm},\mathit{\chi }}={\sum }_i{\displaystyle \frac{F_{\mathrm{new},\mathrm{C},i}}{R_{\mathit{\chi },i}}}-F_{\mathrm{res},\mathit{\chi }},\text{(25)}& {\displaystyle }{F}_{\mathrm{new},\mathrm{C},i}={\mathrm{NPP}}_{\mathrm{a}}\cdot b_{\mathrm{C},i}^{\prime },\end{array}$$

where NPP_a represents the net primary productivity derived from actual plant growth (Eq. 1), $b_{\mathrm{C},i}^{\prime }$ denotes the b_C,i (Eqs. 9–11) specifically influenced by the leaf phenology (Weng and Luo, 2008).

The nutrients supplied to plants from the soil depend not only on the amount of P in the soil but also on soil conditions and the root uptake capacity. We implemented the function of F_sp,χ as described by Goll et al. (2017), and it is calculated by the function of root biomass (BM_root), and root uptake capacity (uc_root,χ), soil temperature scalar (f_T) and the nutrient balance scalar (f_χ,ratio) as follows:

$$\begin{array}{}\text{(26)}& F_{\mathrm{sp},\mathit{\chi }}={\mathrm{BM}}_{\mathrm{root}}\cdot {\mathrm{uc}}_{\mathrm{root},\mathit{\chi }}\cdot f_{\mathrm{T}}\cdot f_{\mathit{\chi },\mathrm{ratio}}.\end{array}$$

The linear index scalar f_χ,ratio regulates the balance between C, N, and P by constraining nutrient uptake rates based on prescribed maximum ratios (Eqs. 27–28), thereby preventing resource overconsumption (Goll et al., 2017). Experiments have shown that N addition enhances the uptake of both N and P, suggesting a benefit for P uptake when more N is available (Zhu et al., 2021). Thus, we assume that the dependence of P uptake on the plant P : N ratio is modeled as a function of the P : N ratio of both the plant and its leaves, thereby capturing the essential N-P interaction through stoichiometric regulation. This regulatory mechanism helps prevent excessive P uptake, which would constitute luxury consumption for the plant (Schachtman et al., 1998). Similarly, if nitrogen uptake exceeds the plant's requirements, it also constitutes luxury consumption. Therefore, to avoid luxury absorption and nutrient accumulation, the uptake of N (or P) by roots needs to be regulated based on the N : C (or P : N) ratios within plant tissues (Goll et al., 2017). The maximum uptake occurs when the leaf N : C (or P : N) ratio is equal to the minimum leaf N : C (or P : N) ratio, which is calculated using a minimum function:

$$\begin{array}{}\text{(27)}& {\displaystyle }{f}_{\mathrm{P},\mathrm{ratio}}=min\left(max\left({\displaystyle \frac{{\mathrm{pn}}_{\mathrm{plant}}-{\mathrm{pn}}_{\mathrm{leaf},max}}{{\mathrm{pn}}_{\mathrm{leaf},min}-{\mathrm{pn}}_{\mathrm{leaf},max}}},\mathrm{0.0}\right),\mathrm{1.0}\right),\text{(28)}& {\displaystyle }{f}_{\mathrm{N},\mathrm{ratio}}=min\left(max\left({\displaystyle \frac{{\mathrm{nc}}_{\mathrm{plant}}-{\mathrm{nc}}_{\mathrm{leaf},max}}{{\mathrm{nc}}_{\mathrm{leaf},min}-{\mathrm{nc}}_{\mathrm{leaf},max}}},\mathrm{0.0}\right),\mathrm{1.0}\right),\end{array}$$

where ${\mathrm{pn}}_{\mathrm{leaf},max}$ and ${\mathrm{pn}}_{\mathrm{leaf},min}$ are prescribed maximum and minimum values of leaf P : N ratios, ${\mathrm{nc}}_{\mathrm{leaf},max}$ and ${\mathrm{nc}}_{\mathrm{leaf},min}$ are prescribed maximum and minimum values of leaf N : C ratios.

The root nutrient-uptake capacity function (uc_root,χ) incorporates both linear and Michaelis-Menten components to accurately represent the uptake process, considering the low-affinity and high-affinity transporter systems operating in parallel for a given solute concentration (Goll et al., 2017). Notably, the root uptake capacity for soil labile P (u_root,P) considers the replenishment of P from soil around the roots to root surfaces (Goll et al., 2017) rather than the total labile P in soil volume (Schachtman et al., 1998; Johnson et al., 2003). Hence, the calculation of root uptake capacity for N and P can be expressed as follows:

$$\begin{array}{}\text{(29)}& {\displaystyle }{u}_{\mathrm{root},\mathrm{P}}=v_{\max ,\mathrm{P}}\cdot P_{\mathrm{lab}}^{\prime }\left({\displaystyle \frac{k_{{\mathrm{P}}_{m\mathrm{1}}}}{c_k}}+{\displaystyle \frac{\mathrm{1}}{P_{\mathrm{lab}}+c_k{k}_{{\mathrm{P}}_{m\mathrm{2}}}}}\right),\text{(30)}& {\displaystyle }{u}_{\mathrm{root},\mathrm{N}}=v_{\max ,\mathrm{N}}\cdot N_{min}\left({\displaystyle \frac{k_{{\mathrm{N}}_{m\mathrm{1}}}}{c_k}}+{\displaystyle \frac{\mathrm{1}}{N_{min}+c_k{k}_{{\mathrm{N}}_{m\mathrm{2}}}}}\right),\end{array}$$

where $v_{max,\mathit{\chi }}$ is the maximum uptake capacity (Table 4). N_min and P_lab is the soil mineral N pool and labile P pool. $P_{\mathrm{lab}}^{\prime }$ represents the dissolved labile P concentration at the root surface and depends on the diffusion of soil labile P from the soil surrounding the roots to the root surface (Table 3; Eq. 54). c_k is a unit conversion factor using the soil-type specific parameter for soil moisture content at saturation as an approximation of pore space following Smith et al. (2014). $k_{{\mathit{\chi }}_{m\mathrm{1}}}$ was chosen to match the observed rate of increase in overall P uptake at high dissolved labile P concentration (low-affinity transporter), and $k_{x_{m\mathrm{2}}}$ is a parameter for Michaelis-Menten constants, dissolved phosphorus concentration at which uptake equals $\frac{{\mathit{\nu }}_{max}}{\mathrm{2}}$.

Mineralization and immobilization processes occur concurrently. Nutrient mineralization fluxes are estimated from the decomposition of litter and soil organic matter, assuming that C, N, and P mineralize at similar rates (Wang et al., 2010; Yang et al., 2014). The mineralization rate is determined by multiplying the litter and soil C pool turnover fluxes by the nutrient-to-carbon ratio. This can be mathematically represented by the following equations:

$$\begin{array}{}\text{(31)}& {\displaystyle }{F}_{\mathit{\chi },min,j}=Q_{\mathrm{C}}\left(t\right){\mathit{\tau }}_j\mathit{\xi }\left(t\right)R_{\mathit{\chi },j}^{-\mathrm{1}},\text{(32)}& {\displaystyle }{F}_{\mathit{\chi },min,m}=Q_{\mathrm{C}}\left(t\right){\mathit{\tau }}_m\mathit{\xi }\left(t\right)R_{\mathit{\chi },m}^{-\mathrm{1}},\end{array}$$

where Q_C(t)τξ(t) estimates the carbon decomposition rate under environmental stress for litter or soil pool. The total mineralization ($F_{\mathit{\chi },min,\mathrm{total}}$) is estimated as the sum of mineralization rate for each pool, which can be expressed as follows:

$$\begin{array}{}\text{(33)}& F_{\mathit{\chi },min,\mathrm{total}}=\sum _j{F}_{\mathit{\chi },min,j}+\sum _m{F}_{\mathit{\chi },min,m}.\end{array}$$

Nutrients are immobilized during the decomposition process of litter and SOM, ultimately entering the SOM pools. Consequently, only three SOM pools can be the receiving pools. The dependency of immobilization rates on the ratios of the receiving pools, under the assumption of approximately constant stoichiometry ratios of SOM pools (Tian et al., 2010; McGroddy et al., 2004), is described as:

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{F}_{\mathit{\chi },\mathrm{imm},m}& =\sum _{jm}{f}_{\mathrm{L}\mathrm{2}\mathrm{S},jm}\mathit{\xi }\left(t\right){\mathit{\tau }}_j{X}_j\left(t\right)R_{\mathit{\chi },m}^{-\mathrm{1}}\\ & +\sum _{mm}{f}_{\mathrm{S}\mathrm{2}\mathrm{S},mm}\mathit{\xi }\left(t\right){\mathit{\tau }}_m{X}_m\left(t\right)R_{\mathit{\chi },m}^{-\mathrm{1}},\end{array}\end{array}$$

where $R_{\mathit{\chi },m}^{-\mathrm{1}}$ represent the N : C ratio (χ = N) or the P : C ratio (χ = P) of the existing SOM. The total amount of immobilization is then calculated as follows:

$$\begin{array}{}\text{(35)}& F_{\mathit{\chi },\mathrm{imm},\mathrm{total}}=\sum _m{F}_{\mathit{\chi },\mathrm{imm},m}.\end{array}$$

Therefore, the net nutrient mineralization is calculated by the difference of total mineralization and total immobilization:

$$\begin{array}{}\text{(36)}& F_{\mathit{\chi },\mathrm{net}}=F_{\mathit{\chi },min,\mathrm{total}}-F_{\mathit{\chi },\mathrm{imm},\mathrm{total}}.\end{array}$$

When net mineralization is negative, the decomposition rate is limited by nutrient availability, L_NP. Since the N : C ratio of the soil pool is higher than that of the litter pool, microbes extract inorganic nitrogen from the soil mineral N pool, leading to negative net mineralization and a L_NP. value less than one. A similar approach has been applied in the CASA-CNP model (Wang et al., 2007).

Plants and microorganisms utilize dissolved inorganic N and P from the soil to fulfill their growth requirements (Vitousek et al., 2010). We assume microbial processes modulate nutrient availability for plants (Jiang et al., 2024b; Pellitier and Jackson, 2023; Jonasson et al., 1999), i.e., the nutrient limitation on plant growth will be alleviated if the net mineralization is positive. Furthermore, the competition between plants and microorganisms for nutrients can be simplified by emphasizing the sequence of immobilization and plant uptake (Achat et al., 2016). In the TECO-CNP model, immobilization takes precedence in nutrient access through the decomposition of litter and soil organic matter. A similar method was used in many models, e.g., models of the CENTURY family (e.g., Parton et al., 1988, 1998); O-CN (Zaehle and Friend 2010); ORCHIDEE (revision 4520; Goll et al., 2017). This also aligns with recent findings regarding the competition between plants and microbes under elevated CO₂ (Ben Keane et al., 2023). Specifically, in the acidic grassland, aboveground productivity and P uptake declined by 11 % and 20 %, respectively, while P immobilization into microbial biomass increased by 36 %.

2.2.3 Distinct processes in N and P cycle

The dynamic of soil inorganic N (N_min) is described as:

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\mathrm{d}t}{\mathrm{d}}}{N}_{min}=F_{\mathrm{fix}}+F_{\mathrm{N},\mathrm{fert}}+F_{\mathrm{N},\mathrm{dep}}-F_{\mathrm{N},\mathrm{leach}}-F_{\mathrm{N},\mathrm{gas}},\end{array}$$

where F_fix, F_N,fert and F_N,dep represent the biological N₂-fixation, atmospheric N deposition, and biological N fixation (Tables 2, 3). F_N,leach and F_N,gas represent the N leaching and gaseous N loss.

Biological nitrogen fixation, a dominant source of new nitrogen in terrestrial ecosystems (Chapin et al., 2011; Vitousek et al., 2013), is performed by N₂-fixing symbionts in plant roots (i.e., symbiotic N₂-fixation; Vitousek et al., 2002; Augusto et al., 2013). This process enhances nitrogen availability when carbon is sufficient for additional nutrient acquisition (Fisher et al., 2010), which is given by:

$$\begin{array}{}\text{(38)}& F_{\mathrm{fix}}=v_{\mathrm{fix}}\cdot f_{\mathrm{nsc}}\cdot \mathrm{NSC}\cdot f_{\mathrm{N}},\end{array}$$

where v_fix=0.00167 (g N g C⁻¹ m⁻² h⁻¹) is the maximum N fixation rate. v_fix is chosen based on estimates ranging from 58 Tg N yr⁻¹ (Vitousek et al., 2013) to 100 Tg N yr⁻¹ (Wiltshire et al., 2021) for a global NPP of 60 Pg C yr⁻¹. The term f_nsc⋅NSC represents the limitation of NSC on nitrogen fixation, implicitly capturing the carbon constraint on this process (Chou et al., 2018; Taylor and Menge, 2021). To prevent unrealistic nitrogen fixation, a scaling function (f_N) is applied, as nitrogen fixation is an energy-intensive process (Gutschick, 1981; Goll et al., 2017). The f_N is calculated as:

$$\begin{array}{}\text{(39)}& f_{\mathrm{N}}=\left\{\begin{array}{ll}\frac{N_{max}-N_{min}}{N_{max}}& N_{max}<N_{min}\\ \mathrm{0}& \text{otherwise}\end{array}\right..\end{array}$$

The carbon cost for biological N fixation is calculated by a function of soil temperature (T_soil) with the observed carbon cost range (Fisher et al., 2010):

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{C}_{\mathrm{fix}}& =-\mathrm{6.25}\cdot \left(\exp \right(-\mathrm{3.62}+\mathrm{0.27}\cdot T_{\mathrm{soil}}\\ & \cdot (\mathrm{1}-\mathrm{0.5}({\displaystyle \frac{T_{\mathrm{soil}}}{\mathrm{25.15}}}\left)\right))-\mathrm{2}).\end{array}\end{array}$$

Nitrogen loss occurs in two pathways: gaseous loss (FN_gas), and leaching (F_N,leach). Losses from denitrification and volatilization are not distinguished separately. Both are proportional to the availability of soil mineral N (N_min). The expression of N leaching is:

$$\begin{array}{}\text{(41)}& F_{\mathrm{N},\mathrm{leach}}=f_{\mathrm{N},\mathrm{leach}}{\displaystyle \frac{V_{\mathrm{runoff}}}{D_{\mathrm{soil}}}}{N}_{min},\end{array}$$

where $f_{\mathrm{N},\mathrm{loss}}=\mathrm{0.001}$ and $f_{\mathrm{N},\mathrm{leach}}=\mathrm{0.5}$, V_runoff is the soil surface runoff and D_soil is the soil depth. Moreover, the gaseous loss is dependent on the soil temperature and soil mineral N. The equation is:

$$\begin{array}{}\text{(42)}& F_{\mathrm{gas}}=f_{\mathrm{N},\mathrm{loss}}{e}^{\frac{(T_{\mathrm{soil}}-\mathrm{25})}{\mathrm{10}}}{N}_{min}.\end{array}$$

The specific processes of the P cycle include biochemical mineralization, weathering, the dynamics of different inorganic soil P components, and the diffusion pathways of soil labile P. In addition to biological mineralization, organic P can be mineralized through direct cleavage by extracellular enzymes produced by plant roots and other organisms (McGill and Cole, 1981). This process decouples the P cycle from the C and N cycles, serving as an adaptive mechanism that can be enhanced under P-limited conditions (Lambers et al., 2006). This decoupling allows for phosphorus acquisition from organic matter without releasing carbon dioxide. We consider this process an N-consuming one, aiming to represent the chemical characteristic that phosphatases are N-rich enzymes and their production in plants can be N-limited (Treseder and Vitousek, 2001; Wassen et al., 2013). The biochemical mineralization of P can be expressed by:

$$\begin{array}{}\text{(43)}& {\mathrm{FP}}_{\mathrm{biomin},m}={\displaystyle \frac{{\mathit{\upsilon }}_{max}\left({\mathit{\lambda }}_{\mathrm{up}}-{\mathit{\lambda }}_{\mathrm{ptase}}\right)}{{\mathit{\lambda }}_{\mathrm{up}}-{\mathit{\lambda }}_{\mathrm{ptase}}+{\mathit{\kappa }}_m}}\sum _m{K}_m{Q}_{\mathrm{P},m},\end{array}$$

where υ_max is maximal specific rate of biochemical P mineralization. λ_up is N cost of plant root P uptake. λ_ptase is the N cost of phosphatase production, κ_m is the Michaelis-Menten constant for biochemical P mineralization. K_m and Q_P,m represent turnover rate and phosphorus pool size of slow (m=8) and passive pools (m=9). Phosphatase production is activated when λ_up>λ_ptase, reflecting N regulation of P acquisition strategy by plants, similar to the cost-benefit approach established in existing coupled carbon-nutrient models (Wang et al., 2007; Houlton et al., 2008). This modeling approach aligns with findings that nitrogen addition significantly increases phosphatase activity (Schleuss et al., 2020; Marklein and Houlton, 2012), potentially through enhanced phosphorus limitation and elevated plant nitrogen status, which favor investment in the phosphatase enzyme. While direct field quantification of biochemical mineralization rates is not yet possible, this mechanistic representation becomes particularly important for predicting ecosystem responses to elevated CO₂ and enhanced N deposition, where enhanced biochemical mineralization of soil organic P may facilitate additional plant growth (Jiang et al., 2024a).

The soil P loss from soil organic P pools can be simulated by the following equations:

$$\begin{array}{}\text{(44)}& F_{\mathrm{P},\mathrm{out},m}\left(t\right)=Q_{\mathrm{C},m}\cdot {\mathit{\tau }}_m\cdot R_{\mathit{\chi },m}^{-\mathrm{1}}+{\mathrm{FP}}_{\mathrm{biomin},m}.\end{array}$$

The term FP_biomin,m equals 0 when m=7 as organic P losses through biochemical mineralization only occur in two soil pools with slow turnover rates (slow and passive pools; Wang et al., 2010).

The external phosphorus input (F_P,in) is modeled as:

$$\begin{array}{}\text{(45)}& F_{\mathrm{P},\mathrm{in}}=F_{\mathrm{wea}}+F_{\mathrm{P},\mathrm{fert}}+F_{\mathrm{P},\mathrm{dep}},\end{array}$$

where F_wea, F_P,fert, and F_P,dep represent phosphorus input rates from weathering, fertilization, and deposition. Based on the soil texture at the Tiantong site (Song and Wang, 1995), the weathering rate is set to 0.005 (g P m⁻² yr⁻¹) (Wang et al., 2010). The deposition rate of phosphorus has been set to 0.06 (g P m⁻² yr⁻¹) (Zhu et al., 2016a).

Labile phosphorus (P_lab) can be directly utilized by plants or microorganisms and adsorbed onto soil particles, organic matter, and other minerals as adsorbed phosphorus (P_S) (Vitousek et al., 2010). The assumption is made that the rapid equilibration of P_lab with P_S occurs within a timestep of less than one hour (Olander and Vitousek, 2005). For the 1 h time step used in our study, we therefore assume that P_lab and P_S are in a state of equilibrium. The equilibrium assumption is applied extensively (e.g., Wang et al., 2007; Yang et al., 2014). The relationship between them is described by a Langmuir isotherm (Barrow, 2008):

$$\begin{array}{}\text{(46)}& P_{\mathrm{S}}={\displaystyle \frac{S_{max}{P}_{\mathrm{lab}}}{K_s+P_{\mathrm{lab}}}},\end{array}$$

where S_max is the maximum amount of sorbed P in the soil, and K_s is the empirical constant representing the tendency of soil labile P to be sorbed. S_max and K_s is set as 133 and 64 (Wang et al., 2010), respectively, according to the soil sorption capacity and substrate age (Olander and Vitousek, 2005) at the Tiantong site. The differential form of Eq. (46) is:

$$\begin{array}{}\text{(47)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{S}}}{\mathrm{d}t}}={\displaystyle \frac{S_{max}{P}_{\mathrm{lab}}}{(K_s+P_{\mathrm{lab}}{)}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{lab}}}{\mathrm{d}t}}.\end{array}$$

Assuming equilibrium between P_lab and P_S, we can model the simultaneous changes in P_lab and P_S as follows:

$$\begin{array}{}\text{(48)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left(P_{\mathrm{S}}+P_{\mathrm{lab}}\right)}{\mathrm{d}t}}& =F_{\mathrm{P},\mathrm{net}}+F_{\mathrm{P},\mathrm{in}}+F_{\mathrm{P},\mathrm{biomin}}-F_{\mathrm{up},\mathrm{P}}\\ & -F_{\mathrm{P},\mathrm{leach}}-{\mathit{\nu }}_m{P}_{\mathrm{S}},\end{array}\text{(49)}& {\displaystyle }{F}_{\mathrm{P},\mathrm{net}}=F_{\mathrm{P},\min ,\mathrm{total}}-F_{\mathrm{P},\mathrm{imm},\mathrm{total}},\end{array}$$

where F_P,net is the net mineralization of litter and soil phosphorus pool, F_P,biomin is the P flux from biochemical mineralization, U_P represents the plant uptake of P, F_P,leach represents the loss of labile P from leaching (Eq. 53), and ν_m is the rate constant for the transformation of sorbed P to secondary P. Based on Eqs. (48) and (49), the dynamics of labile phosphorus can be expressed as follows:

$$\begin{array}{}\text{(50)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{lab}}}{\mathrm{d}t}}& =(F_{\mathrm{P},\mathrm{net}}+F_{\mathrm{P},\mathrm{in}}+{\mathrm{FP}}_{\mathrm{biomin}}-F_{\mathrm{up},\mathrm{P}}-F_{\mathrm{P},\mathrm{loss}}\\ & -{\mathit{\nu }}_m{P}_{\mathrm{S}}){\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{S_{max}{P}_{\mathrm{lab}}}{{\left(K_s+P_{\mathrm{lab}}\right)}^{\mathrm{2}}}}},\end{array}\end{array}$$

The use of solution P would be theoretically more appropriate, as previous studies have shown that models operating at very fine temporal resolutions (hourly or finer) may require distinction between labile and solution phosphorus pools (Reed et al., 2015; Yang et al., 2013). However, implementing this simulation approach is currently not feasible due to limited data availability. Some synthesis studies (Yang et al., 2013; Hou et al., 2018) have indicated that most experimental measurements report total labile P, without separating it into distinct fractions. Additionally, previous studies have demonstrated strong correlations between these fractions. For example, strip- and NaHCO₃-extracted inorganic P are positively correlated and exhibit similar temporal patterns during experimental periods (Hou et al., 2019). Due to these reasons, we adopted labile P as our primary plant-available phosphorus pool in our model.

Secondary mineral phosphorus (P_SS) can be dissolved and enter the labile P pool or encapsulated by iron oxides to form closed-P (P_O; Walker and Syers, 1976; Vitousek et al., 2010). The dynamics of P_SS and P_O can be modeled as:

$$\begin{array}{}\text{(51)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{SS}}}{\mathrm{d}t}}={\mathit{\nu }}_m{P}_{\mathrm{S}}-{\mathit{\nu }}_{\mathrm{dis}}{P}_{\mathrm{SS}}-{\mathit{\nu }}_{\mathrm{o}}{P}_{\mathrm{SS}},\text{(52)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{O}}}{\mathrm{d}t}}={\mathit{\nu }}_{\mathrm{o}}{P}_{\mathrm{SS}}-{\mathit{\nu }}_{\mathrm{re}}{P}_{\mathrm{O}},\end{array}$$

where ν_dis and ν_o is the rate constant for the conversion of secondary P to labile and sorbed P, and occluded P, respectively. ν_re is the rate constant for occluded P re-entering the cycle as bioavailable phosphorus, indicating that occluded phosphorus can transition back into available forms (Huang et al., 2014; Schubert et al., 2020). In this study, we assume that the formation of occluded P pool and loss of occluded P can be considered negligible within the short timescale of simulations (Weihrauch and Opp, 2018). P leaching from the soil inorganic labile pool is proportional to the availability of soil labile P. Description of P leaching below:

$$\begin{array}{}\text{(53)}& F_{\mathrm{P},\mathrm{leach}}=f_{P,\mathrm{leach}}{\displaystyle \frac{V_{\mathrm{runoff}}}{D_{\mathrm{soil}}}}{P}_{\mathrm{lab}},\end{array}$$

where V_runoff is the value of runoff, D_soil is the soil depth. f_P,leach is an empirical parameter for P leaching, representing the fraction of soil mineral P for leaching.

Notably, due to the low mobility of phosphorus in the soil (Vitousek et al., 2010), the actual P concentration that roots can absorb depends on the diffusion of P from the surrounding soil to the root surface ($P_{\mathrm{lab}}^{\prime }$). This finding is consistent with the experimental evidence that roots primarily acquire most inorganic phosphorus through diffusion along concentration gradients (Laliberté et al., 2015). Thus, the root uptake capacity for soil labile P (u_root,P) considers the replenishment of P from soil around the roots to root surfaces (Schachtman et al., 1998) rather than the total labile P in soil volume (Johnson et al., 2003). Thus, the root surface P concentration is calculated by the following equation:

$$\begin{array}{}\text{(54)}& P_{\mathrm{lab}}^{\prime }=a_{\mathrm{root}}\cdot {\displaystyle \frac{P_{\mathrm{lab}}}{\mathrm{\Theta }}},\end{array}$$

where Θ is the volumetric soil water content and a_root representing the fraction of the reduction in P concentration surrounding the roots relative to the initial concentration. a_root is updated after plant uptake as:

$$\begin{array}{}\text{(55)}& {\displaystyle \frac{\mathrm{d}{a}_{\mathrm{root}}}{\mathrm{d}t}}={\displaystyle \frac{{\mathrm{FP}}_{\mathrm{diff}}-F_{\mathrm{up},\mathrm{P}}}{P_{\mathrm{lab}}}},\end{array}$$

where FP_diff is the diffusion of P from the surroundings to the root surface, which is the function of the permeability of the soil to P (K) and the difference in the P concentrations between the soil solution at the root surface and the labile P in the surrounding soil volume outside the diffusive zone around the root (ΔP_lab)

$$\begin{array}{}\text{(56)}& {\mathrm{FP}}_{\mathrm{diff}}=-K\cdot \mathrm{\Delta }{P}_{\mathrm{lab}}.\end{array}$$

ΔP_lab can be described as:

$$\begin{array}{}\text{(57)}& \mathrm{\Delta }{P}_{\mathrm{lab}}=(a_{\mathrm{root}}-\mathrm{1}){\displaystyle \frac{P_{\mathrm{lab}}}{\mathrm{\Theta }}}.\end{array}$$

The K has been calculated analogously to the diffusion coefficient of phosphorus in soils following Barraclough and Tinker (1981), which accounts for the increased path length in soil using a tortuosity factor (f_t), and it is a broken-line function of the volumetric soil water content (Θ). The K and f_t can be calculated based on the following equations:

$$\begin{array}{}\text{(58)}& {\displaystyle }K=K_{\mathrm{0}}{c}_{\mathrm{\Theta }}\mathrm{\Theta }tf{\displaystyle \frac{\mathrm{1}}{r_{\mathrm{diff}}}},\text{(59)}& {\displaystyle }{f}_{\mathrm{t}}=\left\{\begin{array}{ll}{f}_{\mathrm{1}}\mathrm{\Theta }+f_{\mathrm{2}}& \text{for}\mathrm{\Theta }\ge {\mathrm{\Theta }}_{\mathrm{1}}\\ \frac{\mathrm{\Theta }\left(f_{\mathrm{1}}\mathrm{\Theta }+f_{\mathrm{2}}\right)}{{\mathrm{\Theta }}_{\mathrm{1}}}& \text{otherwise}\end{array}\right.,\end{array}$$

where Θ₁ is soil water content at which the two functions intersect according to Barraclough and Tinker (1981), f₁ and f₂ are empirical parameters (Barraclough and Tinker, 1981), D₀ is diffusion coefficient in free water, c_Θ is a unit conversion factor, r_diff is diffusion path, which can be calculated from the function of root length density (RLD, Bonan et al., 2014):

$$\begin{array}{}\text{(60)}& r_{\mathrm{diff}}=min\left(\mathrm{0.1},{\left(\mathit{\pi }\mathrm{RLD}\right)}^{\mathrm{0.5}}\right).\end{array}$$

We assume that the diffusion path can be approximated as half the average distance between roots. We limit the diffusion path length to 0.1 m because the influence of active P uptake by roots on soil P concentrations is negligible beyond a distance of 10 cm (Li et al., 1991). RLD is given by:

$$\begin{array}{}\text{(61)}& \mathrm{RLD}={\displaystyle \frac{B_{\mathrm{root}}^{\ast }}{r_{\mathrm{d}}\mathit{\pi }{r}_{\mathrm{r}}^{\mathrm{2}}}},\end{array}$$

where r_d is the root-specific density and $\mathit{\pi }{r}_{\mathrm{r}}^{\mathrm{2}}$ is the cross-sectional area calculated from the fine root radius, r_r, and $B_{\mathrm{root}}^{\ast }$ is the root biomass density per unit soil volume.

2.3 Model validation

2.3.1 Study site

The tension between high carbon sink capacity and nutrient limitations in subtropical forests warrants detailed investigation to understand the role of nutrients in carbon cycling processes in these regions. To this end, we selected a mature subtropical evergreen broadleaf forest in eastern China, located at the Zhejiang Tiantong Forest Ecosystem National Observation and Research Station (Tiantong, 29°48^′ N, 121°47^′ E, Fig. 3) for the newly developed model. The Tiantong forest has been preserved free from human disturbance since the mid-twentieth century. The average annual temperature in Tiantong is 17 °C, and the annual precipitation is 1600 mm (Cui et al., 2022). The soil type is mainly mountainous yellow-red soil, with the parent material primarily composed of Mesozoic sedimentary rocks, acidic igneous rocks, and residual weathering products of granite (Song and Wang, 1995).

Figure 3 Schematic diagram of the observation system at Tiantong subtropical evergreen forest (29°5^′ N, 121°5^′ E). The system comprises: (1) a forest dynamic plot for monitoring ecosystem state variables, including stoichiometric ratios, plant traits, and C, N, P pools and fluxes. These measurements were conducted in a 5 ha subplot of the whole plot. The asterisk (*) indicates manually measured periodic fluxes. And (2) an eddy covariance (EC) flux tower providing half-hourly NEE measurements, from which GPP and ER were derived. These observations were used for parameterizing and evaluating the TECO-CNP model. Detailed measurement protocols are described in the Methods section, and specific variable applications are listed in Tables S1–S3. The 3D visualization of the study site was created in Blender (v4.2.1) using topographic data sourced from OpenTopography (https://opentopography.org, last access: 15 March 2025).

Research conducted at this site revealed the dominant role of soil phosphorus in driving variations in plant functional traits (Cui et al., 2022), indicating phosphorus deficiency in this mature forest. Consequently, this phosphorus-limited mature subtropical forest, with abundant field observations, can contribute to the development of carbon-nutrient coupling models and further explore phosphorus-limited carbon cycling processes within the ecosystem through the integration of modeling and experiments.

2.3.2 Data collection and site parameterization

The data used for model calibration and validation were primarily derived from our field measurements and literature focusing on the same site (Fig. 3). The forcing data for TECO are collected at 1 h intervals from field-based meteorological observations at the study site, include precipitation (mm), relative humidity ( %), air and soil temperatures (°C), vapor pressure deficit (Pa), wind speed (m s⁻¹), and photosynthetically active radiation (µmol m⁻² s⁻¹). Forcing data from 2001 was used for model spin-up.

Site-specific parameters that can be empirically measured are derived from field observations at the study site, including both our measurements and values reported in previous studies. Plant traits, including specific leaf area (SLA, cm² g⁻¹), leaf area index (LAI, m² m⁻²), plant height (H, m), maximum rate of carboxylation (V_cmax, µmol m⁻² s⁻¹), maximum rate of electron transport (J_max, µmol m⁻² s⁻¹) and leaf P concentration (Leaf P, g m⁻²), were measured at the species level in the forest dynamic plot, and scale up to community-level traits using the community-weight mean method (for detailed sampling methods, refer to Cui et al., 2022). Plant stoichiometry ratios were derived from area-based C, N, and P concentrations from both our measurements and previous studies at Tiantong (Zhou, 2020). N and P resorption efficiencies were determined based on dominant species (i.e., Schima superba, Lithocarpus glaber) at the Tiantong site (Xu et al., 2020). The observed data used for model parameterization are presented in Tables S1–S3.

Parameters not readily accessible through field measurements are estimated using standard procedures that have been extensively validated in other modeling studies, with appropriate selection based on the characteristics of Tiantong. For example, Tiantong forest soils are classified as Ultisols (Song and Wang, 1995), which directly informed our selection of phosphorus weathering rates and inorganic P dynamics parameters (e.g., K_s and S_max, Table 4). Similarly, the subtropical evergreen broadleaf forest vegetation type guided our parameterization of phosphorus mineralization (Wang et al., 2010) and allocation processes (Arora and Boer, 2005). External inputs of N and P, including deposition and weathering, were assumed to occur at constant rates. Deposition rates for N and P were prescribed based on the observed range (Zhu et al., 2016a). Specific parameterization settings are described in Table 4, along with the accompanying process descriptions.

The observed pool sizes and fluxes primarily serve as a basis for model evaluation and as references for model initialization. Soil inorganic pools of mineral N and labile P were determined from 0–10 cm soil samples collected in 2023 from a nearby forest stand of similar age (∼ 200 years) dominated by the same species (Schima superba and Castanopsis fargesii) as the Tiantong forest dynamic plot. Labile P is the soil inorganic phosphorus fraction that can be extracted by resin and NaHCO₃. Sampling employed a five-point design with three replicates per point. Additionally, secondary P and occluded P refer to the measured moderately labile inorganic phosphorus (extracted by NaOH) and moderately stable inorganic phosphorus (extracted by HCl), respectively, at the Tiantong site (Wang, 2022).

Soil C, N, and P were measured using systematic sampling across 185 grid points (each 20 × 20 m) within the permanent Tiantong forest dynamic plot (Fig. 3). At each grid point, soil samples were collected at three depth intervals (0–20, 20–40, 40–60 cm) in 2017 using a 5 cm diameter auger, with three replicates per depth. Additionally, observed plant pools and fluxes, including fluxes from plant to litter and soil respiration, used for model evaluation and their sources are listed in Tables S1–S3. Quality-controlled hourly eddy covariance measurements of gross primary productivity (GPP), ecosystem respiration (ER), and net ecosystem exchange (NEE) were obtained from the on-site flux tower for the year 2021.

All model configurations used identical site-specific parameter sets obtained according to the methods described above. Although a previous study has highlighted the necessity for model-specific reparameterization (Wang et al., 2022), we adopted a consistent parameterization approach across all configurations. This follows common practice in land surface model development studies, where uniform parameterization is essential for isolating the effects of different nutrient coupling schemes.

2.3.3 Data assimilation

We specifically optimized carbon-related parameters for the CNP configuration only, utilizing GPP, ER, and NEE data from 2021 at the study site, to evaluate the effectiveness of the CNP structure coupled with a data assimilation algorithm. Based on the initial carbon pool sizes from the spin-up process, a preliminary sensitivity analysis was first conducted to support the selection of target parameters for data assimilation. We focused on parameters that determine carbon input and retention (Table 6), including SLA, V_cmax, and temperature sensitivity (Q₁₀), which showed high sensitivity in the analysis (Table S6). Additionally, our parameter selection strategy included all carbon pool turnover parameters (T₁–T₉), as these govern carbon residence times and are crucial for matching observed pool dynamics, regardless of their sensitivity indices. The prior range of parameters was prescribed according to the situ measurement or assumed as the range of the distribution to be [${\mathit{\theta }}_{\mathrm{0}}/\mathrm{3},\mathrm{3}{\mathit{\theta }}_{\mathrm{0}}$], where θ₀ is the default value. Using the Bayesian probabilistic inversion approach, we estimated the posterior distribution of model parameters based on prior knowledge of the parameters.

Bayesian probabilistic inversion approach is based on Bayes' theorem:

$$\begin{array}{}\text{(62)}& p\left(\mathit{\theta }\mathrm{|}Z\right)\propto {\displaystyle \frac{p\left(Z\mathrm{|}\mathit{\theta }\right)\times p\left(\mathit{\theta }\right)}{p\left(Z\right)}},\end{array}$$

where p(θ|Z) is the posterior distribution of the parameters θ given the observations Z. Here, we assume that the prior knowledge of parameter distribution p(θ) is uniformly distributed. p(Z|θ) is the likelihood function for a parameter set calculated with the assumption that each parameter is independent from all other parameters and has a normal distribution with a zero mean:

$$\begin{array}{}\text{(63)}& p\left(Z\mathrm{|}\mathit{\theta }\right)\propto \exp \left\{-\sum _{t\in Z_i}{\displaystyle \frac{{\left[Z_i\left(t\right)-X\left(t\right)\right]}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\left(t\right)}}\right\},\end{array}$$

where Z_i(t) is the observations of carbon fluxes at time t, X(t) is the simulated corresponding variable, and σ(t) is the standard deviation of the observation set.

Posterior probability distributions of the parameters were obtained using a Metropolis-Hastings (M-H) algorithm within the Markov Chain Monte Carlo (MCMC) framework. The posterior parameter distribution represents our updated knowledge about parameter values after incorporating observational data through Bayesian inference, quantifying both the most likely parameter estimates and their associated uncertainties. The detailed description of the M-H algorithm can be found in Xu et al. (2006). In brief, the M-H algorithm consists of iterations that alternate between a proposing step and a moving step. In the proposing step, a new parameter set θ^new is proposed based on the previously accepted parameter set θ^old and a proposal distribution ($r\times ({\mathit{\theta }}_{max}-{\mathit{\theta }}_{min})/D$):

$$\begin{array}{}\text{(64)}& {\mathit{\theta }}^{\mathrm{new}}={\mathit{\theta }}^{\mathrm{old}}+r\times ({\mathit{\theta }}_{max}-{\mathit{\theta }}_{min})/D,\end{array}$$

where θ_max and θ_min corresponding to the upper and lower values of prescribed ranges, r is a random variable between −0.5 and 0.5, and D is used to control the proposed step size and was set to 5 (Xu et al., 2006). The new set of parameter values would be accepted when $\frac{p\left({\mathit{\theta }}^{\mathrm{new}}\mathrm{|}Z\right)}{p\left({\mathit{\theta }}^{\mathrm{old}}\mathrm{|}Z\right)}$ is equal or greater than a uniform random number from 0 to 1 (Xu et al., 2006).

We get 10 000 accepted samples from the MCMC chain. The first 5000 accepted samples were discarded, considering the burn-in period. We randomly selected 1000 parameter sets from the accepted space to run the simulations in 2021. The mean and maximum likelihood estimations are calculated to compare the parameters.

2.3.4 Model performance evaluation

The state variables estimations from three nutrient coupling configurations of TECO-CNP: (1) carbon-only (C-only), (2) carbon-nitrogen coupled (CN), and (3) carbon-nitrogen-phosphorus coupled (CNP) are evaluated against observations. Model initialization involved a spin-up process using 2001 meteorological forcing data until a quasi-equilibrium state was reached, defined as inter-annual variations of less than 0.05 g C m⁻² yr⁻¹ in the slowest pools. Following initialization, we conducted transient simulations from 2002 to 2021 using the tuned parameter set. To evaluate model performance, we compared pool sizes from different nutrient coupling configurations (C-only, CN, and CNP) in 2021 with observed data (Tables S1–S3), assuming that our mature forest study site was at a quasi-steady state, where interannual changes in major pool sizes were negligible. The configuration that produced pool sizes closest to observations was selected to determine the initial state for subsequent simulations. Model performance was further evaluated by comparing simulated carbon fluxes in 2021 against observational data using both manually tuned and optimized parameters. The model evaluation metrics for carbon fluxes included the Root Mean Square Error (RMSE) and concordance correlation coefficient (CC), which quantify the absolute errors and the agreement between simulated and observed values. All statistical analyses and data visualizations were implemented in R (version 4.3.1).

3 Results & Discussion

3.1 Evaluate the carbon-nutrient configurations

3.1.1 Carbon cycle

The CNP configuration accurately reproduced carbon pool sizes across ecosystem components. In contrast, the C-only and CN configurations tended to overestimate these pools (Figs. 4, 5a). In this P-limited site, the introduction of phosphorus limitations in CNP configurations progressively reduced carbon pool sizes compared to the C-only and CN configurations (Fig. 4a). This reduction reflects a fundamental assumption in carbon-nutrient coupled models that nutrient availability constrains carbon sequestration (Wieder et al., 2015; Sun et al., 2017) through various physiological processes (Jiang et al., 2019). At the ecosystem level (Fig. 4b), the C-only and CN configurations substantially overestimated total carbon stocks by 73.7 % and 57.5 %, respectively. In contrast, the CNP configuration produced estimates that were much closer to the observed values, with only a slight overestimate of 1.9 %. The partitioning between plant and soil pools (Fig. 4b) showed that this overestimation occurred in both compartments, with the CNP configuration providing the closest match to observations.

Figure 4 Comparison of carbon pools among different nutrient coupling configurations. (a) Trajectories of ecosystem carbon pools during model spin-up for carbon-only (C-only), coupled carbon-nitrogen (CN), and coupled carbon-nitrogen-phosphorus (CNP) simulations. The ecosystem carbon pool comprises nine pools within the components of plant, litter, and soil organic matter. (b) Comparison of simulated and observed (OBS) carbon pools in plant biomass and soil organic matter. Plant carbon pools comprise leaf, wood, and root carbon (excluding reproductive organs due to data unavailability), and soil carbon pools include fast, slow, and passive soil organic carbon components. The error bar for observation represents the standard deviation of the sum of plant and soil pools.

Figure 5 Comparison of simulated and observed ecosystem pools across different nutrient coupling configurations (C-only, CN, and CNP). (a) Carbon pools in vegetation components (leaf, wood, root) and soil, with values for leaf and root scaled by 10¹. (b) Nitrogen pools in vegetation components, soil (scaled by 10⁻²), and mineral nitrogen (N_min). (c) Phosphorus pools in vegetation components, soil organic P (scaled by 10⁻²), labile P (P_lab), and sorbed P (P_S). Error bars on observed data (OBS) indicate standard deviations. Numbers in parentheses indicate scaling factors applied to improve visualization. For example, the soil P value marked with 10⁻² indicates that this value has been scaled down, and the actual value is $\mathrm{1.58}/{\mathrm{10}}^{-\mathrm{2}}=\mathrm{158}$ g P m⁻². Shaded areas indicate inorganic nutrient pools.

A more detailed examination of individual carbon pools (Fig. 5a) revealed that the overestimation was mainly contributed by wood and soil pools for C-only and CN configurations, which represent the major carbon stocks in the ecosystem. For plant components, wood carbon stocks were substantially overestimated by approximately 122.2 % and 89.6 % in the C-only and CN configurations, respectively. In contrast, the CNP configuration showed remarkable agreement with observations, with only a 5 % deviation. Leaf carbon pools showed similar patterns of overestimation (C-only: 82.7 %, CN: 59.1 %, CNP: 3.6 %). This improvement in leaf carbon estimation by CNP was further confirmed by better LAI prediction: the CNP configuration (3.94 m² m⁻²) showed only 5 % deviation from observations (3.75 ± 0.15 m² m⁻²), while C-only and CN configurations overestimated by 85 % and 61 %, respectively.

The observed reduction in LAI represents a decrease in photosynthetic capacity achieved through nutrient limitation of plant growth, which reduces the photosynthetic leaf area rather than directly affecting leaf-level photosynthetic physiological parameters. The relationships between leaf nitrogen and phosphorus concentrations and photosynthetic traits (e.g., V_cmax, J_max) are well established (Walker et al., 2014; Ellsworth et al., 2022) and have been incorporated into some land surface models (e.g., JULES-CNP). However, these large-scale emergent relationships significantly overestimated photosynthetic parameters at our study site (Table S5). At the same time, our site-specific dataset was insufficient to derive robust empirical relationships between nutrient concentrations and photosynthetic capacity. Future studies with more comprehensive site-level measurements could enhance this aspect of the model to represent nutrient-carbon interactions better.

Additionally, the CNP configuration better captured observed carbon fluxes compared to the C-only and CN configurations (Table 5). Although the total plant carbon litterfall rate was moderately overestimated by 22.7 %, this still reflects improved simulation of aboveground carbon dynamics and could be further refined by incorporating reproductive pool measurements in future studies.

Table 5 Observed and simulated carbon, nitrogen, and phosphorus fluxes with C-only, CN, and CNP configurations. The plant litterfall rate is the sum of the litterfall of leaves, woody, and reproductive parts.

^* Jiang et al. (2024b).

In contrast, root carbon pools showed an overestimation across all configurations, with CNP exhibiting the lowest bias (34.2 %) and falling within one standard deviation of the observed values (Table S1), while the C-only and CN configurations showed larger deviations (68.8 % and 65.1 %, respectively). The relatively higher root carbon estimation in CNP may be attributed to its dynamic allocation strategy, which preferentially allocates carbon to roots under nutrient-limited conditions. While our model successfully reproduced the enhanced belowground carbon allocation under nutrient limitation, consistent with experimental evidence (Wu et al., 2025; Gill and Finzi, 2016), the overestimated root carbon suggests additional constraints are needed. Indeed, the nutrient-dependent allocation scheme remains a significant source of uncertainty in terrestrial biosphere models (Zaehle et al., 2014; Jiang et al., 2024a). Although dynamic allocation schemes have been demonstrated to be significantly influenced by nutrient availability (Xia et al., 2023), explicit nutrient controls on allocation remain underrepresented in many ecosystem models (De Kauwe et al., 2014; but see Knox et al., 2024). Our model presents a practical approach for representing the nutrient regulation of carbon allocation processes. These results highlight the necessity of improved observational constraints on root turnover and carbon allocation patterns for more accurate process-based simulations.

For soil carbon pools, while C-only and CN configurations showed significant overestimations of 59.1 % and 52.1 %, respectively, the CNP configuration demonstrated the closest agreement with observations, with a slight overestimation of 1.06 %. Despite the considerable observational uncertainty in soil carbon stocks (Table S1), the substantial overestimation by C-only and CN configurations was clearly beyond the reasonable range. This distinct improvement in soil carbon estimation by CNP configuration suggests that proper representation of nutrient limitations is crucial for realistic soil carbon predictions (Cui et al., 2024; Wei et al., 2022; Achat et al., 2016). In conclusion, the CNP model consistently shows better alignment with observed carbon pools, particularly in reducing the systematic overestimation seen in the C-only and CN models.

3.1.2 N cycle

For nitrogen cycling properties, the CNP configuration exhibited superior performance in simulating nutrient pools compared to CN configurations (Fig. 5b). Regarding plant nitrogen pools, the CN configuration demonstrated substantial overestimations for leaf (59.2 %), woody tissue (89.9 %), and root N (55.9 %). In contrast, the CNP configuration showed markedly improved accuracy, with only slight overestimations of 3.3 %, 5.0 % for leaf and wood N, and 28.8 % for root N. The patterns of plant organic N across model configuration simulations were consistent with the carbon simulation results in both CN and CNP configurations, reflecting the constraints of plant tissue stoichiometry on coupled C-nutrient dynamics (Knox et al., 2024; Wang et al., 2010). For soil N pools, the CNP simulation (16.74 g N m⁻²) fell within the range of observed values (18.6 ± 5.5 g N m⁻²), whereas the CN configuration substantially overestimated soil N (28.75 g N m⁻²). The slight underestimation of soil N in CNP relative to observations may be attributed to the flexible soil C : N ratios, as these ratios can vary within specific ranges due to complex microbial processes and dynamics of organic matter decomposition (Tian et al., 2010, 2021). The introduction of P cycling into the model resulted in reduced carbon allocation to both plant and soil pools, which consequently led to proportional reductions in organic N pools compared to the CN configuration, ultimately capturing the observed N pools more accurately.

For soil mineral N content, the CN configuration underestimated soil mineral N content by 33.3 % despite simulating higher net N mineralization rates (Table 5). This discrepancy likely reflects the absence of phosphorus constraints in the CN model. While the CN model simulated higher net N mineralization than the CNP model (Table 5), this enhanced nitrogen input was offset by excessive plant N uptake. This is consistent with the substantial overestimation of plant carbon pools in the CN configuration (Fig. 5a) and the correspondingly lower soil mineral N reserves (Fig. 5b). In contrast, the CNP configuration showed a moderate overestimation (15.9 %) of soil mineral N content, demonstrating better agreement with observations compared to CN. The elevated soil mineral N levels in CNP could be attributed to the higher plant N litterfall rates (10.4 % above observed rates, Table 5), which compensated for the underestimated net N mineralization rates.

The incorporation of P cycling constraints in the CNP configuration substantially improved the simulation of N pools and fluxes compared to the CN configuration, demonstrating the importance of considering N-P interactions in ecosystem modeling. This improvement reflects the fundamental interconnectedness of nitrogen and phosphorus cycles, where phosphorus availability directly regulates plant nitrogen demand and uptake efficiency, while nitrogen status influences phosphorus acquisition strategies (Elser et al., 2007; Peñuelas et al., 2013). In our model, these interactions are primarily captured through the tight coupling between soil nutrient availability, plant stoichiometry, and plant growth processes, which prevents unrealistic carbon and nitrogen accumulation when phosphorus becomes limiting. Notably, our model has limitations in capturing the full complexity of N-P interactions, reflecting broader challenges in coupled CNP modeling (Achat et al., 2016). For example, the absence of linkages between nitrogen fixation processes and phosphatase enzyme activity (Batterman et al., 2018), as well as the simplified representation of plant-microbe competition for nutrients and the lack of explicit mycorrhizal associations, suggest areas for future model refinement (Wu et al., 2025; Braghiere et al., 2021; Zhu et al., 2019).

3.1.3 P cycle

The CNP model showed good overall performance in simulating phosphorus pools across ecosystem compartments (Fig. 5c). For plant components, the model accurately reproduced organic P pools, with slight overestimations of 5.0 %, 2.8 %, and 10.0 % for leaf, wood, and root compartments, respectively. For the soil P, the CNP simulated a lower value (1.58 g P m⁻²) than observed, but within its range (1.8 ± 0.6). Those organic P pools have the same pattern as organic N pools for CNP simulations, as C-N-P is coupled through stoichiometry.

The simulated inorganic P content (0.8 g P m⁻²) fell within the observed range (0.48–1.6 g P m⁻²). Additionally, the simulated net P mineralization rate (0.54 g P m⁻² yr⁻¹) was comparable to observations from tropical forests (0.67 ± 0.14 g P m⁻² yr⁻¹; Jiang et al., 2024b). The model successfully reproduced the observed levels of various P pools overall; however, it significantly underestimated plant P litterfall rates by 69 % (Table 5). This discrepancy suggests potential limitations in the model's representation of nutrient-related processes, such as plant nutrient resorption mechanisms. Nutrient resorption is a crucial physiological process through which plants adapt to varying N and P availability in ecosystems. In our model, we implemented a fixed resorption coefficient (Table 4), which may oversimplify the dynamic nature of nutrient resorption. Additionally, our model does not account for the reciprocal effects of nitrogen and phosphorus availability on nutrient resorption dynamics, where N availability influences P resorption efficiency and vice versa (See et al., 2015; Li et al., 2019). This simplified representation likely contributes to the contrasting patterns observed in plant nutrient litterfall rates, which overestimate N litterfall while underestimating P litterfall. Plants typically adjust their nutrient resorption efficiency in response to both internal nutrient status and external resource availability (Mao et al., 2015; Reed et al., 2012; Aerts and Chapin, 2000; Aerts, 1996). The fixed resorption coefficients in the current model structure may not capture these adaptive responses, potentially leading to unrealistic nutrient cycling patterns, especially under varying environmental conditions.

The CNP configuration successfully captured the steady-state P distributions across ecosystem pools despite some discrepancies in P cycling processes. Further refinements in P cycling processes, particularly in plant-soil P transfer mechanisms and plant internal P recycling, would be valuable for improving model performance (Jiang et al., 2019, 2024a). However, these improvements are currently constrained by limited observational data, as data scarcity remains a significant challenge for C-nutrient coupled modeling (Achat et al., 2016; Reed et al., 2015). Future research should prioritize comprehensive field measurements of P cycling processes, including plant P resorption efficiency, soil P transformation rates, and plant-soil P transfer dynamics. Such empirical data would not only help validate and improve model performance but also enhance our understanding of terrestrial P cycling and its interactions with C and N cycles in terrestrial ecosystems.

3.2 Evaluate the model-data fusion module

To evaluate the efficiency of the integrated data assimilation module, we compared the carbon fluxes from CNP simulations with default and optimized parameters (Figs. 6 and 7). The optimization showed varied improvements across different carbon flux components. For gross primary productivity (GPP), both default and optimized simulations captured the seasonal patterns well, with only a minor improvement in RMSE from 10.94 to 10.69 and a slightly increased correlation coefficient from 0.53 to 0.57 after optimization (Fig. 6a, e).

Figure 6 Comparison of weekly observed and simulated carbon fluxes using default parameters and optimized parameters for the Tiantong site in 2021. (a–c) Time series of observed (black dots) and simulated values with default parameters (blue line) and optimized parameters (red line), where the optimized results are derived from 1000 parameter sets randomly selected from 10 000 accepted parameter sets during the data assimilation process (shaded areas represent standard deviation). (d–f) Scatter plots of simulated versus observed values corresponding to the time series above, where the dashed line represents the 1 : 1 line. CC, correlation coefficient; RMSE, root mean square error.

Figure 7 Diurnal patterns of hourly net ecosystem exchange (NEE) across different months simulated by the CNP model configuration before (default) and after data assimilation (MCMC) compared with observations. Black lines represent observational data with shaded areas indicating ±1 standard deviation (SD). Colored lines indicate model simulations with shaded areas showing their respective ±1 SD. Root mean square errors (RMSE) between model outputs and observations are colored in blue for simulations with default parameters and in red for simulations with accepted parameters.

The photosynthetic capacity per unit area and photosynthetic surface area, indicated by V_cmax and SLA respectively, are key determinants of GPP. Both V_cmax and SLA were adjusted within their reference ranges during data assimilation (Fig. 8). Although these parameters showed compensatory effects in their adjustments, their combined effect still demonstrated a tendency to enhance GPP (Fig. 6a, e). Notably, the systematic underestimation of GPP, particularly during the growing season, suggests the need for improving current carbon cycle process representations. These improvements should include (1) the soil moisture control on stomatal conductance specific to evergreen broadleaf forests (Weng and Luo, 2008) and (2) the calculation of sunlit and shaded leaf proportions through more accurate clumping index parameterization in the two-leaf model (Wang et al., 2024a; Bi et al., 2022; Yan et al., 2017).

Ecosystem respiration (ER) showed more substantial improvement with data assimilation, with RMSE decreasing from 11.03 to 6.72 g C m⁻² per week, particularly in reducing the high-frequency fluctuations present in the default simulation (Fig. 6b). This improvement in ER led to a notable improvement in NEE, where the RMSE decreased from 14.21 to 8.83 g C m⁻² per week, and the correlation coefficient improved dramatically from −0.03 to 0.51. The significantly improved representation of carbon exchange dynamics with parameter optimization is further confirmed by the diurnal patterns across months (Fig. 7), with reduced RMSE in most months (7 out of 12). However, certain limitations persist, notably the underestimated NEE during midday hours in the growing season, primarily attributed to underestimated GPP, which requires further investigation.

The enhancement in ER and NEE primarily resulted from the efficiently constrained key parameters (Table 6, Fig. 8) based on the validated state variables (Fig. 5). While the default parameters achieved reasonable state variables, the response of state variables to new meteorological forcing conditions required adjustment (Ma et al., 2021). For instance, the Q₁₀ and soil carbon residence time (T₆–T₈) are well-constrained in our case. The temperature sensitivity parameter represents microbial responses to soil temperature, and carbon residence times serve as a proxy for microbial accessibility to carbon substrates, rather than just soil carbon properties, both of which are related to heterotrophic respiration. Through the optimization of these parameters, the CNP model effectively reduced the high-frequency fluctuations present in the default simulation and better captured the observed temporal dynamics.

Table 6 Target parameters, their ranges, mean values and maximum likelihood estimation (MLE) of the posterior distribution. Q₁₀ represents temperature sensitivity; SLA, specific leaf area; and V_cmax, maximum carboxylation rate. T₁–T₉ indicate turnover times for individual pools.

Figure 8 Posterior distributions of model parameters derived from Bayesian calibration. Grey shaded areas represent parameter posterior distributions, with red and blue vertical lines indicating posterior means and default values, respectively. The parameters (listed in Table 6) include Q₁₀, SLA, V_cmax, and carbon residence time parameters (T₁–T₉). The corresponding numerical values are shown in matching colors.

Data assimilation substantially improved CNP model performance in carbon flux simulation, highlighting the potential for applying our developed model to other flux sites without tedious manual calibration procedures. Given that parameter optimization can potentially compensate for structural deficiencies in models (e.g., the equifinality issue; Luo et al., 2016, 2009; Sierra et al., 2015), it's understandable that models with different nutrient coupling schemes can generate similar performance with optimized parameters (Fig. S1, Text S1 in the Supplement). However, while parameter optimization can help the C-only model fit historical data, it may result in unrealistic parameter values (Fig. S2) and essentially “bakes in” current nutrient conditions without representing the underlying processes, thereby compromising its predictive capacity for future scenarios.

4 Conclusions

In this study, we developed and evaluated a process-based CNP-coupled model for subtropical evergreen broadleaf forest. The CNP configuration demonstrated superior performance compared to C-only and CN models across most biogeochemical pools and fluxes, effectively addressing the overestimation issues prevalent in models with simplified biogeochemical processes. The incorporation of phosphorus cycling mechanisms proved crucial for capturing ecosystem dynamics in these phosphorus-limited systems, providing an essential foundation for predicting subtropical evergreen broadleaf forest responses to climate change. Beyond mechanistic improvements, site-scale models like TECO-CNP can fully leverage rich, localized datasets, including forest inventory records, experimental manipulations, and eddy covariance measurements, to constrain model parameters and processes. This integration is crucial because unobserved or weakly observed processes cannot be reliably constrained through data assimilation alone (Luo et al., 2011). TECO-CNP is designed to facilitate the fusion of such multi-process information, thereby enabling more mechanistic and robust representations of ecosystem C-N-P dynamics. Furthermore, we implemented and evaluated a model-data fusion framework using the MCMC algorithm, which significantly improved the simulation of carbon fluxes. The optimization of key parameters, including those that control photosynthetic capacity, temperature sensitivity, and carbon turnover rate, effectively reduced simulation uncertainties and enhanced model performance. The success of the data assimilation approach not only demonstrates its effectiveness in current model optimization but also provides a promising path for future model improvement and applications across diverse ecosystems. More importantly, integrating data assimilation frameworks with site-level biogeochemical models facilitates a synergistic loop between experimental findings and model development, enhancing our understanding of the nutrient cycle processes and our ability to make reliable predictions. This integrated approach provides a robust framework for improving ecosystem models and advancing our understanding of nutrient cycling in response to environmental changes.