Coupling of Phosphorus Processes With Carbon and Nitrogen Cycles in the Dynamic Land Ecosystem Model: Model Structure, Parameterization, and Evaluation in Tropical Forests

The biogeochemical processes of carbon (C), nitrogen (N), and phosphorous (P) are fully coupled in the Earth system, which shape the structure, functioning, and dynamics of terrestrial ecosystems. However, the representation of P cycle in terrestrial biosphere models (TBMs) is still in an early stage. Here we incorporated P processes and C‐N‐P interactions into the C‐N coupled Dynamic Land Ecosystem Model (DLEM‐CNP), which had a major feature of the ability in simulating the N and P colimitation on vegetation C assimilation. DLEM‐CNP was intensively calibrated and validated against daily or annual observations from four eddy covariance flux sites, two Hawaiian sites along a chronosequence of soils, and other 13 tropical forest sites. The results indicate that DLEM‐CNP significantly improved simulations of forest gross and net primary production (R2: 0.36–0.97, RMSE:1.1–1.49 g C m−2 year−1 for daily GPP at eddy covariance flux sites; R2 = 0.92, RMSE = 176.7 g C m−2 year−1 for annual NPP across 13 tropical forest sites). The simulations were also consistent with field observations in terms of biomass, leaf N:P ratio and plant response to fertilizer addition. A sensitivity analysis suggests that simulated results are reasonably robust against uncertainties in model parameter estimates and the model was very sensitive to parameters of P uptake. These results suggest that incorporating P processes and N‐P interaction into terrestrial biosphere models is of critical importance for accurately estimating C dynamics in tropical forests, particularly those P‐limited ones.


Introduction
In many terrestrial ecosystems, the availability of soil nutrients shapes their structures and functions, including vegetation productivity, biodiversity, succession, and interactions of plant, animal, and microbial communities (Vitousek, 2004). As one of the most important nutrients involving with cellular energy cycles and being responsible for building the molecules DNA and RNA, phosphorus (P) is fundamental to all living organisms (Ruttenberg, 2003;Vitousek et al., 2010). In terrestrial ecosystems, P comes primarily from parent material weathering and atmospheric deposition and is lost through nutrient leaching, soil erosion, and fire (Filippelli, 2008;Newman, 1995;Wang et al., 2015). Ecosystems with long-time development can reach a final steady state with P limitation, in which the ecosystem becomes depleted in parent P and a large fraction of the P should be bound up in occluded fraction (Buendía et al., 2010;Porder et al., 2007;Walker & Syers, 1976). P limitation on primary productivity at the level of individual species is recognized to be widespread in tropical forests because of the absence of glaciation and the relatively rapid weathering of parent material P as a result of warmer and wetter climate (Davidson et al., 2004;Tanner et al., 1998;Turner et al., 2018;Vitousek et al., 2010). Previous studies of soil P fractions and fertilization addition experiments also demonstrate that P limitation on plant productivity widely exists in tropical areas (Vitousek, 2004;. P limitation on primary productivity, therefore, plays a key role in the global carbon (C) cycle (Lloyd et al., 2001;Wieder et al., 2015). Moreover, P limitation on primary productivity may intensify in the future as a result of rising atmospheric CO 2 and nitrogen (N) deposition, which can weaken the carbon sequestration capability of terrestrial ecosystems (Goll et al., 2012;Zhang et al., 2014).
The interactions of C, N, and P play critical roles in regulating carbon uptake and storage and N, P biogeochemical processes in terrestrial ecosystems (Ågren et al., 2012;Guignard et al., 2017;Wang et al., 2017). P and N are both involved in critical C processes such as photosynthesis and decomposition. Leaf N and P can determine leaf photosynthetic capacity by regulating the maximum rate of carboxylation (V cmax ) and the electron transport capacity (J max ) (Domingues et al., 2010). A large amount of N is needed during the production of photosynthesis enzymes Rubisco and light-capturing proteins of thylakoid (Evans, 1989). N limitation could limit carboxylation capacity and electron transport rates (Hikosaka, 2004;Wang et al., 2018). P is required for many transformations of P rich molecules (ATP, NADP, and sugar phosphates from the Calvin cycle) and the regeneration of ribulose-1,5-bisphosphate (RUBP) (Farquhar et al., 1980). P limitation could reduce light-use efficiency, electron transport rates, regeneration of ribulose bisphosphate, and the allocation of leaf N to photosynthetic processes (Wang et al., 2018). Field experiments demonstrate that N and P fertilizer addition can increase forest trunk growth and biomass increment (Harrington et al., 2001;Newbery et al., 1999;Tanner et al., 1990;Vitousek et al., 1993). Additionally, multiple lines of evidence show that N and P interactions play important roles in affecting available N and P in soil (Marklein & Houlton, 2012;Reed et al., 2007;Treseder & Vitousek, 2001). The supply of P can influence N fixation rates, such that more supply of P can increase the N inputs and N availability in terrestrial ecosystems (Reed et al., 2007). Additionally, extra N is required by plants and microbes to produce more N-rich extracellular phosphatase enzymes that cleave ester-P bonds in soil organic matter increasing availability of P (Mcgill & Cole, 1981;Olander & Vitousek, 2000;Wang et al., 2007).
Since the industrial revolution in the 1850s, the global P cycle has been substantially disturbed and altered by anthropogenic activities (Cordell et al., 2009). Human has unprecedentedly accelerated the P cycle by fertilizer production and application in the agricultural lands to increase food production (Elser & Bennett, 2011). The P scarcity was exacerbated in the context of increasing atmospheric CO 2 concentration and N deposition indirectly induced by human activities (Jiang et al., 2019). Thus, there is an important need to represent P processes in the terrestrial biosphere models (TBMs) for better understanding and simulating C, N, and P interactions. This endeavor can help to improve the predictions of the C cycle in terrestrial ecosystems under a changing global environment. Meanwhile, an in-depth understanding of the P cycle will be important for decision makers to develop feasible policies for managing P nutrient sustainability in the context of the enhanced anthropogenic disturbances.
The incorporation of P cycle is an important task in the development of terrestrial biosphere models, particularly in tropical forests where P availability is often presumed to limit ecosystem primary production (Reed et al., 2015). Within the past decades, TBMs have evolved from the first-generation C only models (Houghton et al., 2001;Lieth, 1975) to C-N interactions models (Gerber et al., 2010;Thornton et al., 2007;Yang et al., 2009;Zaehle et al., 2010;Zaehle & Dalmonech, 2011). In recent years, the P cycle has been incorporated into several TBMs, such as Community Land Model (CLM-CNP) , CABLE-CNP (Wang et al., 2007), CASA-CNP (Wang et al., 2010), and Organizing Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) . There is mounting evidence that the representation of P in TBMs is important for our understanding and prediction of ecosystem dynamics, particularly for tropical forests with strong P limitation (Fleischer et al., 2019;Goll et al., 2012Goll et al., , 2018. Goll et al. (2012), Wang et al. (2010), and , and Yang,  highlighted the importance of incorporating the P cycle in TBMs, particularly in tropical ecosystems. Several studies also demonstrated that coupling P dynamics into TBMs could reduce the simulated terrestrial C sink in Amazon basin and other tropical regions due to increasing atmospheric CO 2 concentrations (Fleischer et al., 2019;Yang et al., 2016Yang et al., , 2019Zhang et al., 2011). However, it is still challenging for terrestrial ecosystem modelers to provide more realistic C-N-P models because some P processes and their interactions on C and N processes are still not well understood (Achat et al., 2016;Jiang et al., 2019;Reed et al., 2015).
In this study, we developed a process-based P module on the platform of the Dynamic Land Ecosystem Model (DLEM) by considering P impacts on vegetation and soil biogeochemical processes, which upgraded the coupled C-N model (DLEM-CN) into the coupled C-N-P model (DLEM-CNP). Although some existing models seem to work, the development of the DLEM-CNP model is still worthwhile because there may be more than one suitable model representation for one process (Beven, 2008). Here we first describe the structure of the P module in the DLEM-CNP. Then we evaluate the performance of the model to reproduce daily and seasonal C fluxes of four eddy covariance fluxes at tropical forest sites. After that, we evaluate model performance against field data collected from long-term fertilization experiments in a soil formation chronosequence in Hawaii (Harrington et al., 2001;Ostertag, 2001). Additionally, the simulated annual mean net primary production (NPP) is evaluated against NPP estimates developed by Clark et al. (2001) based on existing field measurements, which have been used as benchmarks for validating biogeochemical models. Finally, we test the robustness of the model by systematically varying key parameters related to P processes.

The DLEM
The DLEM is a highly integrated process-based Terrestrial Biosphere Model aiming to provide predictive understanding of terrestrial ecosystem patterns and processes at multiple spatial and temporal scales in the context of global change. DLEM simulates the consequences of natural and anthropogenic disturbances on the structure and functioning of terrestrial ecosystems and their feedbacks to human and natural systems on Earth's land surface and terrestrial-aquatic continuum. More specifically, the model has integrated biophysical, hydrological, and biogeochemical processes (C, N cycling), vegetation dynamics, disturbances including natural and anthropogenic stresses (e.g., changes in climate, atmospheric composition, land use land cover patterns, intensive management on crops and forests, and wildfire), and the synergy effects with the human system ( Figure 1). The DLEM model is capable of making daily and spatially explicit estimates on fluxes of water, greenhouse gases (including CO 2 , CH 4 , and N 2 O), dynamics of soil C, N, and water pools in terrestrial ecosystems, and the associated river discharge, riverine export of C and N from land to the ocean (Liu et al., 2013;Lu et al., 2012;Lu & Tian, 2013;Pan et al., 2020;Ren et al., 2007;Tian et al., 2011. In the DLEM, the simulation time step for soil thermal and hydrological processes (e.g., evapotranspiration, water interceptions, and water movement between grids) is 30 min, while physiological and biogeochemical processes (such as photosynthesis, plant respiration, organic matter decomposition, nitrogen mineralization and immobilization, etc.) are daily. Land use/land cover change and the associated carbon emissions are simulated annually. Spatial resolution for the DLEM is flexible ranging from hundreds of meters to tens of kilometers depending on the area of the study domain.
DLEM is driven by various input data sets mainly in four major types: (1) time invariant data (such as topography, soil properties, DEM, river network, and plant biophysical characteristics), (2) climate (including average temperature, minimum temperature, maximum temperature, precipitation, and shortwave radiation), (3) atmospheric chemistry (e.g., atmospheric CO 2 and nitrogen deposition), and (4) land use history (including land cover conversion and land management). Currently, part of the time-series data sets, such as climate is required as daily time step, while the other driving data sets, for example, atmospheric CO 2 concentration, nitrogen deposition, land use history, and land management (irrigation and N fertilizer application) are applied with monthly or annual time step.
As a spatially-explicit terrestrial ecosystem model, the basic simulation unit of DLEM is a single grid with corresponding coverage area. In the DLEM, we use a cohort structure to represent multiple plant functional types (PFTs) and land covers in each grid. Totally seven types of land cover and water bodies can be specified in each grid, including vegetation cover, impervious surface, lake, stream, sea, bare ground, and glacier. DLEM builds on the concept of plant functional types to describe vegetation distributions. Many different PFTs who adapt to local climate can coexist in the same grid, competing for light, water, and nutrient resources. The vegetation area consists of up to five PFTs, among which four types are reserved for natural PFTs and one for crops. All land cover types in each grid share a common soil water column, while physiological and soil biogeochemical processes for each PFT is simulated independently. In the DLEM, C enters the ecosystems mainly through vegetation CO 2 uptake during photosynthesis and leaves the ecosystems via ecosystem respiration, various land disturbances and harvest, and lateral transport to a water body. N enters natural ecosystems through biological N fixation and atmospheric N deposition. N leaves ecosystem through different pathways, including NH 3 volatilization, emissions of N 2 O, NO, and N 2 during nitrification and denitrification, N leaching from root zone to groundwater, and lateral N transport with surface runoff. The DLEM simulates the linkage between C, N, and water cycles across the plantsoil-atmosphere continuum, in which C and N processes had been included; however, P processes have not been incorporated.

Phosphorus Processes in the DLEM-CNP
In this study, we incorporated the P processes into the DLEM-CN, and the cycles of C, N, and P are fully coupled in the major processes including photosynthesis, allocation, turnover, nutrient uptake, and decomposition in the newly developed DLEM-CNP ( Figure 2). In the DLEM-CNP, organic P transfers into inorganic forms through mineralization in the soil, while inorganic P converts to organic form through immobilization and vegetation uptake.

Vegetation P Processes
In the DLEM-CNP, woody vegetation is composed of six pools (1: leaf, 2: sapwood, 3: fine root, 4: coarse root, 5: reproduction are living pools, and 6: heartwood is a dead pool), and herbaceous vegetation is composed of five pools (leaf, stem, fine root, coarse root, and reproduction). The stoichiometric relationship in each pool is determined by P uptake from soil, P allocation, tissue turnover, the minimum C:N ratio and C:P ratio (see Parameter Table in supporting information Table S1). Here we describe leaf nutrient limitations on photosynthesis, P uptake, and allocation processes.

10.1029/2020MS002123
Journal of Advances in Modeling Earth Systems Figure 2. Structure of (a) C-N-P cycles in DLEM-CNP:C enters the ecosystems through vegetation CO 2 uptake during photosynthesis. The plant biomass box consists of six C, N, P pools: (leaf, heartwood, sapwood, fine root, coarse root, and production). Litters are grouped into two added organic matter pools (AOM1 and AOM2) with different residence time. Soil organic matter has six pools: three microbial pools, namely, soil microbial 1(SMB1), soil microbial 2 (SMB2), and soil microbial residues (SMR), and three slow soil organic matter pools, namely, native organic matter (NOM), passive soil organic matter (PSOM), one dissolved organic matter (DOM). All these pools have corresponding N, P pools with specific C:N:P ratios. Soil available N box includes ammonium (NH 4 + ), nitrate (NO 3 − ). P in soil includes dissolve inorganic P pool, secondary mineral P pool, and occluded P pool. NP limitation has impacts on photosynthesis, carbon allocation, soil mineralization and immobilization processes. (b) P module: P enters ecosystems in the form of dissolved inorganic P from weathering of minerals, atmospheric deposition, and fertilizer. Dissolved inorganic P is the sole source for plants and microbes and can be reversibly adsorbed (secondary mineral P) onto soil particles or lost through leaching. Secondary mineral P can transform into occluded P, which is irreversibly lost to biota. When plants take up P from soils, it enters the plant allocating to growing plant tissues. When plant tissue is shed, part of the phosphorus is resorbed, while the rest enters the litter pools, from where it is either transformed into soil organic matter or mineralized.

P and N Colimitation on Photosynthesis
We adopted the N-P colimitation equation from the MIKE 3 ecohydrodynamic three-dimensional model developed by the Danish Hydraulics Institute, which is originally applied to modeling N and P limitations on phytoplankton growth (Lessin et al., 2007). This is a biochemical substitution approach, which presumes a shortage of one element can be alleviated by the substitution of another. Compared with most models incorporating both N and P simply describe nutrient limitation by using the Liebig's law of minimum (e.g., ORCHIDEE and CLM), this N and P interactive colimitation scheme has the capability of simulating the interactive effects of N and P limitation on photosynthesis.
In the DLEM-CNP, vegetation photosynthesis is simulated using the two-big-leaf model that separates plant canopy into sunlit and shaded leaves (Pan et al., 2014. More detailed equations about photosynthesis can be found in the supporting information. Accordingly, leaf P concentration is calculated separately for sunlit and shaded leaves to estimate P limitation on photosynthesis. To account for the interactive colimitation effect of N and P, the joint nutrient scalar f(NP)is calculated as follows: where f sun (N) and f shade (N) are sunlit and shaded leaves N limitation coefficients, and f sun (P) and f shade (P) are sunlit and shaded leaves P limitation coefficients. Leaf nutrient limitation coefficients are applied to impact V cmax , through which J max is also affected.
Specifically, the leaf P limitation coefficients f sun (P) and f shade (P) are estimated in a similar way as the leaf N limitation coefficients in the DLEM (Pan et al., 2014), where Pcon sun _ leaf and Pcon shade _ leaf are the P concentration (g P/g dry matter) for the sunlit and shaded leaves and Pcon leaf _ max is the maximum leaf P concentration (g P/g dry matter). The formulas of the three variables are as follows: where LAI sun and LAI shade are leaf area index for sunlit and shaded leaves, LAI is the total leaf area index of plant canopy, P leaf is leaf P content (g P m −2 ), C leaf is leaf carbon content (g C m −2 ), 0.45 is the conversion factor from carbon to dry biomass, k is the canopy extinction coefficient, and CP min,leaf is the minimum C:P ratio for leaves.

P Uptake
Plant P uptake rate (P up , g P m −2 day −1 ) is determined by the potential uptake rate (P pot , g P m −2 day −1 ), vegetation P demand (P demand , g P m −2 day −1 ), and available dissolved soil inorganic P (P dip , g P m −2 day −1 ):

Journal of Advances in Modeling Earth Systems
P up _ pot is the maximum vegetation P uptake rate (P up _ max , g P m −2 day −1 ) downregulated by volumetric water content (W soil , %), and P dip : where k up (g P m −2 day −1 ) is P dip concentration at which vegetation P uptake is half of its maximum rate and θ fc is the volumetric water content at field capacity (%). The demanded P for vegetation to keep the minimum C:P ratio, P demand (g P m −2 day −1 ), is estimated as where i denotes the vegetation living pools (1: leaf, 2: sapwood, 3: fine root, 4: coarse root, and 5: reproduction root); C i and P i are organic C and P pool of the vegetation (g C m −2 , g P m −2 ); CP min,i is a prescribed parameter representing the minimum C:P ratio in each plant pool.

P Allocation
The P uptaken by vegetation from soil is allocated to vegetation pools. The total allocable vegetation P (P alloc , g P m −2 ) is calculated as the sum of P in each pool after resorbing from shed leaf and P uptaken from soil. Part of vegetation P is transferred to litter organic P pools according to leaf litter C:P ratio (supporting information Figure S1) during vegetation turnover processes (such as leaf shedding, tree mortality, etc.). The nutrients resorption from leaf is simulated through the difference C:P ratio between leaf and leaf litter. For simplification purposes, we assumed that the P in plant is redistributed at daily allocation time step after P resorbing from shed tissue and being uptaken from soil.
P alloc is allocated to vegetation living pools (leaf, sapwood, coarseroot, fineroot, and reproduction) according to the stoichiometric relationship of C and P in each vegetation pool. P content in each vegetation pool after P allocation (P i * , g P m −2 ) is calculated as follows: where CP min,i is the minimum C:P ratio and i represents the ith vegetation pool.

Soil P Processes
In soil, we define three inorganic P pools, including dissolved inorganic P (P dip , g P m −2 ), secondary mineral P (P sm , g P m −2 ), and occluded P (P oc , g P m −2 ) ( Figure 2). We assume that soil P dip is the only readily available P for plant uptake: 10.1029/2020MS002123

Journal of Advances in Modeling Earth Systems
P dip is replenished by apatite rock weathering (P w , g P m −2 day −1 ), atmospheric deposition (P dep , g P m −2 day −1 ), fertilizer application (P fer , g P m −2 day −1 ), desorption of phosphorus from secondary mineral P pool (P des , g P m −2 day −1 ), and mineralization of organic phosphorus (P min , g P m −2 day −1 ). Meanwhile, P dip is consumed by plant root uptake (P up , g P m −2 day −1 ), microbial immobilization (P imb , g P m −2 day −1 ), adsorption by fine soil partials or cations (P ads , g P m −2 day −1 ), occlusion (P ocl , g P m −2 day −1 ), and lose by leaching and erosion (P lose_p , g P m −2 day −1 ).
For the soil organic matter (SOM), DLEM includes three microbial pools (SMB1, SMB2, and SMR), one dissolved organic matter pool (DOM), one native organic matter pool (NOM), and one passive organic matter pool (PSOM) (Figure 2a). Each soil organic pool has a prescribed stoichiometric relationship for C, N, and P (Table S1). Soil organic P accumulates input from litter P, transfer to soil inorganic P through mineralization processes, and leave soil column during leaching and erosion processes.

Weathering
Parent material weathering is the primary source of soil inorganic P, which is controlled by the quantity of apatite rock substrate and physical soil conditions (Guidry & Mackenzie, 2003). In the DLEM-CNP, the P weathering algorithm is adopted from the Century model (Parton et al., 1988), which considers the influences of phosphate rock availability, soil temperature and moisture, and soil texture: where P par is the parent P in a soil column (g P m −2 ), obtained from soil P synthesis data , f weather (T soil , W soil ) represents the effect of soil temperature and moisture on weathering, f weather (clay+silt) is the soil texture factor on weathering, clay and silt are the fractions of clay and silt in soil. The effect of soil temperature and moisture is estimated as follows: where Θ wilt is the volumetric water content at wilting point (%), Θ fc is the volumetric water content at field capacity, and Θ is the volumetric water content.

P Exchange Between Dissolved Inorganic P and Secondary Mineral P
The dissolved inorganic P can be adsorbed by fine soil particles and cations to form secondary mineral P Smeck, 1985). Meanwhile, adsorbed P can be released to the dissolve inorganic P pool when the balance of adsorption and desorption is disturbed Williams et al., 1980). DLEM-CNP simulates the bidirectional P fluxes (R des−ads , g P m −2 day −1 ) between soil dissolve inorganic P pool and secondary mineral P pool using the algorithm in Knisel (1993), The secondary mineral P pool is initialized as follows: where positive values of P des−ads represent adsorption (P ads ) and negative values represent desorption (P des ), f ad − de (T soil ) and f ad−de (W soil ) are soil temperature and soil moisture coefficients for inorganic P adsorption and desorption processes, and PAI is the soil chemical coefficient as a function of soil pH and clay fraction.

P Occlusion
In soil, iron (Fe), aluminum (Al), and calcium (Ca) compounds have a large binding capacity for P adsorption (Filippelli, 2008). The occluded P by these compounds is unavailable for plant uptake. In DLEM-CNP, the occlusion of P sources from P sm and the P occlusion rate (P ocl ) are calculated following Chaubey et al. (2007) and Knisel (1993): where ras is the ratio between the occluded P pool and the secondary mineral P pool and K as is a flow coefficient parameterized as follows:

P Mineralization and Immobilization
The competition between microbes and plants for dissolved inorganic P is handled analogously to the competition for soil mineral N in DLEM. Gross P immobilization, gross P mineralization, as well as plant P uptake are calculated daily. In each timestep, microbes are given a higher priority in accessing nutrients than plants.
Soil P can either be mineralized (organic P to inorganic P) or immobilized (inorganic P to organic P), depending on carbon decomposition rate (C dec ) and C:P ratios of the source and destination organic carbon pools (C : P source and C : P destination ): where P min − imm is the P mineralization or immobilization rate (g P m −2 ) with a positive value representing mineralization (P min ) and a negative value representing immobilization (P imb ), f CO2 is the fraction of the decomposed carbon emitted from the soil as CO 2 . The decomposition rate of organic pools is simulated following the classic first-order decay algorithm (Parton et al., 1987(Parton et al., , 1988: where C dec is the decomposed carbon (g C m −2 day −1 ), C pool is the size of carbon pool (g C m −2 ), and k dec is the decomposition rate (day −1 ) of each pool regulated by soil temperature, soil moisture, nutrient availability, and soil texture:

10.1029/2020MS002123
Journal of Advances in Modeling Earth Systems f(NP) is the nutrient scalar which is controlled by N limitation and P limitation: where K max is the maximum decay rate (day −1 ); f(PM) and f(PI) are P scalars in mobilization and immobilization, respectively; P imm is the potential P immobilization rate; avp is the available soil inorganic P (g P m −2 ), and avp opt is the optimum available soil P (g P m −2 ). The N related equations and detailed equations in soil decomposition can be found in the supporting information.

P Loss
The loss of P in dissolved forms may occur with surface runoff or leaching, and the loss of P in particulate forms is usually associated with the erosion of soil mineral or organic particles (McDowell et al., 2004). Dissolved inorganic P leaching (P lch¯dip , g P m −2 day −1 ) is calculated as follows: where f flow is the runoff and drainage coefficient for leaching, P dipc is the total inorganic P content in soil (g P g −1 soil), and lchb dip is the desorption coefficient for inorganic P. f flow is calculated as follows: where q srun is the surface runoff (mm); q drain is the drainage runoff (mm); θ is soil water content (mm); q srun , q drain , and θ are simulated by the DLEM soil hydrology procedures (Liu et al., 2013). Wt soil is soil mass in the top 50 cm (g m −2 ).
Similarly, dissolved organic P leaching (P lch¯dop , g P m −2 day −1 ) is simulated as follows: where P dop is the total amount of soil dissolved organic P (g P m −2 ); lchb dop is the desorption coefficient for dissolved organic P (DOP, g P g −1 soil); P dopc is the concentration of DOP concentration in the soil column (g P g −1 soil): Export of particular organic P (P lch¯pop , g P m −2 day −1 ) is assumed to occur with soil erosion (R erosion , g soil m −2 day −1 ), which is calculated using the

Journal of Advances in Modeling Earth Systems
where P opc is the total organic P content in soil column (g P g −1 soil) and P op (g P m −2 ) is the total organic P in the soil P pools at a depth of 1 m.

Hawaii Sites
The two forest sites (Thurston and Kokee) at Hawaiian Islands ( Figure 3) provide an excellent opportunity for nutrient cycle research as they have similar climate but different soil development stages (Vitousek, 2004). The soil at the Thurston site (latitude = 19.41°N, longitude = 155.24°W) is 300 years old (0.3 ky), and vegetation growth is limited mainly by N, while soil at the Kokee site (latitude = 22.14°N, longitude = 159.62°W) is 4 million years old (4,100 ky) and vegetation growth is limited mainly by P (Harrington et al., 2001;Ostertag, 2001). The annual mean temperature at the two sites is around 16°C, and the mean annual precipitation is about 2,500 mm. Both sites have similar parent material, and Metrosideros polymorpha is the major tree species (Crews et al., 1995). Field observations, including plant production, vegetation biomass, soil organic matters, and leaf N:P ratio, are available. In addition, long-term (6-11 years) factorial fertilization experiments were conducted at the two sites (at the Thurston site since 1985 for 12 years fertilization addition and at the Kokee site since 1991 for 6 years fertilization addition) (Harrington et al., 2001;Ostertag, 2001;Vitousek, 2004). The fertilization experiments used a factorial design that consisted of all four combinations of zero or 100 kg ha −1 year −1 N and zero or 100 kg ha −1 year −1 P. These fertilization experiments can help validate model performance in representing the N and P limitation interactive effects. Clark et al. (2001) synthesized the data in the literature on NPP in old-growth tropical forests to produce a consistent data set on NPP. They developed consistent, documented estimates of the upper and lower bounds around total NPP for 39 diverse tropical forest sites, to serve as benchmarks for calibrating and evaluating biogeochemical models. These estimates based on existing field measurements and judgment that belowground production is bounded by the range of 0.2-1.2 ANPP (aboveground NPP). As total NPP has frequently been estimated by assuming that BNPP equals ANPP (Esser et al., 1997)

Journal of Advances in Modeling Earth Systems
on longitude and latitude information. For the sites with the same longitude and latitude, we averaged the NPP value. We excluded sites without longitude and latitude information. We also took off the Hawaii sites from the Clark dataset as these sites have been included in section 2.3.1.2. A summary of site information can be found in Table S3 (site name, location, rainfall, and temperature).

Input Data 2.3.2.1. Input Data for FLUXNET Sites
At the four FLUXNET sites, climate variables (daily precipitation, daily shortwave radiation, daily maximum temperature, daily average temperature, and daily minimum temperature) and elevation were extracted from the site observation (https://fluxnet.fluxdata.org/). Atmospheric CO 2 concentration data was obtained from NOAA (https://www.esrl.noaa.gov/gmd/ccgg/trends/data.html). The data of elevation, slope, and aspect were extracted from Global 30 Arc-Second Elevation product (GTOPO30; https://lta.cr. usgs.gov/GTOPO30). Soil property data including soil texture, pH, and bulk density were extracted from Global Soil Data Task (http://www.daac.ornl.gov). Soil parent P data was from Global Gridded Soil Phosphorus Distribution Maps at a half-degree resolution .

Input Data for Hawaii Chronosequence Sites
At the two Hawaii sites, climate variables (daily precipitation, daily shortwave radiation, daily maximum temperature, daily average temperature, and daily minimum temperature) were extracted from the Daily Surface Weather Data on a 1 km Grid for North America, Version 3 ( Thornton et al., 2017). Atmospheric CO 2 concentration data was obtained from NOAA (https://www.esrl.noaa.gov/gmd/ccgg/trends/data. html). The data of elevation, slope, and aspect were extracted from USGS 10-meter DEM product (https:// gis.ess.washington.edu/data/raster/tenmeter/hawaii/). The other input data, such as soil bulk density, soil texture, and soil pH, were obtained from Goll et al. (2017) (Table S2 in the supporting information).

Clark Tropical Sites Forcing Data
For these tropical sites, climate data was obtained from the National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) product (https://globalweather.tamu.edu/), which has 0.3°resolution and has been widely used in watershed modeling research (Dile & Srinivasan, 2014;Fuka et al., 2014). Other forcing data are from the same data source as FLUXNET and Hawaii sites.

Simulation Setup
Model simulation at each site followed a three-step procedure: an equilibrium run, a spin-up run, and a transient run. The equilibrium run was carried out to obtain an equilibrium state, which is assumed to be reached when the net carbon exchange is less than 0.1 g C m −2 year −1 , the change in soil water pool is less than 0.1 mm year −1 , the change in total N content is less than 0.1 g N m −2 year −1 , and the change in total P content is less than 0.01 g P m −2 between two consecutive 20 year periods. The aim of the model spin-up is to smooth the transition between the equilibrium run and transient run. After the spin-up run, we conducted the transient run with all the driving forces changed year by year.
At the two Hawaii sites, we implemented three fertilizer addition experiments using the DLEM-CNP: (1) adding 10 g N m −2 year −1 N fertilizer only, (2) adding 10 g P m −2 year −1 P fertilizer only, and (3) adding 10 g N m −2 year −1 N and 10 g P m −2 year −1 P fertilizer together. The fertilizer was applied on every day with 1/365 of the annual rate.
For tropical forest sites, the simulations of both DLEM-CN and DLEM-CNP were driven by the same input data from 1980 to 2013. The simulated annual mean NPP was compared with field observations.
A one-at-a-time (OAT) sensitivity analysis was performed for the four FLUXNET sites by modifying Prelated model parameters (see Table S1 for definition of the parameters), including the maximum vegetation P uptake rate (P up_max ), the C:P ratio in leaflitter (CP leaflitter ), concentration of P dip at which vegetation P uptake is half of the maximum rate (k up ), the minimum C:P ratio in sapwood (CP min,sapwood ), the minimum C:P ratio in leaf (CP min,leaf ), and the desorption coefficient for inorganic P (lchb dip ). In one simulation, we only modified one parameter by ±10% with respect to the reference value and held other parameters constantly.

Statistical Analysis
Nutrient use efficiency is a measurement of how well plants use the available mineral nutrients.

Journal of Advances in Modeling Earth Systems
Response ratio (RR) represents the measured or modeled plant production in the fertilizer treatment divided by its value under unfertilized conditions . We use RR to evaluate the response of vegetation to nutrient addition. The NPP with fertilizer and without fertilizer used here as multiple year average during the experiment period.
The coefficient of determination (R 2 ), root-mean-square error (RMSE), and normalized root-mean-square error (NRMSE) are used to evaluate the model performance. NRMSE is of the same order of magnitude as the coefficient of variance of the measured variables. The simulation is considered excellent with NRMSE < 10%, good if 10-20%, acceptable if 20-30%, and poor if >30% (Jamieson et al., 1991).

Simulated and Observed GPP at FLUXNET Sites
The capacity of DLEM-CNP to reproduce observed daily and seasonal cycles of GPP is tested at the four FLUXNET sites. To evaluate the effects of the representation of P-related processes in DLEM-CNP, a reference simulation was performed using the DELM-CN model.  Figure 4 demonstrates this for the daily step GPP stimulation. Explicitly accounting for the effects of P limitation on photosynthesis does not change the pattern of either the simulated daily or seasonal cycle of carbon fluxes. The simulated leaf nutrient limitation factor at the four sites shows that with including P limitation, the leaf N-P limitation factor is stronger than the leaf N limitation factor in DLEM-CN, thus lower the simulated GPP. At these sites, we output daily average foliar N content, foliar P content, V cmax , P mineralization, P weathering rate, and plant P uptake rate, available PO 4 to show more information about P pools and fluxes (Table 1). The leaf N:P ratio is between 20.3 and 25.9, which is higher than the optimum N:P ratio, implying that P limitation exists at these four sites. The V cmax ranges from 39.1 to 43.3 μmol CO 2 m −2 s −1 . Walker (2014) compiled global data set of photosynthetic rates and leaf nutrient traits including estimates of V cmax , leaf nitrogen content, leaf phosphorus content data from both experimental and ambient field conditions. We filtered Walker's data set to get tropical forests V cmax , leaf N, and leaf P then compared with our daily simulated results at FLUXNET sites. Walk's data set indicates the V cmax of tropical forests ranges from 33.66 to 140.88 μmol/m 2 s −1 , the leaf N range from 1.1 to 2.81 g/m 2 , and the leaf P range from 0.06 to 0.18 g/m 2 . Our simulation results of daily V cmax and leaf nutrient concentration fell within their ranges.

Carbon Stocks and Fertilizer Addition Experiments at Hawaiian Sites
The NPP and carbon stocks simulated by the DLEM-CNP at the Thurston (0.3 ky) and the Kokee (4,100 ky) sites are consistent with field observations for both sites (Table 2 and Figure 5). The simulated NPP is 785.8 ± 14.8 g C m −2 year −1 at Thurston and 725.7 ± 40.1 g C m −2 year −1 at Kokee, and the observed NPP is 789 ± 63.0 g C m −2 year −1 and 757 ± 73.0 g C m −2 year −1 at corresponding sites. The woody biomass (8,535.2 ± 37.9 g C m −2 ) was reasonably reproduced by the model with a slightly higher (4.6%) than observation at the Thurston site. While the simulated woody biomass (8,875.9 ± 74.4 g C m −2 ) showed 9.4% lower than observation at the Kokee site ( Figure 5). The soil organic carbon simulated by the DLEM-CNP is 13,301.6 g C m −2 at Thurston and 18,075.1 g C m −2 at Kokee, which is 12.9% and 24.5% lower than field measurements, respectively.
Nutrient use efficiencies (NPP divided by plant nutrient uptake) is an implicit plant property that depends on the tissue stoichiometry . The DLEM-CNP simulated NUE is 184.8 ± 8.8 g C g −1 N at Thurston site and 124.5 ± 20.5 g C g −1 N at Kokee site, which is underestimated by 19.6% and 44.6%, respectively, implying lower carbon productivity per plant nitrogen simulated by our model. The N uptake (4.26 ± 0.22 g N m −2 year −1 at Thurston site and 5.94 ± 0.90 g N m −2 year −1 at Kokee site) is overestimated by the model at both sites. The simulated PUE was comparable to observations (2.07 ± 0.04 g C mg −1 P vs. 3.22 ± 0.23 g C mg −1 P; 3.47 ± 0.08 g C mg −1 P vs. 3.86 ± 0.53 g C mg −1 P). The simulations of DLEM-CNP captured the pattern of higher NUE at N limited site and higher PUE at P limited site, which is consistent with Vitousek (2004) reporting that forests could acquire and use nutrients more efficiently at nutrient scarcity sites. Besides, the simulated P uptake (0.38 ± 0.01 g P m −2 year −1 at Thurston and 0.21 ± 0.02 g P m −2 year −1 at Kokee) is consistent with observations across two sites (Table 2).

Journal of Advances in Modeling Earth Systems
The simulated leaf N:P ratio at Thurston and Kokee is 13.3 ± 0.17 g N g −1 P and 17.9 ± 0.44 g N g −1 P, respectively. The values are close to the field observations, which are 12.6 ± 1.6 g N g −1 P and 17.3 ± 2.7 g N g −1 P, indicating a good simulation regarding leaf nutrient condition. Tissue N:P ratios are widely used as an indicator of nutrient availability (Koerselman & Meuleman, 1996;McGroddy et al., 2004). Commonly, foliage N: P ratio of less than 14 indicates N limitation and above 16 indicates P limitation (Koerselman & Meuleman, 1996). The modeled leaf N:P ratio shows that N limits plant growth at the 0.3 ky site (leaf N:P ratio < 14), and P limits plant growth at the 4,100 ky site (leaf N:P ratio > 16).
For fertilizer addition experiments, at Thurston site (0.3 ky), which is N limited, adding N significantly increased plant production (RR = 1.34), meanwhile adding P also increased plant production with a much smaller amount than adding N (RR = 1.15) ( Figure 6). At the Kokee site (4,100 ky), which is P limited, adding P results in a greater increase in plant production (RR = 1.46) and adding N slightly enhanced the plant production (RR = 1.10). At both sites, the simultaneous N and P addition results in the maximum increase of plant production, with RR is 1.48 at the 0.3 ky site (Thurston) and 1.45 at the 4,100 ky site (Kokee).
Leaf N content and leaf P content can reflect the nutrient status of the plant and can have direct effects on leaf photosynthesis and vegetation productivity. At both sites, model simulations show leaf N content increase when the N fertilizer added and leaf P content increase when the P fertilizer added, which agrees with the observation. Overall, leaf N, P contents are realistically simulated by the DLEM-CNP. The variation of Leaf N content and leaf P content is roughly captured by the DLEM-CNP under all nutrient treatment conditions (add N, add P, and add N-P). However, the model cannot catch the amplitude of leaf P content Note. NPP is Net Primary Production; N uptake is plant nitrogen uptake rate; P uptake is plant phosphorus uptake rate; NUE is nitrogen use efficiency; PUE is phosphorus use efficiency; leaf N:P ratio with g/g.

10.1029/2020MS002123
Journal of Advances in Modeling Earth Systems increment when P and N-P fertilizers were added to P limited site (Figure 7b). At the P limited Kokee site, after P and N-P fertilizers were added, the observation of leaf P content was 0.260 ± 0.1000 (%, dry weight) and 0.190 ± 0.040 (%, dry weight), but model results are only 0.062 ± 0.008 and 0.058 ± 0.004 (%, dry weight).

Simulations at Tropical Sites by DLEM-CNP and DLEM-CN
We expanded the model validation to more tropical forest sites and examined model performance on the multiyear average NPP. We compiled NPP records at 13 tropical forest sites in Clark et al. (2001), ranging from 800 to 1,600 g C year −1 (Figure 8). The NPP simulated by the DLEM-CNP is consistent with the observational data sets (R 2 = 0.92). The introduction of the P module improved NPP simulation at the tropical forest sites, decreasing the overestimated NPP by the DLEM-CN model across all the 13 sites. The R 2 is 0.92 for DLEM-CNP and 0.87 for DLEM-CN. RMSE is 176.7 g C m −2 year −1 for DLEM-CNP and 364.1.5 g C m −2 year −1 for DLEM-CN. NRMSE is 38.5% for DLEM-CNP and 79.4% for DLEM-CN. All three statistical indicators (R 2 , RMSE, and NRMSE) demonstrate the DLEM-CNP has better performance in modeling tropical forest NPP by including P processes and considering C-N-P interactions. Figure 9 shows the sensitivity of simulated GPP, foliage P, available PO 4 , leaf N:P ratio to P-related parameters, and potential uncertainties in model parameterization. Simulated GPP varies between −2.7% and 2.5% in relative to the simulated values in the default model parameterization (Figure 9a). GPP is most sensitive to the P uptake (P up _ max , k up ), reflecting the important role of plant P uptake on P limitation and GPP. Leaf P content (Figure 9b) is predominantly controlled by the P uptake (P up _ max , k up ) and P allocation (CP min,sapwood , CP min,leaf ) processes. P up _ max lead to −5.4% to 5.1% changes in Leaf P content. And leaf P content varies with change of −5.3-6.2% and −5.0-5.9% for parameter k up and CP min,leaf , respectively. Available PO 4 shows slightly response to the different parameter set. Notably, lchb dip is among the most

10.1029/2020MS002123
Journal of Advances in Modeling Earth Systems sensitivity for available PO 4 leading to a relatively small change −1.1-0.7%. Leaf N:P ratio shows a similar pattern with leaf P content of sensitivity to parameters. The sensitivity analysis demonstrates that changes in parameters result in small changes in model outcomes, and model is most sensitive to parameters on P uptake.

Model Evaluation at the Site Scale
The results presented here demonstrate the incorporation of key P dynamics within the DLEM. The model simulations are consistent with the observed ranges of GPP, NPP, P tissue concentrations, NUE, net N mineralization, and the response to the fertilizer addition experiments at the FLUXNET sites, Hawaiian sites, and other tropical sites.
Some TBMs with coupled C-N-P cycles have been reported previously (Goll et al., 2012Wang et al., 2010;. The DLEM-CNP model distinguishes itself from other models in the explicit consideration of the N-P colimitation (Equations 1 and 2) on photosynthesis. With this N-P colimitation scheme, DLEM-CNP has better potential for simulating N-P interaction effects on C assimilation in contrast to models without considering P effect on photosynthesis, only considering leaf P to calculate GPP or describing nutrients limitation on plant growth by using the Liebig's law of minimum (e.g., CLM-CNP and ORCHIDEE). In the Hawaiian site fertilizer addition experiments, the addition of P at the 0.3 ky N-limited Thurston site alleviated N limitation on vegetation growth. Also, the N addition at P-limited site slightly increased NPP. These responses demonstrated the innovative improvements of introducing the N and P interactive colimitation effect. Both theory and field empirical relationship support this N and P interactive colimitation needs to be considered (Domingues et al., 2010;Jiang et al., 2019;Walker et al., 2017;Wang et al., 2018). Considering N, P colimitation may increase the capability of models to predict future conditions in P-limited tropical forests, especially when combined with other environmental drivers (Domingues et al., 2010;Jiang et al., 2019). However, the mechanisms are not fully clear. Therefore, more data and experiments on foliar N, P concentrations and photosynthetic parameters are needed to help modelers develop more robust relationships.
DLEM-CNP diagnoses the leaf N, P concentration to support optimal leaf C:P, C:N ratios and GPP is regulated through the N-P limitation on V cmax when there is no sufficient available P and N to maintain the optimal leaf C:P, C:N ratios. Our simulation results of daily V cmax are at the lower end of the observation range. However, there are a range of V cmax values from different studies indicating our simulated V cmax is reasonable. Kattge et al. (2009) reported the observed V cmax for tropical broadleaf evergreen trees was 41 μmol CO 2 m −2 s −1 (Bonan et al., 2012). Some modeling studies for Amazonia also reported various V cmax values: 43 μmol CO 2 m −2 s −1 (Carswell et al., 2000); 58 μmol CO 2 m −2 s −1 (Mercado et al., 2006); 64 μmol CO 2 m −2 s −1 , based on the Lloyd et al. (2010) analysis of Domingues et al. (2005); and 68 μmol CO 2 m −2 s −1 (Lloyd et al., 1995). Mercado et al. (2009)   10.1029/2020MS002123

Journal of Advances in Modeling Earth Systems
At Hawaiian sites, the underestimated NUE could be due to the overestimation of N uptake, which is related to C:N ratio of biomass and biological N fixation rate. The C:N ratio of biomass is a very species dependent parameter, so it may tend to introduce bias at site level. Moreover, the biological N fixation rate is described as a plant-type-specific parameter in DLEM, which may not reflect the reality at Hawaiian sites.
Our model underestimated the absolute value of RR in nutrient fertilization experiments. Also, the model did not catch the amplitude of leaf P content increase when P or N-P fertilizer was added to P limited site (Figure 7b). Leaf N and P concentrations determine the leaf N-P limitation on GPP, so the smaller increase in leaf P content than observation resulted in the underestimated RR. Two reasons explain the underestimated leaf P content increase. DLEM-CNP simulated fertilizer application by assuming that annual fertilizer input was applied evenly each day. Specifically, in our simulation, the N and P fertilizers were added at a rate

10.1029/2020MS002123
Journal of Advances in Modeling Earth Systems of 0.027 (10/365) gN/gP m −2 per day. These fertilizers directly entered soil available N pool and soil dissolve inorganic P pool. The 0.027 gN/gP m −2 per day fertilizers did not lead to a substantial increase in soil available N and dissolve inorganic P pool, which impeded the nutrient uptake by plant. The other reason is that Harrington et al. (2001) observed P limitation results in disproportionally significant increases in P uptake after fertilization during the fertilization experiment. However, this mechanism strongly depends on species; it has not been included in our model. In the future, we will improve our model to simulate fertilizer with monthly or daily data to more accurately capture the fertilization effect.
The mismatch between the modeled SOM and observed values at Hawaiian sites can be attributed to the uncertainties in modeling decomposition process. In contrast to first-order decomposition models (Parton et al., 1987), SOC decomposition rates should depend not only on the SOC stock size but also on the size and composition of the decomposer microbial pool (Schimel & Weintraub, 2003) as well as carbon nutrients interactions (Six et al., 2002). Therefore, the SOC stocks are still far from certain in TBMs. Some research demonstrated P could have impacts on the microbial decomposition rate, for example, Qualls and Richardson (2000) suggested P enrichment influences the litter decomposition rate in the Everglades. Tian, Lu, et al. (2015) also emphasized that nutrient limitation on microbial activities are important structural uncertainties but being ignored or poorly represented by most models leading to uncertainties in the modeled SOC. All these studies implying a more realistic decomposition process are needed.

Uncertainty and Future Research
Although DLEM-CNP could well reproduce the observed forest productivity and biomass, several issues need to be addressed in the future. First, our model does not consider biochemical mineralization, a process that plants and microbes can produce phosphatase enzymes to mineralize P in soil organic matter but without mineralizing C and N (Mcgill & Cole, 1981). Wang et al. (2007) showed that neglecting biochemical P mineralization tended to overestimate the fraction of soil organic P and underestimated the fractions of the labile P. Neglect of biochemical mineralization may underestimate available P, causing stronger P limitation on plant productivity. Second, DLEM-CNP did not consider P effects on some N processes, such as the biological N fixation process. Several studies indicate that P availability may have impacts on N fixation (Edwards et al., 2006;Reed et al., 2007). Edwards et al. (2006) examined biological N fixation of swards response to elevated CO 2 under both high and low P availability, showing high P increased biological N fixation. Reed et al. (2007) suggested that P availability, possibly via regulation of N fixation, strongly influence N availability in recovering prairie soils. However, the mechanisms of P influence on N fixation have not yet been fully understood. Third, plants' adaptation to P limitation can change vegetation physiological characteristics (Vance et al., 2003). Plants evolve strategies for P acquisition and use in P-limiting environments, including decreased growth rate, increased growth per unit of P uptake, remobilization of internal P, modifications in carbon metabolism that bypass P-requiring steps, and alternative respiratory pathways (Vance et al., 2003). All these adaption strategies are supposed to lead to vegetation parameters change (e.g., V cmax , minimum leaf C:N ratio and minimum leaf C:P ratio).
In tropical areas, forests in different places have different P limitation levels and highly variable NPP (Castanho et al., 2013;Fyllas et al., 2009;Tanner et al., 1998). Our results at various tropical forest sites demonstrated that using a uniform value for key ecological parameters of tropical forests would decrease the simulation accuracy. In future work, a more detailed forest classification scheme could be developed to represent the spatial heterogeneity of tropical forests.
As P limitation on productivity is recognized to be widespread in tropical forests, at this stage, we only focus our model development in tropical forest areas. In the future, we will expand the model application to other terrestrial ecosystems.

Conclusions
The P cycle has been successfully coupled with C and N cycles in the DLEM. The DLEM-CNP fully incorporates C-N-P cycles in all pools (plant, litter, and soil organic/inorganic pools) and key biogeochemical processes. Evaluation of model performance at FLUXNET sites, Hawaiian long-term fertilization experiment sites, and comparison of NPP simulated by the DLEM-CNP and DLEM-CN with benchmarks in tropical forests demonstrate the improvement of the model's performance. Furthermore, we simulated the response of plants to fertilizer addition for evaluating the mechanisms of nutrient limitation on plants. Results imply a significant P impact on ecosystem C dynamics and highlight the innovative improvements of introducing the N and P interactive colimitation effect. Our results reveal interactions between C, N, P processes, indicating that the inclusion of the P cycle in the current TBMs is essential to better understand the impacts of global change on terrestrial ecosystems. With adequate parameterization, the DLEM-CNP model can be applied to simulate and predict the productivity of terrestrial ecosystems and C, N, P dynamics across the global land surface.

Data Availability Statement
The relevant data sets of this study are archived in the box site of International Center for Climate and Global Change Research at Auburn University (https://auburn.box.com/v/dlemcnpfolder).