Representation of the phosphorus cycle in the Joint UK Land Environment Simulator (vn5.5_JULES-CNP)

Most land surface models (LSMs), i.e. the land components of Earth system models (ESMs), include representation of nitrogen (N) limitation on ecosystem productivity. However, only a few of these models have incorporated phosphorus (P) cycling. In tropical ecosystems, this is likely to be important as N tends to be abundant, whereas the availability of rock-derived elements, such as P, can be very low. Thus, without a representation of P cycling, tropical forest response in areas such as Amazonia to rising atmospheric CO₂ conditions remain highly uncertain. In this study, we introduced P dynamics and its interactions with the N and carbon (C) cycles into the Joint UK Land Environment Simulator (JULES). The new model (JULES-CNP) includes the representation of P stocks in vegetation and soil pools, as well as key processes controlling fluxes between these pools. We develop and evaluate JULES-CNP using in situ data collected at a low-fertility site in the central Amazon, with a soil P content representative of 60 % of soils across the Amazon basin, to parameterize, calibrate, and evaluate JULES-CNP. Novel soil and plant P pool observations are used for parameterization and calibration, and the model is evaluated against C fluxes and stocks and those soil P pools not used for parameterization or calibration. We then evaluate the model at additional P-limited test sites across the Amazon and in Panama and Hawaii, showing a significant improvement over the C- and CN-only versions of the model. The model is then applied under elevated CO₂ (600 ppm) at our study site in the central Amazon to quantify the impact of P limitation on CO₂ fertilization. We compare our results against the current state-of-the-art CNP models using the same methodology that was used in the AmazonFACE model intercomparison study. The model is able to reproduce the observed plant and soil P pools and fluxes used for evaluation under ambient CO₂. We estimate P to limit net primary productivity (NPP) by 24 % under current CO₂ and by 46 % under elevated CO₂. Under elevated CO₂, biomass in simulations accounting for CNP increase by 10 % relative to contemporary CO₂ conditions, although it is 5 % lower compared to CN- and C-only simulations. Our results highlight the potential for high P limitation and therefore lower CO₂ fertilization capacity in the Amazon rainforest with low-fertility soils.

1 Introduction

Land ecosystems currently take up about 30 % of anthropogenic CO₂ emissions (Friedlingstein et al., 2020), thus buffering the anthropogenic increase in atmospheric CO₂. Tropical forests play a major role in the land C cycle, account for about half of global net primary production (NPP) (Schimel et al., 2015), and store the highest above-ground carbon among all biomes (Pan et al., 2011; Mitchard, 2018).

The C sink capacity of tropical forests may be constrained by nutrient availability for plant photosynthesis and growth (Vitousek and Howarth, 1991; Elser et al., 2007; LeBauer and Treseder, 2008) via P-related (Nordin et al., 2001; Shen et al., 2011) and/or N-related processes (DeLuca et al., 1992; Perakis and Hedin, 2002). Global process-based models of vegetation dynamics and function suggest a continued land C sink in the tropical forests, largely attributed to the CO₂ fertilization effect (Sitch et al., 2008; Schimel et al., 2015; Koch et al., 2021). However, many of these models typically do not consider P constraints on plant growth (Fleischer et al., 2019), which is likely to be an important limiting nutrient in tropical ecosystems, characterized by old and heavily weathered soils. The importance of nutrient cycling representation in Earth system models (ESMs) (and the lack thereof) was highlighted by Hungate et al. (2003) and Zaehle and Dalmonech (2011), showing the significance of nitrogen inclusion in ESMs for generating more realistic estimations of the future evolution of the terrestrial C sink. However, in the Coupled Climate C Cycle Model Inter-comparison Project (C4MIP), none of the participating ESMs included N dynamics (Friedlingstein et al., 2006). Following this, for the update of CMIP5 7 years later (Anav et al., 2013), 3 models out of the 18 used included N dynamics (Ji et al., 2014; Long et al., 2013; Bentsen et al., 2013). Although much progress has been made in the inclusion of an N cycle in ESMs so far, none of the CMIP5 models included P cycling, and in the most recent CMIP6 only one model includes P (ACCESSESM1.5 model) (Arora et al., 2020).

The long history of soil development in tropical regions, which involves the loss of rock-derived nutrients through weathering and leaching on geologic timescales (Vitousek et al., 1997, 2010), results in highly weathered soils. Soil P is hypothesized to be among the key limiting nutrients to plant growth in tropical forests (Hou et al., 2020; Vitousek et al., 1997, 2010), which is unlike temperate forest, where N is hypothesized to be the main constraint (Aerts and Chapin, 1999; Luo et al., 2004). Low P availability in tropical soils is related to the limited amount of unweathered parent material and organic compounds (Walker and Syers, 1976), active sorption (Sanchez, 1977), and high occlusion (Yang and Post, 2011), which all contribute to further reduce plant-available P. Although N limitation can impact the terrestrial C sink response to increasing atmospheric CO₂ by changing plant C fixation capacity (Luo et al., 2004), this can be partially ameliorated over time by input of N into the biosphere via the continuous input of N into ecosystems from atmospheric deposition and biological N fixation (Vitousek et al., 2010). P limitation is pervasive in natural ecosystems (Hou et al., 2020), and the lack of large P inputs into ecosystems, especially those growing on highly weathered soil, may make P limitation a stronger constraint on ecosystem response to elevated CO₂ (eCO₂) than N (Sardans et al., 2012; Gentile et al., 2012). This causes considerable uncertainty in predicting the future of the Amazon forest C sink (Yang et al., 2014).

There is evidence to suggest P limitation on plant productivity in the Amazon rainforest (Malhi, 2012), where it has been shown that the younger, more fertile western and southwestern Amazonian soils have higher tree turnover (Phillips et al., 2004; Stephenson and Van Mantgem, 2005) and stem growth rates (Malhi et al., 2004) and lower above-ground biomass (Malhi et al., 2006; Baker et al., 2004) compared to their central and eastern counterparts. Total soil P has been found to be the best predictor of stem growth (Quesada et al., 2010) and total NPP (Aragão et al., 2009) across this fertility gradient, and foliar P is positively related to plant photosynthetic capacity (V_cmax and J_cmax) in these forests (Mercado et al., 2011).

However, modelling studies are unable to reproduce observed spatial patterns of NPP and biomass in the Amazon, with one possible reason being the lack of inclusion of soil P constraints on plant productivity and function (Wang et al., 2010; Vicca et al., 2012a; Yang et al., 2014). Nevertheless, some modelling studies have focused on improving process and parameter representation using the observational data of spatial variation in woody biomass residence time (Johnson et al., 2016), soil texture, and soil P to parameterize the maximum carboxylation capacity (V_cmax) (Castanho et al., 2013). Results from these studies successfully represent observed patterns of Amazon forest biomass growth increases with increasing soil fertility. However, the full representation of these interactions and the impact of the soil nutrient availability on biomass productivity is still missing in most of ESMs.

So far, several dynamic global vegetation models have been developed to represent P cycling within the soil (Yang et al., 2013; Haverd et al., 2018) and between plant and soils for tropical forests in particular (Yang et al., 2014; Zhu et al., 2016; Goll et al., 2017). Furthermore, a comprehensive study included several models with C–N–P cycling and their feedbacks on the atmospheric C fixation and biomass growth in Amazon forests under ambient and elevated CO₂ conditions (eCO₂) (Fleischer et al., 2019). Despite these developments, data to underpin them and their projections, particularly for the tropics, is sparse and remains challenging, especially for the Amazon forest (Reed et al., 2015; Jiang et al., 2019). Moreover, due to the lack of detailed measurements, the P-related processes such as adsorption, desorption, and uptake represented in these models are under-constrained and likely oversimplified, and thus the future predictions of Amazon rainforest responses to eCO₂ and climate change are uncertain. To fill this gap, in this study we use data collected as part of the Amazon Fertilization Experiment (AFEX), the first project that focuses on experimental soil nutrient manipulation in the Amazon, with a comprehensive data collection programme covering plant ecophysiology, C stocks and fluxes, and soil processes (including P stocks). Thus, as opposed to previous P modelling studies, our model parameterization includes detailed P process representation using site measurements.

Here, we describe the development and implementation of the terrestrial P cycle in the Joint UK Land Environment Simulator (JULES) (Clark et al., 2011), the land component of the UK Earth System Model (UKESM), following the structure of the prior N cycle development (Wiltshire et al., 2021) and utilizing already tested and implemented state-of-the-art descriptions of P cycling of other land surface models (Wang et al., 2007; Zhu et al., 2016; Goll et al., 2017).

The model (JULES-CNP) is parameterized and calibrated using novel in situ P soil and plant data from a well-studied forest site in the central Amazon near Manaus, Brazil, with soil P content representative of 60 % of soils across the Amazon basin. The newly developed P component estimates the sorption of the soil organic and inorganic P based on the saturation status of the adsorbed P pools, which is unique compared to the other existing P models, and enables more realistic estimation of P adsorption and desorption processes. We first evaluate the model at our study site and at additional five test sites across the Amazon and in Panama and Hawaii. We then apply the model under ambient and eCO₂ following the protocol of Fleischer et al. (2019) to predict nutrient limitations on land biogeochemistry under these conditions. Predictions of the CO₂ fertilization effect in JULES-CNP are compared to those in current versions of the model with coupled C and N cycles (JULES-CN) and with the C cycle only (JULES-C).

2 Material and methods

2.1 JULES

JULES is a process-based model that integrates water, energy, C cycling (JULES-C) (Clark et al., 2011), and N cycling (JULES-CN) (Wiltshire et al., 2021) between the atmosphere, vegetation, and soil (Best et al., 2011; Clark et al., 2011). Vegetation dynamics are represented in JULES using the TRIFFID model, with nine distinct plant functional types (PFTs) (tropical and temperate broadleaf evergreen trees, broadleaf deciduous trees, needle-leaf evergreen and deciduous trees, C₃ and C₄ grasses, and evergreen and deciduous shrubs), as well as height competition (Harper et al., 2016). Leaf-level photosynthesis (Collatz et al., 1991, 1992) is scaled to estimate canopy-level gross primary productivity (GPP) using a multilayer approach that accounts for vertical variation of radiation interception and partition of sunlit and shaded leaves and associated vertical variation of leaf N and P exponential decrease through the canopy (Clark et al., 2011; Mercado et al., 2007, 2009), while the C:P and N:P ratios remain the same. NPP is estimated as the difference between GPP and autotrophic respiration for each living tissue (leaf, wood, root). NPP is then allocated to increase tissue C stocks and to spread, i.e. expand the fractional coverage of the PFT. The resultant PFT fractional coverages also depend on competition across PFTs for resources, e.g. light. Tissue turnover and vegetation mortality add C into the litter pools. Representation of soil organic C (SOC) follows the Rothamsted carbon model (RothC) equations (Jenkinson et al., 1990; Jenkinson and Coleman, 2008) defining four C pools: decomposable plant material (DPM) and resistant plant material (RPM), which receive direct input from litterfall, and microbial biomass (BIO) and humified material (HUM), which receive a fraction of decomposed C from DPM and RPM that is not released to the atmosphere. The limitation of N on SOC is applied to the vegetation and soil components using a dynamic C:N ratio to modify the mineralization and immobilization processes as described in Wiltshire et al. (2021). Note that the soil component of JULES-CN can either be run as a single box model or vertically resolved over soil depth (JULES-CN layered), and in this paper we build upon the vertically resolved version described in Wiltshire et al. (2021).

2.2 JULES-CNP

JULES-CNP includes the representation of the P cycle in JULES version (vn5.5), and it is built on existing and well-tested representations of P cycling in other global land surface models (Yang et al., 2014; Sun et al., 2021; Wang et al., 2007; Goll et al., 2017). It includes P fluxes within the vegetation and soil components and the specification of P pools and processes related to P cycling within the soil column (Fig. 1). A parent material pool is introduced to consider the input of weathered P. The adsorbed, desorbed, and occluded fractions of P for both organic and inorganic P are also represented. However, except for parent material and occluded P pools, all other pools are estimated at each soil layer. The description of changes in pools and associated relative fluxes are explained in detail in the next sections. Although JULES-CN includes N leaching and deposition, P leaching and deposition are not included in the current version of JULES-CNP.

Figure 1 JULES-CNP model scheme including P pools (grey boxes) and fluxes (arrows).

2.2.1 P pools

JULES represents eight P pools comprising organic and inorganic P: plant P (P_p) and soil P pools (in each soil layer (n)), litter P ($P_{{\mathrm{O}}_{\mathrm{l}}}$), soil organic P ($P_{{\mathrm{O}}_{\mathrm{s}}}$), soil inorganic P (P_in), organic sorbed P (P_org-sorp), inorganic sorbed P (P_inorg-sorp), parent material (P_pm), and occluded (P_occ) P comprised of both organic and inorganic P (all pools are in units of kg P m⁻²; see Fig. 1 and Tables 1 and 2).

Plant P pool is composed of leaf (P_leaf), fine root (P_root), and stem together with coarse root (P_stem), which are related to their associated C pools (C_leaf,C_root, C_stem, in units of kg C m⁻²) and fixed C to P ratios (${\mathrm{C}:\mathrm{P}}_{\mathrm{leaf}}{\mathrm{C}:\mathrm{P}}_{\mathrm{root}}{\mathrm{C}:\mathrm{P}}_{\mathrm{stem}}$) as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }{P}_{\mathrm{leaf}}={\displaystyle \frac{C_{\mathrm{leaf}}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{leaf}}}},\text{(2)}& {\displaystyle }{P}_{\mathrm{root}}={\displaystyle \frac{C_{\mathrm{root}}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{root}}}},\text{(3)}& {\displaystyle }{P}_{\mathrm{stem}}={\displaystyle \frac{C_{\mathrm{stem}}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{stem}}}}.\end{array}$$

Therefore, the plant P pool (P_p) is the sum of all vegetation P pools as follows:

$$\begin{array}{}\text{(4)}& P_{\mathrm{p}}=P_{\mathrm{leaf}}+P_{\mathrm{root}}+P_{\mathrm{stem}}.\end{array}$$

Description of the plant P pool (P_p) follows Zhu et al. (2016) and is estimated as the difference between the input, plant uptake $F_{\mathrm{P}}^{\mathrm{Up}}$ (Eq. 26) and output of this pool, and plant litter flux $F_{\mathrm{P}}^{\mathrm{lit}}$ (Eq. 28) with both fluxes expressed as follows (in units of kg P m⁻² yr⁻¹):

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{p}}}{\mathrm{d}t}}=F_{\mathrm{P}}^{\mathrm{up}}-F_{\mathrm{P}}^{\mathrm{lit}}.\end{array}$$

The litter P pool ($P_{{\mathrm{O}}_{\mathrm{l}}}$) is estimated as a sum of P_DPM and P_RPM pools over soil layers (n). Each pool is formed by the fluxes of plant litter input $(F_{\mathrm{P}}^{\mathrm{lit}}$) and the outgoing decomposed P (${\mathrm{dec}}_{\mathrm{P}}^{\mathrm{lit}}$), which expressed in Eqs. 28 and 29, respectively (in kg P m⁻² yr⁻¹). Furthermore, the plant litter input is modified based on the plant type material ratio α (in order to distribute the litter input based on the DPM $/$ RPM fraction) as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{DPM}}}{\mathrm{d}t}}& {\displaystyle }=F_{{\mathrm{P}}_n}^{\mathrm{lit}}\times \mathit{\alpha }-{\mathrm{dec}}_{{\mathrm{P}}_{\mathrm{DPM},n}},\text{(7)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{RPM}}}{\mathrm{d}t}}& {\displaystyle }=F_{{\mathrm{P}}_n}^{\mathrm{lit}}\times (\mathrm{1}-\mathit{\alpha })-{\mathrm{dec}}_{{\mathrm{P}}_{\mathrm{RPM},n}},\text{(8)}& {\displaystyle }{P}_{{\mathrm{O}}_{\mathrm{l}}}& {\displaystyle }=\sum _{n=\mathrm{1}}^N{P}_{{\mathrm{DPM}}_n}+\sum _{n=\mathrm{1}}^N{P}_{{\mathrm{RPM}}_n}.\end{array}$$

The soil organic pool ($P_{{\mathrm{O}}_{\mathrm{s}}}$) is represented as the sum of P_BIO and P_HUM. These pools are estimated from the difference between P inputs from total immobilized P $\left(F_{{\mathrm{immob}}_{\mathrm{P}}}\right)$ distributed between BIO and HUM based on a fixed fraction (0.46 for BIO, 0.54 for HUM) (Jenkinson et al., 1990; Jenkinson and Coleman, 2008), desorbed P $\left(F_{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}^{\mathrm{desorp}}\right)$, and P outputs from mineralized $(F_{{\mathrm{minl}}_{\mathrm{P}}}$), and adsorbed P fluxes $(F_{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}^{\mathrm{sorp}}$) (adsorption is shown in Eq. 40, and desorption is shown in Eq. 41) with the fluxes expressed (in kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{BIO}}}{\mathrm{d}t}}& {\displaystyle }=\mathrm{0.46}\times F_{{\text{immob}}_{{\mathrm{P}}_n}}+F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_{\mathrm{BIO},n}}^{\mathrm{desorp}}-F_{{\mathrm{minl}}_{{\mathrm{P}}_{\mathrm{BIO},n}}}\\ \text{(9)}& {\displaystyle }& {\displaystyle }-F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_{\mathrm{BIO},n}}^{\mathrm{sorp}},{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{HUM}}}{\mathrm{d}t}}& {\displaystyle }=\mathrm{0.54}\times F_{{\text{immob}}_{{\mathrm{P}}_n}}+F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_{\mathrm{HUM},n}}^{\mathrm{desorp}}-F_{{\mathrm{minl}}_{{\mathrm{P}}_{\mathrm{BIO},n}}}\\ \text{(10)}& {\displaystyle }& {\displaystyle }-F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_{\mathrm{HUM},n}}^{\mathrm{sorp}},\text{(11)}& {\displaystyle }{P}_{{\mathrm{O}}_{\mathrm{s}}}& {\displaystyle }=\sum _{n=\mathrm{1}}^N{P}_{{\mathrm{BIO}}_n}+\sum _{n=\mathrm{1}}^N{P}_{{\mathrm{HUM}}_n}.\end{array}$$

The description of the inorganic sorbed P pool (P_inorg-sorp) follows Wang et al. (2007) and is represented as the difference between the input flux of inorganic sorption ($F_{{\mathrm{P}}_{\mathrm{in}}}^{\mathrm{sorp}}$) (Eq. 37) and output fluxes of inorganic desorption $(F_{{\mathrm{P}}_{\mathrm{in}}}^{\mathrm{desorp}}$) (Eq. 38) and occluded P$(F_{\mathrm{P}}^{\mathrm{occ}}$) (Eq. 39), with the fluxes expressed (in kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\mathrm{d}{P}_{\text{inorg-sorp}}}{\mathrm{d}t}}=\sum _{n=\mathrm{1}}^N{F}_{{\mathrm{P}}_{{\mathrm{in}}_n}}^{\mathrm{sorp}}-\sum _{n=\mathrm{1}}^N{F}_{{\mathrm{P}}_{{\mathrm{in}}_n}}^{\mathrm{desorp}}-\sum _{n=\mathrm{1}}^N{F}_{{\mathrm{P}}_n}^{\mathrm{occ}}.\end{array}$$

The description of the occluded (P_occ) P pool follows Wang et al. (2007) and Hou et al. (2019) and is represented as the sum of input fluxes of occluded P from both organic $(F_{\mathrm{P}}^{\text{or-occ}}$) (Eq. 42) and inorganic P pools $(F_{\mathrm{P}}^{\mathrm{occ}}$) expressed (in kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{occ}}}{\mathrm{d}t}}=\sum _{n=\mathrm{1}}^N{F}_{{\mathrm{P}}_n}^{\mathrm{occ}}+\sum _{n=\mathrm{1}}^N{F}_{{\mathrm{P}}_n}^{\text{or-occ}}.\end{array}$$

The description of the organic sorbed P pool (P_org-sorp) follows Wang et al. (2007) and is represented as the difference between the input flux of organic sorption $(F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{sorp}}$) and output fluxes of organic desorption $(F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{desorp}}$) and occluded P$(F_{{\mathrm{P}}_n}^{\mathrm{occ}}$) with the fluxes expressed (in kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}{P}_{\text{org-sorp}}}{\mathrm{d}t}}=\sum _{n=\mathrm{1}}^N{F}_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{sorp}}-\sum _{n=\mathrm{1}}^N{F}_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{desorp}}-\sum _{n=\mathrm{1}}^N{F}_{{\mathrm{P}}_n}^{\text{or-occ}}.\end{array}$$

P from the parent material (P_pm) pool follows Wang et al. (2007) and depends on the weathering flux ($F_{\mathrm{P}}^{\mathrm{w}}$) (Eq. 43) expressed (in kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{d}{P}_{\mathrm{pm}}}{\mathrm{d}t}}=-\sum _{n=\mathrm{1}}^N{F}_{{\mathrm{P}}_n}^{\mathrm{w}}.\end{array}$$

2.2.2 C and P fluxes

NPP in JULES is calculated as the difference between GPP and autotrophic respiration. In JULES-CNP, potential NPP represents the amount of C available for tissue growth (C density increase) on a unit area and the associated spreading (vegetation cover increase as a result of reproduction and recruitment), i.e. increase to the area covered by the vegetation type, assuming no nutrient limitation. The reported NPP in the literature often includes other C fluxes related to the exudates, production of volatiles, and non-structural carbohydrates (Walker et al., 2021; Chapin et al., 2011; Malhi et al., 2009), which are challenging to measure (Malhi et al., 2011). Therefore, actual NPP is for our purposes equal to biomass production (BP) and is calculated as potential NPP minus excess C (lost to the plant through autotrophic respiration), with the latter being the C that cannot be used to grow new plant tissue due to insufficient plant nutrient supply. Hence, if the system is limited by the availability of N and/or P, NPP will be adjusted to match the growth that can be supported with the limited N or P supply, with any excess carbohydrate lost through excess C.

The total excess C term (ψ_t) (kg C m⁻² yr⁻¹) is calculated as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{\psi }}_{\mathrm{t}}={\mathit{\psi }}_{\mathrm{g}}+{\mathit{\psi }}_{\mathrm{s}},\end{array}$$

where ψ_g and ψ_s are the excess C fluxes due to growth (g) and spread (s) and are assumed to be rapidly respired by plants.

Therefore, BP is calculated as the difference between potential NPP (Π_c) and total excess C:

$$\begin{array}{}\text{(17)}& \mathrm{BP}={\mathrm{\Pi }}_c-{\mathit{\psi }}_{\mathrm{t}}.\end{array}$$

The litter production in JULES before limitation is estimated as follows:

$$\begin{array}{}\text{(18)}& F_{{\mathrm{C}}_n}^{\mathrm{lit}}={\mathit{\gamma }}_{\mathrm{leaf}}{C}_{\mathrm{leaf}}+{\mathit{\gamma }}_{\mathrm{root}}{C}_{\mathrm{root}}+{\mathit{\gamma }}_{\mathrm{wood}}{C}_{\mathrm{wood}},\end{array}$$

where γ is a temperature-dependent turnover rate representing the phenological state (Clark et al., 2011). P limitation is applied to the C litter production similar to the N scheme of JULES (JULES-CN) (Wiltshire et al., 2021). In JULES-CN the N limitation effect on the litter production is captured by estimating the available C for litter production as a difference between the NPP and excess C (Wiltshire et al., 2021).

Similar to other P-enabled models (Goll et al., 2017; Yang et al., 2014), JULES-CNP follows the same structure as its N model component. Descriptions of the plant P and N demand follow Wang et al. (2007) and are represented by the sum of demand (∅_t) to sustained growth (P-related: (${\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}}$); N-related: $\left({\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}}\right))$ and to sustain vegetation spreading (to increment PFT fractional coverage) (P-related: $({\mathrm{\varnothing }}_{{\mathrm{S}}_{\mathrm{P}}}$); N-related: $\left({\mathrm{\varnothing }}_{{\mathrm{S}}_{\mathrm{N}}}\right))$ and are expressed in different units (P-related values are expressed in kg P m⁻² yr⁻¹, while N-related values are expressed in kg N m⁻² yr⁻¹). The total demand for growth (∅_g) and spreading (∅_s) is controlled by the dominant demand between P (${\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}}$) and N $\left({\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}}\right)$ as follows:

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathrm{\varnothing }}_{\mathrm{t}}& {\displaystyle }={\mathrm{\varnothing }}_{\mathrm{g}}+{\mathrm{\varnothing }}_{\mathrm{s}},\text{(20)}& {\displaystyle }{\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}}& {\displaystyle }={\displaystyle \frac{P_{\mathrm{p}}}{C_V}}\left({\mathrm{\Pi }}_c-{\displaystyle \frac{\mathrm{d}{C}_v}{\mathrm{d}t}}-{\mathit{\psi }}_{\mathrm{g}}\right),\text{(21)}& {\displaystyle }{\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{P}}}& {\displaystyle }={\displaystyle \frac{P_{\mathrm{p}}}{C_V}}\left({\mathrm{\Pi }}_c-{\displaystyle \frac{\mathrm{d}{C}_v}{\mathrm{d}t}}-{\mathit{\psi }}_{\mathrm{s}}\right),\text{(22)}& {\displaystyle }{\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}}& {\displaystyle }={\displaystyle \frac{N_v}{C_V}}\left({\mathrm{\Pi }}_c-{\displaystyle \frac{\mathrm{d}{C}_v}{\mathrm{d}t}}-{\mathit{\psi }}_{\mathrm{g}}\right),\text{(23)}& {\displaystyle }{\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{N}}}& {\displaystyle }={\displaystyle \frac{N_v}{C_V}}\left({\mathrm{\Pi }}_c-{\displaystyle \frac{\mathrm{d}{C}_v}{\mathrm{d}t}}-{\mathit{\psi }}_{\mathrm{s}}\right),\text{(24)}& {\displaystyle }{\mathrm{\varnothing }}_{\mathrm{g}}& {\displaystyle }=\left\{\begin{array}{ll}{\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}};& {\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}}\times \frac{C_V}{P_{\mathrm{p}}}>{\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}}\times \frac{C_V}{N_v}\\ {\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}};& {\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}}\times \frac{C_V}{N_v}>{\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}}\times \frac{C_V}{P_{\mathrm{p}}}\end{array},\right.\text{(25)}& {\displaystyle }{\mathrm{\varnothing }}_{\mathrm{s}}& {\displaystyle }=\left\{\begin{array}{ll}{\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{P}}};& {\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{P}}}\times \frac{C_V}{P_{\mathrm{p}}}>{\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{N}}}\times \frac{C_V}{N_v}\\ {\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{N}}};& {\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{N}}}\times \frac{C_V}{N_v}>{\mathrm{\varnothing }}_{{\mathrm{s}}_{\mathrm{P}}}\times \frac{C_V}{P_{\mathrm{p}}}\end{array},\right.\end{array}$$

where $\frac{P_{\mathrm{p}}}{C_V}$ is the inverse of whole plant C:P ratio, $\frac{N_v}{C_V}$ is inverse plant C:N ratio, $\frac{\mathrm{d}{C}_v}{\mathrm{d}t}$ is rate of change in plant C (see Clark et al., 2011, for more detail), Π_c is nutrient-unlimited (or potential) NPP (kg C m⁻² yr⁻¹), ψ_g is excess C due to either P or N limitation for plant growth (kg C m⁻² yr⁻¹), and ψ_s is excess C due to either P or N limitation for vegetation spreading (kg C m⁻² yr⁻¹).

Equations 20 and 22 are solved by first setting ψ_g=0.0 to find the total plant P (Eq. 20) and N demand (Eq. 22). If the P and N demand for growth are less than the available P and N and fractional coverage (λ) (NPP fraction used for fractional cover increment; for more details, see Wiltshire et al., 2021) at the considered time step Δt, there is no limitation to growth (i.e. ${\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}}<\frac{\left(\mathrm{1}-\mathit{\lambda }\right)P_{\mathrm{avail}}}{\mathrm{\Delta }t};{\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}}<\frac{\left(\mathrm{1}-\mathit{\lambda }\right)N_{\mathrm{avail}}}{\mathrm{\Delta }t}$). Where there is limited P and/or N availability, the uptake is equal to the available P and N, i.e. $({\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{P}}}=\frac{\left(\mathrm{1}-\mathit{\lambda }\right)P_{\mathrm{avail}}}{\mathrm{\Delta }t};{\mathrm{\varnothing }}_{{\mathrm{g}}_{\mathrm{N}}}=\frac{\left(\mathrm{1}-\mathit{\lambda }\right)N_{\mathrm{avail}}}{\mathrm{\Delta }t})$, and the plant growth that cannot be achieved due to nutrient constraints will be deducted from potential NPP, here termed excess C term (ψ_g), to give an actual NPP value. Following Wiltshire et al. (2021), we assume that excess C is respired by the plant.

Similarly, in order to estimate the P and N demand for spreading (Eq. 21 and 23), initially the excess C from spreading is set to 0.0 ($\mathit{\psi }{}_{\mathrm{s}}=\mathrm{0}.$0), i.e. under the assumption that there is no nutrient limitation. If the P and N demand for spreading are lower than the available P and N and fractional coverage (λ) $({\mathrm{\varnothing }}_{{\mathrm{S}}_{\mathrm{P}}}<\frac{\left(\mathrm{1}-\mathit{\lambda }\right)P_{\mathrm{avail}}}{\mathrm{\Delta }t};{\mathrm{\varnothing }}_{{\mathrm{S}}_{\mathrm{N}}}<\frac{\left(\mathrm{1}-\mathit{\lambda }\right)N_{\mathrm{avail}}}{\mathrm{\Delta }t})$, then there is no limitation on spreading, and in cases of limited P and N availability, the uptake is equal to the available P and N, i.e. $({\mathrm{\varnothing }}_{{\mathrm{S}}_{\mathrm{P}}}=\frac{\left(\mathrm{1}-\mathit{\lambda }\right)P_{\mathrm{avail}}}{\mathrm{\Delta }t};{\mathrm{\varnothing }}_{{\mathrm{S}}_{\mathrm{N}}}=\frac{\left(\mathrm{1}-\mathit{\lambda }\right)N_{\mathrm{avail}}}{\mathrm{\Delta }t})$, and the excess C for spread (ψ_s) is subtracted from potential NPP.

Plant P uptake ($F_{\mathrm{p}}^{\mathrm{up}}$) (arrow a in Fig. 1) is estimated based on the P demand for growth and spreading (∅_t) and the root uptake capacity (u^max) (kg P kg⁻¹ C yr⁻¹) as follows:

$$\begin{array}{}\text{(26)}& F_{\mathrm{p}}^{{\mathrm{up}}_n}=\left\{\begin{array}{ll}{\mathrm{\varnothing }}_{\mathrm{t}};& {\mathrm{\varnothing }}_{\mathrm{t}}\le u^{max}\\ u^{max};& {\mathrm{\varnothing }}_{\mathrm{t}}>u^{max}\end{array}.\right.\end{array}$$

Plant P uptake ($F_{\mathrm{p}}^{\mathrm{up}}$) varies spatially depending on the root uptake capacity (u^max), following the method of Goll et al. (2017). Therefore, in regions with limited P supply, the plant P uptake is limited to the u^max and consequently impacts the excess C and BP.

The root uptake capacity depends on the maximum root uptake capacity (v_max) (kg P kg⁻¹ C yr⁻¹), root depth (d_root), the concentration of inorganic P at different soil depths (P_in), and a half-saturation term at which half of the maximum uptake capacity is reached using inorganic P at different soil depths (P_in), a scaling uptake ratio (K_p) (µmol P L⁻¹), unit conversion (C_f) (1 kg P⁻¹), and soil moisture (θ) (L m⁻²) as follows:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{u}^{max}& =v_{\max }\times d_{\mathrm{root}}\times \sum _{n=\mathrm{1}}^N{P}_{{\mathrm{in}}_n}\\ & \times \left({\displaystyle \frac{\mathrm{1}}{\sum _{n=\mathrm{1}}^N{P}_{{\mathrm{in}}_n}+c_{\mathrm{f}}\times K_{\mathrm{p}}\times {\mathit{\theta }}_n}}\right),\end{array}\end{array}$$

Description of the litter production of P $(F_{{\mathrm{P}}_n}^{\mathrm{lit}}$) (arrow b in Fig. 1) follows JULES-CN, as in Wiltshire et al. (2021), and is calculated based on the litter flux of C (kg C m⁻² yr⁻¹) using leaf, root, and wood turnovers (yr⁻¹) and through the vegetation dynamics due to large-scale disturbance and litter production density as follows:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{F}_{{\mathrm{P}}_n}^{\mathrm{lit}}& =\left({\displaystyle \frac{\left(\mathrm{1}-k_{\mathrm{leaf}}\right){\mathit{\gamma }}_{\mathrm{leaf}}{C}_{\mathrm{leaf}}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{leaf}}}}\right)\\ & +\left({\displaystyle \frac{\left(\mathrm{1}-k_{\mathrm{root}}\right){\mathit{\gamma }}_{\mathrm{root}}{C}_{\mathrm{root}}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{root}}}}\right)+\left({\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{wood}}{C}_{\mathrm{wood}}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{stem}}}}\right)\end{array},\end{array}$$

where λ is the leaf, root, and stem re-translocation (at daily time step) coefficient (Zaehle and Friend, 2010; Clark et al., 2011) and the related C:P ratios for P fraction and γ are a temperature-dependent turnover rate representing the phenological state (Clark et al., 2011).

The decomposition of litter (dec^lit) (arrow c in Fig. 1) depends on soil respiration (R) (kg C m⁻² yr⁻¹) and the litter C:P ratio (C:P_lit) at each soil layer (n) as follows:

$$\begin{array}{}\text{(29)}& {\mathrm{dec}}_{\mathrm{P}}^{\mathrm{lit}}={\displaystyle \frac{\sum _{n=\mathrm{1}}^N{R}_n}{{\mathrm{C}:\mathrm{P}}_{\mathrm{lit}}}},\end{array}$$

where the C:P_lit is calculated based on the litter C pool (DPM and RPM) (lit^C) (kg C m⁻² yr⁻¹) and litter P pool ($P_{{\mathrm{O}}_{\mathrm{l}}}$) as follows:

$$\begin{array}{}\text{(30)}& {\mathrm{C}:\mathrm{P}}_{\mathrm{lit}}={\displaystyle \frac{\sum _{n=\mathrm{1}}^N{\mathrm{lit}}_n^{\mathrm{C}}}{P_{{\mathrm{O}}_{{\mathrm{l}}_n}}}}.\end{array}$$

The mineralized ($F_{{\mathrm{minl}}_{\mathrm{P}}}$) (arrow d in Fig. 1) and immobilized ($F_{{\mathrm{immob}}_{\mathrm{P}}}$) (arrow e in Fig. 1) P fluxes are calculated based on C mineralization and immobilization, C:P ratios of plant (i) (DPM $/$ RPM) (C:P_plant) and soil (HUM $/$ BIO) (C:P_soil), soil pool potential respiration ($R_{{\mathrm{POT}}_i}$) (kg C m⁻² yr⁻¹), and the respiration partitioning fraction (resp_frac) as follows:

$$\begin{array}{}\text{(31)}& {\displaystyle }{F}_{{\mathrm{minl}}_{{\mathrm{P}}_n}}& {\displaystyle }={\displaystyle \frac{\sum _{n=\mathrm{1}}^N{R}_{{\mathrm{POT}}_{i,n}}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{plant}}}},\text{(32)}& {\displaystyle }{F}_{{\text{immob}}_{{\mathrm{P}}_n}}& {\displaystyle }={\displaystyle \frac{\sum _{n=\mathrm{1}}^N{R}_{i,n}\times \text{resp\_frac}}{{\mathrm{C}:\mathrm{P}}_{\mathrm{soil}}}}.\end{array}$$

The soil respiration from each soil layer (R_i,n) is estimated from potential soil respiration ($R_{{\mathrm{POT}}_{i,n}}$) for the DPM and RPM pools and the litter decomposition rate modifier ($F_{{\mathrm{P}}_n}$) as follows:

$$\begin{array}{}\text{(33)}& R_{i,n}=R_{{\mathrm{POT}}_{i,n}}\times F_{{\mathrm{P}}_n},\end{array}$$

where the description of $F_{{\mathrm{P}}_n}$ for P pools ($F_{P_{{\mathrm{P}}_n}}$) follows Wang et al. (2007) and is estimated based on the soil pool (BIO $/$ HUM) mineralization (${\mathrm{minl}}_{{\text{P-BIO}}_n}$, ${\mathrm{minl}}_{{\text{P-HUM}}_n}$) and immobilization (${\mathrm{immob}}_{{\text{P-BIO}}_n}$, ${\mathrm{immob}}_{{\text{P-HUM}}_n}$) (in kg P m⁻² yr⁻¹), soil inorganic P ($P_{{\mathrm{inorg}}_n}$) (in kg P m⁻²), and litter pool (DPM $/$ RPM) demand (in kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(34)}& \begin{array}{rl}& F_{P_{{\mathrm{P}}_n}}=\\ & {\displaystyle \frac{({\mathrm{minl}}_{{\text{P-BIO}}_n}+{\mathrm{minl}}_{{\text{P-HUM}}_n}-{\mathrm{immob}}_{{\text{P-BIO}}_n}-{\mathrm{immob}}_{{\text{P-HUM}}_n})+P_{{\mathrm{inorg}}_n}}{{\mathrm{DEM}}_{{\mathrm{DPM}}_n}+{\mathrm{DEM}}_{{\mathrm{RPM}}_n}}}\end{array}.\end{array}$$

The net demand associated with decomposition of litter pools (DEM_k,n) represents the P required by microbes that convert DPM and RPM into BIO and HUM. The limitation due to insufficient P availability is estimated based on the potential mineralization (minl_p-pot) and immobilization (immob_p-pot) (in kg P m⁻² yr⁻¹) of pools (k) as follows:

$$\begin{array}{}\text{(35)}& {\mathrm{DEM}}_{k,n}={\mathrm{immob}}_{\text{p-pot},k}-{\mathrm{minl}}_{\text{p-pot},k}.\end{array}$$

The $F_{{\mathrm{P}}_n}$ estimated for N pools ($F_{{\text{P}}_{{\mathrm{N}}_n}}$) follows the same formulation as P (see Wiltshire et al., 2021, for further details), and the $F_{{\mathrm{P}}_n}$ is estimated based on a higher-rate modifier between N and P as follows:

$$\begin{array}{}\text{(36)}& F_{{\mathrm{P}}_n}=\left\{\begin{array}{ll}{F}_{{\mathrm{P}}_{{\mathrm{P}}_n}};& F_{{\mathrm{P}}_{{\mathrm{P}}_n}}>F_{{\mathrm{P}}_{{\mathrm{N}}_n}}\\ F_{{\mathrm{P}}_{{\mathrm{N}}_n}};& F_{{\mathrm{P}}_{{\mathrm{N}}_n}}>F_{{\mathrm{P}}_{{\mathrm{P}}_n}}\end{array}.\right.\end{array}$$

Description of the fluxes of adsorption ($F_{{\mathrm{P}}_{{\mathrm{in}}_n}}^{\mathrm{sorp}}$) (arrow e in Fig. 1) and desorption ($F_{{\mathrm{P}}_{{\mathrm{in}}_n}}^{\mathrm{desorp}}$) (arrow f in Fig. 1) of inorganic P (in kg P m⁻² yr¹) follow Wang et al. (2010) and are calculated based on soil inorganic ($P_{{\mathrm{in}}_n}$) and sorbed inorganic ($P_{{\text{inorg-sorbed}}_n}$) P pools, inorganic adsorption (K_sorp-in) and desorption (K_desorp-in) coefficients (kg P m⁻² yr⁻¹), and maximum sorbed inorganic (P_in-max) (kg P m⁻²) as follows:

$$\begin{array}{}\text{(37)}& {\displaystyle }{F}_{{\mathrm{P}}_{{\mathrm{in}}_n}}^{\mathrm{sorp}}& {\displaystyle }=P_{{\mathrm{in}}_n}\times K_{\text{sorp-in}}\times {\displaystyle \frac{\left(P_{{\text{in-max}}_n}-P_{{\text{inorg-sorbed}}_n}\right)}{P_{{\text{in-max}}_n}}},\text{(38)}& {\displaystyle }{F}_{{\mathrm{P}}_{{\mathrm{in}}_n}}^{\mathrm{desorp}}& {\displaystyle }=P_{{\text{inorg-sorbed}}_n}\times K_{\text{desorp-in}}.\end{array}$$

Description of the occluded inorganic P flux ($F_{{\mathrm{P}}_n}^{\mathrm{occ}}$) (arrow g in Fig. 1) follows Wang et al. (2007) and Hou et al. (2019) and is calculated based on sorbed inorganic P pool and P occlusion rate (K_occ) (kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(39)}& F_{{\mathrm{P}}_n}^{\mathrm{occ}}=P_{{\text{inorg-sorbed}}_n}\times K_{\mathrm{occ}}.\end{array}$$

Description of the fluxes of adsorption ($F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{sorp}}$) (arrow h in Fig. 1) and desorption ($F_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{desorp}}$) (arrow i in Fig. 1) of organic P follow Wang et al. (2010) and are calculated based on soil organic and sorbed organic P pools and organic adsorption (K_sorp-or) (kg P m⁻² yr⁻¹), desorption (K_desorp-or) coefficients (kg P m⁻² yr⁻¹), and maximum sorbed organic (P_org-max) (which corresponds to the sorbed soil P saturation, thus modifying the sorption rate) (kg P m⁻²) as follows:

$$\begin{array}{}\text{(40)}& {\displaystyle }{F}_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{sorp}}& {\displaystyle }=P_{{\mathrm{O}}_{{\mathrm{S}}_n}}\times K_{\text{sorp-or}}\times {\displaystyle \frac{\left(P_{{\text{or-max}}_n}-P_{{\text{org-sorbed}}_n}\right)}{P_{{\text{or-max}}_n}}},\text{(41)}& {\displaystyle }{F}_{{{\mathrm{P}}_{{\mathrm{O}}_{\mathrm{S}}}}_n}^{\mathrm{desorp}}& {\displaystyle }=P_{{\text{org-sorbed}}_n}\times K_{\text{desorp-or}}.\end{array}$$

Description of the occluded organic P flux ($F_{{\mathrm{P}}_n}^{\text{or-occ}}$) (kg P m⁻² yr⁻¹) (arrow j in Fig. 1) follows Wang et al. (2007) and Hou et al. (2019) and is calculated based on sorbed organic P pool ($P_{{\text{org-sorbed}}_n}$) and P occlusion rate (K_occ) (kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(42)}& F_{{\mathrm{P}}_n}^{\text{or-occ}}=P_{{\text{org-sorbed}}_n}\times K_{\mathrm{occ}}.\end{array}$$

Description of the P flux from weathered parent material ($F_{{\mathrm{P}}_n}^{\mathrm{w}}$) (arrow k in Fig. 1) follows Wang et al. (2007) and is calculated based on amount of P in the parent material (P_pm) and P weathering rate (K_w) (kg P m⁻² yr⁻¹) as follows:

$$\begin{array}{}\text{(43)}& F_{{\mathrm{P}}_n}^{\mathrm{w}}=P_{{\mathrm{pm}}_n}\times K_{\mathrm{w}}.\end{array}$$

Description of P diffusion between soil layers ($F_{{\mathrm{D}}_n}$) (expressed in kg P m⁻² yr⁻¹) (arrow l in Fig. 1) follows Goll et al. (2017) and is calculated following Fick's second law as a function of the diffusion coefficient (Dz) in m² s⁻¹, the concentration of inorganic P at different soil depths (P_in) in kg P m⁻², the distance (z) between the midpoints of soil layers in metres, and a seconds to years unit conversion (Yr):

$$\begin{array}{}\text{(44)}& F_{{\mathrm{D}}_n}={\displaystyle \frac{\partial }{\partial z}}\left(D_{z_n}{\displaystyle \frac{\partial P_{{\mathrm{s}}_n}}{\partial z}}\right)\times \mathrm{Yr}.\end{array}$$

Table 1 Model variables.

Table 2 P model parameters.

3 Study sites

This study primarily uses data from two close together sites in the central Amazon rainforest in Manaus, Brazil. The main site, from here on termed “study site” ($\mathrm{2}{}^{\circ }{\mathrm{35}}^{\prime }{\mathrm{21.08}}^{\prime \prime }$ S, $\mathrm{60}{}^{\circ }{\mathrm{06}}^{\prime }{\mathrm{53.63}}^{\prime \prime }$ W) (Lugli et al., 2020), is for model development and evaluation. The second site is the Manaus K34 flux site ($\mathrm{2}{}^{\circ }{\mathrm{36}}^{\prime }{\mathrm{32.67}}^{\prime \prime }$ S, $\mathrm{60}{}^{\circ }{\mathrm{12}}^{\prime }{\mathrm{33.48}}^{\prime \prime }$ W), which provides meteorological station data for running the model but also provides data for model evaluation. Our study site is the main lowland tropical forest site maintained by the National Institute for Amazon Research (INPA). Research at this site focuses on projects combining experimental approaches (Keller et al., 2004; Malhi et al., 2009) with modelling (Lapola and Norby, 2014). We use detailed novel soil and plant P pool data from the study site (Lugli et al., 2021, 2020) for model parameterization and calibration and carbon stock data for model validation. The study site has a very similar forest, geomorphology, soil chemistry, and species composition to the well-known and well-studied K34 flux site (Araújo et al., 2002). The average reported annual precipitation is 2431 mm yr⁻¹, with a monthly range of 95 to 304 mm per month, and the average temperature is 26 ^∘C (Araújo et al., 2002). The soil type at this site is Geric ferralsol with high clay content and weathering activities (Malhi et al., 2004).

In addition to the study site, we use data from other P-limited locations from the Amazon, Panama, and Hawaii (Table 3) for model evaluation. Old-growth forest sites in the Amazon are located across a fertility gradient from west to east (AGP-01, SA3, CAX) where detailed C cycle measurements are available (Aragão et al., 2009). The site in Panama is located on the Gigante Peninsula in the Barro Colorado Nature Reserve and is a 200-year-old semi-deciduous rainforest (Wright et al., 2011) growing on Oxisols developed on Miocene basalt (Dieter et al., 2010) with a topsoil that has a dominant clay texture (Turner and Condron, 2013). It has been the location of a long-term running nutrient fertilization experiment since 1998 (Mirabello et al., 2013). The site in Hawaii (Hawaii Kokee) is a P-limited chronosequence that developed on the 4-million-year-old Oxisols soil (Vitousek, 2004) and has a long-term fertilization experiment. Site information is provided in Table 3.

Table 3 Test site name, location, and climate characteristics.

3.1 Model parameterization, calibration, and evaluation at study site

We use observations from the four control plots of the study site to parameterize, calibrate, and evaluate different processes in JULES (Table 4). The observations were collected at four soil depths and processed using the Hedley sequential fractionation (Hedley et al., 1982; Quesada et al., 2010). Observed leaf mass per area (LMA) and leaf N and leaf P estimated from fresh leaves were used as input parameters to JULES to estimate photosynthetic capacity and respiration parameters. JULES vn5.5 (JULES-CN in this study) estimates V_cmax (µmol m⁻² s⁻²) based on Kattge et al. (2009) using foliar N concentrations in area basis (nleaf), as follows:

$$\begin{array}{}\text{(45)}& V_{\mathrm{cmax}}=v_{\mathrm{int}}+v_{\mathrm{sl}}\cdot \mathrm{nleaf},\end{array}$$

where v_int is the estimated intercept and v_sl is the slope of the linear regression derived for the V_cmax estimation. We incorporated an additional P dependency on the estimation of V_cmax following Walker et al. (2014) as follows:

$$\begin{array}{}\text{(46)}& \begin{array}{rl}\ln \left(V_{\mathrm{cmax}}\right)& =\mathrm{3.946}+\mathrm{0.921}\ln \left(N\right)+\mathrm{0.121}\ln \left(P\right)\\ & +\mathrm{0.282}\ln \left(N\right)\ln \left(P\right),\end{array}\end{array}$$

where N and P are foliar concentrations on an area basis.

Implementation of Eq. (46) resulted in higher V_cmax than in the original version of JULES. A higher V_cmax predicted higher leaf and plant respiration (Eq. 47). Constrained by observations of NPP and plant respiration at the study site, we modified one of the most uncertain parameters in the description of plant respiration (f_dr) (Eq. 47), namely the scale factor for leaf dark respiration (R_d), as follows:

$$\begin{array}{}\text{(47)}& R_{\mathrm{d}}=f_{\mathrm{dr}}{V}_{\mathrm{cmax}}.\end{array}$$

The default value is 0.01 (Clark et al., 2011), and for JULES-CNP simulations at our study site it was modified to 0.005.

Observations of above-ground biomass were used to calibrate the non-PFT-dependent allometric relationships in JULES (Clark et al., 2011) (Eqs. 48–50) for leaf, root, and wood C. Specifically, the a_wl parameter (Eq. 50) was modified from 0.65 to 1.204 to better match tropical forest allometry:

$$\begin{array}{}\text{(48)}& {\displaystyle }{C}_{\mathrm{leaf}}& {\displaystyle }={\mathit{\sigma }}_{\mathrm{l}}{L}_{\mathrm{b}},\text{(49)}& {\displaystyle }{C}_{\mathrm{root}}& {\displaystyle }=C_{\mathrm{leaf}},\text{(50)}& {\displaystyle }{C}_{\mathrm{stem}}& {\displaystyle }=a_{\mathrm{wl}}{L}_{\mathrm{b}}^{b_{\mathrm{wl}}},\end{array}$$

where L_b is the balanced (or seasonal maximum) leaf area index (LAI, m² m⁻²), σ_l is the specific leaf density (kg C m⁻² per unit LAI), a_wl is the allometric coefficient (kg C m⁻²), and b_wl is the allometric exponent.

Note that JULES-CNP uses the C₃ and C₄ photosynthesis model from Collatz et al. (1991, 1992), which does not include estimation of J_max.

JULES-CNP has fixed stoichiometry and C:P ratios of leaf and root (measured) and wood (estimated from fresh coarse wood, Lugli, 2013), which were taken from the study site and prescribed in JULES to simulate P dynamics in the plant. The following belowground data were used to represent various soil P pools: resin and bicarbonate inorganic P (inorganic P: P_in), organic bicarbonate P (organic P: $P_{{\mathrm{O}}_{\mathrm{S}}}$), NaOH organic P (sorbed organic P: P_org-sorp), NaOH inorganic P (sorbed inorganic P: P_inorg-sorp), residual P (occluded P: P_occ), and HCL P (parent material P: P_pm) (Table 4). The measurements were collected between 2017 and 2018 in control plots. All measurements were conducted in four soil layers (0–5, 5–10, 10–20, 20–30 cm). However, to be consistent with the JULES model soil layer discretization scheme, we defined four soil layers (0–10, 10–30, 30–100 and 100–300 cm) and used the average between 0 and 30 cm to compare against the measurement from the same depth for model evaluation.

Vegetation C stocks were derived based on tree diameter measurements at breast height, which are linked to allometric equations and wood density databases to estimate the C stored in each individual tree and then scaled to the plot (Chave et al., 2014).

The organic and inorganic soil P was assumed to be always at equilibrium with the relative sorbed pools (Wang et al., 2010). Thus, in order to cap P sorption and uptake capacity, the maximum sorption capacities ($P_{{\text{in-max}}_n}$, $P_{{\text{or-max}}_n}$ Eqs. 37 and 39) (adopted from Wang et al., 2007) were prescribed using maximum observed sorbed inorganic and organic P. Hence, the maximum sorption capacity defines the equilibrium state of sorbed and free-soil P. Moreover, despite the initial representation of the parent material pool in JULES and its depletion through weathering (Eq. 43), as the magnitude of changes in the occluded and parent material pools are insignificant over a short-term (20 years) simulation period (Vitousek et al., 1997), these two pools were prescribed using observations. The remaining parameters used to describe soil P fluxes (Eqs. 27–44) were prescribed using values from the literature (Table 4).

We used a combination of data from the study site and the nearby K34 site for model evaluation of C fluxes (GPP, NPP) and C pools (soil and vegetation C, leaf, root and wood C) with no calibration of plant and soil organic P pools and soil inorganic P pools included (Table 4).

Table 4 Observations from the study site (taken during 2017–2018) and from the Manaus site K34 used for model parameterization and evaluation.

3.2 Model parameterization and evaluation at test sites

JULES-CNP was parameterized using reported C:P ratios and maximum sorbed organic and inorganic P for each test site (Table 5) as follows.

Table 5 Additional test site data used for model parameterization.

^a C:P ratios from Wang et al. (2010). ^b Maximum sorbed P capacities from Yang et al. (2014). ^c Mirabello et al. (2013). ^d C:P ratios from Vitousek (2004).

Model evaluation at test sites was performed using observed NPP, litterfall, autotrophic respiration, biomass, and soil C pools taken from different sources. We used NPP and litterfall for the Amazon sites from Aragão et al. (2009) and for Gigante Peninsula from Chave et al. (2003), Hawaii Kokee NPP as reported in Goll et al. (2017), and litterfall as reported in Yang et al. (2014). Plant respiration was only available at two of Amazon sites (AGP and CAX) (Malhi et al., 2009). The biomass and soil C pools for two of the Amazon sites (CAX and SA3) are taken from Malhi et al. (2009), and biomass from AGP is taken from Jiménez et al. (2009). The Gigante Peninsula biomass is taken from Chave et al. (2003), soil C is taken from Turner et al. (2015), and the Hawaii Kokee C pools are taken as reported in Yang et al. (2014).

3.3 JULES simulations

JULES was first applied at the K34 flux tower site using observed meteorological forcing data from 1999–2019 (Fleischer et al., 2019) at half-hourly resolution. The following meteorological variables are needed to drive JULES (model inputs) (Best et al., 2011): atmospheric specific humidity (kg kg⁻¹), atmospheric temperature (K), air pressure at the surface (Pa), short and longwave radiation at the surface (W m⁻²), wind speed (m s⁻¹), and total precipitation (kg m⁻² s⁻¹). Furthermore, the averaged measured LAI from the study site was used to initialize the vegetation phenology module but was allowed to vary in subsequent prognostic calculations. Soil organic and inorganic sorbed P pools were initialized with study site observations. The JULES-CNP simulations were initialized following the same methodology as in Fleischer et al. (2019), and the spin-up from 1850 CE resulted in an equilibrium state (Fig. S1 in the Supplement). The spin-up was performed separately for three versions of JULES (C/CN/CNP) following the same procedure. Furthermore, the transient run was performed for the period 1851–1998 using time-varying CO₂ and N deposition fields. Finally, for the extended simulation period (1999–2019) two runs were performed, the first with ambient CO₂ concentrations and the second with elevated CO₂ concentrations.

We evaluate the impact of including a P cycle in JULES using three model configurations (JULES-C, JULES-CN, and JULES-CNP). We apply JULES in all three configurations using present-day climate under both ambient CO₂ and eCO₂. Ambient and eCO₂ were prescribed following Fleischer et al. (2019), with present-day CO₂ based on global monitoring stations and an abrupt (step) increase in atmospheric CO₂ of +200 ppm on the onset of the transient period (i.e. 1999). However, the comparison period is limited to the 2017–2018 period for which the P measurements are available.

We compare simulated C fluxes (GPP, NPP, litterfall C), C stocks (total vegetation, fine root, leaf, wood, soil), and the CO₂ fertilization effect across model configurations. The CO₂ fertilization effect $\left({{\mathrm{CO}}_{\mathrm{2}}}_{\text{fert-eff}}\right)$ (Eq. 51) is calculated based on simulated vegetation C under ambient (VegC (aCO₂)) and eCO₂ (VegC (eCO₂)) as follows:

$$\begin{array}{}\text{(51)}& {{\mathrm{CO}}_{\mathrm{2}}}_{\text{fert-eff}}={\displaystyle \frac{\left(\mathrm{VegC}\left({\mathrm{eCO}}_{\mathrm{2}}\right)-\mathrm{VegC}\left({\mathrm{aCO}}_{\mathrm{2}}\right)\right)\times \mathrm{100}}{\mathrm{VegC}\left({\mathrm{aCO}}_{\mathrm{2}}\right)}}.\end{array}$$

Furthermore, the net biomass increases due to the CO₂ fertilization effect (ΔC_veg) is estimated as follows:

$$\begin{array}{}\text{(52)}& \mathrm{\Delta }{\mathrm{C}}_{\mathrm{veg}}=\mathrm{\Delta }\mathrm{BP}-\mathrm{\Delta }\text{litterfall\hspace{0.17em} C}.\end{array}$$

We studied the water use efficiency (WUE) (Eq. 53) at half-hourly time steps, which was then aggregated per month as one of the main indicators of GPP changes (Xiao et al., 2013), and soil moisture content (SMCL), as one of the main controllers of maximum uptake capacity (Eq. 27), in order to gain a better understanding of the changes in GPP, P demand and uptake, and excess C fluxes.

$$\begin{array}{}\text{(53)}& \mathrm{WUE}=\mathrm{GPP}/\mathrm{Transpiration}\end{array}$$

Moreover, we also estimated the carbon use efficiency (CUE) as an indicator of the required C for the growth (Bradford and Crowther, 2013) as follows:

$$\begin{array}{}\text{(54)}& \mathrm{CUE}=\mathrm{BP}/\mathrm{GPP}.\end{array}$$

We use JULES-CNP to evaluate the extent of P limitation under ambient and eCO₂ at this rainforest site in the central Amazon. P limitation is represented by the amount of C that is not used to grow new plant tissue due to insufficient P in the system (excess C) (Eq. 27). The excess C flux is highly dependent on the plant P and the overall P availability to satisfy demand. We also explore the distribution of the inorganic and organic soil P and their sorbed fraction within the soil layers and under ambient and eCO₂.

3.3.1 Model sensitivity

To test the sensitivity of the P- and C-related processes to individual model P parameters, six sets of simulations were conducted independently with modified plant C:P stoichiometry (Plant C:P: SENS1), P uptake scaling factor (K_P) (Kp: SENS2), inorganic (KP_ sorb_in: SENS3) and organic (KP_sorb_or: SENS4) P adsorption coefficients (K_sorp−or, K_sorp-in), and maximum inorganic (KP_sorb_in_max: SENS5) and organic (KP_sorb_or_max: SENS6) sorbed P (K_or-max, K_in-max). These values were prescribed to vary between ±50 % of the observed values, and their effect on C pools (plant and soil C) and fluxes (NPP and excess C) and P pools (plant, soil, and soil sorbed P) was assessed. As the derived model parameters from measurements have their own level of uncertainty, we took 50 % change to test these parameters at a reasonable degree. However, the occluded and weathered P pools are prescribed for this model application, and the occluded and weather P coefficients (the other two P-related model parameters) were not part of sensitivity tests.

Our model evaluation period is limited to the years 2017–2018 due to the P measurement availability. However, in order to compare it with 15 models studied by Fleischer et al. (2019), we also studied the response of GPP, NPP, and BP to eCO₂ for both initial (1999) and 15-year periods (between 1999 and 2013).

3.3.2 Simulations at test sites

To perform JULES (C, CN, CNP) simulations at test sites we extracted the meteorological input data to drive the model from a global dataset (CRU-NCEP) (Harris et al., 2014) by selecting the closest grid cell to each site when data were not available for a given site (Table 3). Soil texture ancillaries for each site were extracted from a global soil database (HWSD) (Nachtergaele et al., 2010). All simulations were initialized from a global JULES-CN run (Wiltshire et al., 2021) extracted for each site and further spun up for 2000 years over the 1980–2000 period for the three versions of JULES (C/CN/CNP). Finally, the transient (2000–2013) run was performed using the output of the spin-up for each site.

4 Results

4.1 Model application under ambient CO₂

4.1.1 Calibration of simulated soil P pools at study site

The maximum sorption capacities ($P_{{\text{in-max}}_n}$, $P_{{\text{or-max}}_n}$, Eqs. 37 and 40, respectively) were calibrated to the observed P pools. As a result, JULES-CNP could reproduce the measured soil P pools (Fig. 2 and Table 6). Simulated inorganic soil P and sorbed organic and inorganic soil P closely matched the observations (Table 7 and Fig. 2). However, simulated organic soil P overestimates the observations by 60 %.

Figure 2 Modelled vs. measured soil phosphorus pools under ambient CO₂ (for the soil depth of 0–30 cm). The black line represents the standard deviation.

Table 6 Observed and simulated phosphorus pools and fluxes. Occluded and weathered P pools were prescribed using the observed values (for 2017–2018).

4.1.2 Model evaluation

JULES-CNP could reproduce the plant and soil C (Fig. 2 and Table 7) and N pools and fluxes (Fig. S6 and Table 8) under ambient CO₂. Our results show that simulated GPP is within the range of measurement (3.02 kg C m⁻² yr⁻¹ modelled vs. 3–3.5 kg C m⁻² yr⁻¹ observed, Table 7).

Simulated NPP is close to the measured values (NPP: 1.14–1.31 kg C m⁻² yr⁻¹ observed vs. 1.26 kg C m⁻² yr⁻¹ modelled) with autotrophic respiration (RESP) also closely following the observations (1.98 kg C m⁻² yr⁻¹ observed vs. 1.81 kg C m⁻² yr⁻¹ modelled). Biomass production is estimated as the difference between NPP and the amount of C that is not fixed by plants due to the insufficient P in the system (excess C) (Eq. 27). The excess C flux depends on the plant P and the overall P availability to satisfy demand (Table 7). The simulated flux of excess C is 0.3 kg C m⁻² yr⁻¹ under ambient CO₂. In JULES-CNP this flux is subtracted from NPP in order to give the BP (Eq. 17) (Table 7). Our simulated litterfall overestimates the observations by 32 %; however, simulated vegetation and its components (fine root, leaf, and wood) and soil C stocks match the observations well (Table 7).

Table 7 Observed and simulated carbon pools and fluxes with JULES-CNP (between period 2017–2018).

4.1.3 Comparison of JULES-C, JULES-CN, and JULES-CNP under ambient CO₂ at the study site

We compare simulated C pools and fluxes from JULES-C, JULES-CN, and JULES-CNP (Fig. 3). There is no difference between C stocks and fluxes in simulations from JULES-C and JULES-CN, indicating that there is no N limitation at this tropical site in the CN simulations. However, simulated BP and litter flux of C by JULES-C and JULES-CN are higher than in JULES-CNP but also overestimate the observations: the litter flux of JULES-C and JULES-CN is 1.18 kg C m⁻² yr¹ and that of JULES-CNP is 0.91 kg C m⁻² yr¹, while the observations show a value of 0.69 kg C m⁻² yr¹. The BP of JULES-C and JULES-CN is 1.24 and that of JULES-CNP is 0.96, while the observations show values of 1.14–1.31  kg C m⁻² yr⁻¹. By including P cycling in JULES, an excess C flux of 0.3 kg C m⁻² yr⁻¹ is simulated, indicating a 24 % P limitation to BP at this site according to JULES-CNP, which represents a 29 % decrease in BP compared to JULES-C and JULES-CN. Consequently, the total vegetation C stock for models without P inclusion is higher than the CNP version (+3 % difference) due to the lack of representation of P limitation. The simulated soil C stock in JULES-C and JULES-CN is also higher than in the CNP version (JULES-C and JULES-CN at 13.93 kg C m⁻² yr⁻¹ vs. JULES-CNP at 13.18 kg C m⁻² yr⁻¹) and higher than the observations. Moreover, CUE in JULES-C and JULES-CN (Eq. 54) is higher than observations and JULES-CNP (JULES-C and JULES-CN: 0.38; JULES-CNP: 0.31; obs: 0.34 ± 0.1; note that this value is dimensionless).

Figure 3 JULES-C, JULES-CN, and JULES-CNP modelled vs. measured C pools (leaf, root, wood, vegetation, and soil C) (in kg C m⁻²) and fluxes (BP and litter C) (in kg C m⁻² yr⁻¹) and CUE under ambient CO₂. Note that CUE is unitless.

4.1.4 Model evaluation at test sites under ambient CO₂

Evaluation of JULES-C, JULES-CN, and JULES-CNP at five test sites against the observed C pools and fluxes demonstrates that the inclusion of P processes improved the simulation of C pools and fluxes across all test sites (Fig. 4). At all Amazon sites, JULES-C and JULES-CN overestimated BP compared to JULES-CNP, which estimated lower BP values that were closer to the measurements for AGP (JULES-C: +35 %; JULES-CN: +33 %; JULES-CNP: +21 %), CAX (JULES-C: +45 %; JULES-CN: +44 %; JULES-CNP: +7 %), and SA3 (JULES-C: +27 %; JULES-CN: +26 %; JULES-CNP: −23 %). Moreover, at Gigante Peninsula JULES-C and JULES-CN overestimated BP (+42 % and +40 %, respectively) and CNP slightly underestimated BP (−15 %). Furthermore, at the Hawaii Kokee site, all three versions of JULES underestimated the BP (JULES-C: −8 %; JULES-CN: −8 %; JULES-CNP: −32 %). The litterfall and respiration fluxes in JULES-CNP have decreased compared to JULES-C and JULES-CN, which overestimated both fluxes at all the test sites compared to the measurements. The litterfall flux comparisons show a significant overestimation using JULES-C and JULES-CN across all the tested sites. For the Amazon sites inclusion of P limitation reduced the litterfall flux but still overestimated it (AGP: +50 %; CAX: +24 %; SA3: +16 %), and at Gigante Peninsula and Hawaii Kokee litterfall flux was slightly underestimated (Gigante Peninsula: −9 %; Hawaii Kokee −19 %).

The respiration measurements were only available at two Amazon sites (CAX and SA3) at which inclusion of P limitation resulted in a well-estimated flux at both sites compared to JULES-C and JULES-CN (CAX: +38 %, +38 %, and −1 % for C, CN, and CNP, respectively; SA3: +38 %, +38 %, and −2 % for C, CN, and CNP, respectively).

The total vegetation biomass also reduced using JULES-CNP compared to the other versions and yielded closer values to the measurements across all the sites. However, except at the AGP site in which all three versions of JULES slightly underestimated the biomass (C: −1 %; CN: −1 %; CNP: −6 %), at the other test sites JULES-CNP estimated lower biomass pools compared to the other versions, which also overestimated total vegetation biomass.

Similarly, the soil C pool was overestimated prior to inclusion of P limitation in JULES at the test sites, and JULES-CNP estimated a closer value compared to the measurements (slight underestimation at CAX and SA3 sites of −5 % and −18 %, respectively, and close values at Gigante Peninsula and Hawaii Kokee of +3 % and +4 %, respectively).

Figure 4 Observed and simulated (JULES-C, JULES-CN, JULES-CNP) C fluxes and pools (averaged measurements: red points; SD: red arrows) and available observed P (dark red points and lines, reported in ppm) at test sites across the Amazon (AGP, SA03, CAX), Gigante Peninsula (Gig. Pen.), and Hawaii Kokee (Hawaii K.).

4.1.5 Model sensitivity

The results indicate that among all the corresponding C and P pools and fluxes, the excess C flux – which demonstrates P limitation to growth – shows the highest sensitivity to changes in C:P ratios (Fig. 5a), K_P (Fig. 5b), K_or-max (Fig. 5c), and K_in-max (Fig. 5d). A decrease in plant C:P results in a large increase in excess C. This is due to the higher plant P demand as a result of lower plant C:P ratios. An increase in the uptake factor and maximum sorbed organic and inorganic P also results in an increase in excess C. This is due to the higher uptake demand through higher uptake capacity (due to higher K_P) and lower available P for uptake due to higher organic and inorganic sorbed P (due to higher K_or-max, K_in-max). Since the total P in the system is lower than the plant demand, the uptake capacity, and sorbed P, higher P limitation is placed on growth (decreasing BP), which results in an increase in excess C and decrease in plant C, and soil C also decreases as a result of lower litter input (Fig. 5). Total soil P shows low sensitivity to changes in plant C:P and uptake factor but high sensitivity to maximum inorganic sorbed P. Moreover, sorbed P shows medium and high sensitivity to maximum organic and inorganic sorbed P, respectively (Fig. S5). Nevertheless, organic and inorganic P adsorption coefficients (K_sorp-or, K_sorp-in) show no sensitivity to C and P pools and fluxes. This is due to limiting the organic and inorganic P sorption terms so that they are only controlled by maximum sorption capacity, meaning that they are not affected by organic and inorganic adsorption coefficients.

Figure 5 Sensitivity of C and P pools to variation in key model parameters: prescribed tissue C:P (a), Kp (b), Kp_sorb_in (c), and Kp_sorb_or (d) under ambient CO₂.

4.2 Model application under elevated CO₂

4.2.1 Simulated plant and soil C and P pools and fluxes in JULES-CNP: eCO₂ vs. ambient CO₂

The eCO₂ simulation using JULES-CNP yields a higher GPP compared to the ambient CO₂ (0.83 kg C m⁻² yr⁻¹ increase) as a result of CO₂ fertilization. Moreover, due to the GPP increase, NPP and RESP also increased compared to ambient CO₂ (NPP: 0.49 kg C m⁻² yr⁻¹ increase; RESP: 0.3 kg C m⁻² yr⁻¹ increase) (Table 7). The total simulated vegetation C pool increases under eCO₂ compared to ambient CO₂ (0.41 kg C m⁻²), and thus the estimated plant P (estimated as a fraction of C:P ratios) increases as well (+0.45 g P m⁻²) (Fig. 6, Table 6). Thus, the simulated plant P demand is higher, and as the total available soil P for uptake is limited, the simulated excess C flux increases to 0.51 kg C m⁻² yr⁻¹. Moreover, despite the higher NPP under eCO₂ compared to simulated NPP under ambient CO₂, the BP is similar to the ambient CO₂ (2 % difference) due to the substantial increase in simulated excess C.

The simulated organic soil P values under eCO₂ were close to those under ambient CO₂ (1.6 g P m⁻²) (Table 7). This is due to the same parameterization of the output fluxes from this pool being used for eCO₂ and ambient CO₂. The simulated pool of inorganic P under eCO₂ decreases compared to the ambient CO₂ by 0.11 g P m-² due to the increased plant P pools and slight increase in uptake (+0.13 %).

However, the simulated sorbed organic and inorganic soil P values from eCO₂ are similar to those simulated under the ambient CO₂, which is due to the same parameterization of the sorption function (maximum sorption capacity) from the ambient CO₂ run, as explained in Sect. 3.1. Moreover, the modelled occluded and weathered soil P values were similar to those in the ambient CO₂ simulation (Table 7), which is due to the same prescribed observational data being used for this simulation.

4.2.2 Comparison of JULES-C, JULES-CN, and JULES-CNP under elevated CO₂

JULES-C and JULES-CN show higher vegetation and soil C pools, BP, and litter flux compared to JULES-CNP (Table 8, Fig. S2). Under eCO₂, simulated NPP using JULES-C and JULES-CN is 4.5 % higher than JULES-CNP and BP is 96.8 % higher in JULES-C and JULES-CN than in JULES-CNP, which simulates an excess C flux of 0.81 kg C m⁻² yr⁻¹, equivalent to 46 % P limitation under eCO₂. As a result of P limitation and eCO₂, the simulated CO₂ fertilization effect estimated based on changes in biomass under ambient and eCO₂ was reduced from 13 % in JULES-C and JULES-CN to 10 % in JULES-CNP. Moreover, the CUE from JULES-C and JULES-CN is 87.5 % higher than in JULES-CNP as a result of high P limitation over biomass production.

Table 8 Simulated C pools and fluxes with JULES-C and JULES-CN and the difference between them and JULES-CNP (in percent) under eCO₂. A positive percentage means larger respective values simulated with JULES C and JULES CN than with JULES-CNP (between period 2017–2018).

Inter-model comparisons under elevated CO₂

Following Fleischer et al. (2019), we report the simulated response to eCO₂ for the year 1999 (initial: CO₂ effect) and from 1999 to 2013 (15 years: final effect), which are both different from our evaluation period (2017–2018). Using JULES-C and JULES-CN under eCO₂, simulated GPP and NPP during the first year increase by 30 % and 61 % and by 28 % and 52 % after 15 years, respectively (Fig. 6). However, using JULES-CNP, eCO₂ increases simulated GPP, NPP, and BP responses during the first year by 29 %, 51 %, and 20 % and by 28 %, 43 %,and 7 % after 15 years, respectively.

Corresponding simulated CUE values during the first year and over 15 years show an increase of 24 % and 20 % in response to eCO₂ using JULES-C and JULES-CN, respectively. However, using JULES-CNP, simulated CUE values for the first year and after 15 years are reduced by 7 % and 17 %, respectively, in response to eCO₂.

Simulated total biomass (leaf, fine root, and wood C) (ΔC_veg) using JULES-C and JULES-CN for the first year and after 15 years of eCO₂ increased by 9 % and 13 %, respectively. However, using JULES-CNP ΔC_veg only increases this value by 0.5 % and 9 % for first year and after 15 years of eCO₂, respectively.

Figure 6 Relative effect of eCO₂ on simulated GPP, NPP, BP, CUE, ΔC_veg, leaf C, wood C, and fine root C using three versions of the JULES model over the first year (initial response) and after 15 years (final response).

4.3 Plant P demand, uptake, and excess C under ambient and elevated CO₂

To further understand the CP cycle dynamics, we studied the monthly averaged plant P demand and the relative (limited) P uptake (Eq. 26) under both ambient and elevated CO₂ conditions (Fig. 7).

Under ambient CO₂ conditions, the highest GPP is estimated at 0.29 ± 0.016 kg C m⁻² per month in July and the lowest value is 0.17 ± 0.051kg C m⁻² per month in October (Fig. 7a). The estimated WUE and SMCL in October is among the lowest estimated monthly values at 2.3 ± 0.51 kg CO₂ kg⁻¹ H₂O and 526.2 ± 31 kg m⁻², respectively (Fig. 7c). The highest P demand is estimated at 0.4 ± 0.02 g P m⁻² per month in July and the lowest demand is estimated at 0.2 ± 0.08 g P m⁻² per month in October. Consequently, the highest and lowest uptake are 0.32 ± 0.01 and 0.19 ± 0.07 g P m⁻² per month, respectively. The excess C for the highest and lowest GPP and demand periods are estimated at 0.4 ± 15 and 0.04 ± 0.07 kg C m⁻² per month, respectively.

However, similar to ambient CO₂, under eCO₂ conditions the highest estimated GPP is in July at 0.36 ± 0.017 kg C m⁻² per month and the lowest value is in October at 0.25 ± 0.062 kg C m⁻² per month (Fig. 7b). The estimated WUE and soil moisture content (SMCL) for the lowest GPP period are among the lowest monthly estimated values at 3.5 ± 0.74 kg CO₂ kg⁻¹ H₂O and 552 ± 33 kg m⁻² for October, respectively (Fig. 7d). The highest P demand is estimated for July at 0.51 ± 0.02 g P m⁻² per month, with an uptake flux of 0.31 ± 0.02 g P m⁻² per month, and the lowest P demand is estimated for October at 0.32 ± 0.1 g P m⁻² per month, with an estimated uptake flux of 0.26 ± 0.06 g P m⁻² per month. The highest excess C flux is also for July at 1.01 ± 0.17 kg C m⁻² per month and the lowest is for October at 0.27 ± 0.29 kg C m⁻² per month.

However, despite the P limitation in both eCO₂ and ambient CO₂ conditions, the P uptake flux under eCO₂ is higher than the ambient CO₂ condition. This is due to the higher WUE and increased SMCL (controlling uptake capacity, Eq. 27) under eCO₂ conditions, and thus more water availability during the dry season to maintain productivity and critically transport P to the plant (see Eq. 27), compared to ambient CO₂ condition (Fig. 7c and d). Additionally, in JULES both the vertical discretization (Burke et al., 2017) and mineralization terms (Wiltshire et al., 2021) depend on the soil moisture and temperature. Thus, higher P concentration and uptake are found under eCO₂ conditions.

Figure 7 Simulated monthly plant P demand and uptake (g P m⁻² per month), excess C, and GPP (kg C m⁻² per month) under (a) aCO₂ and (b) eCO₂ conditions. Water use efficiency (g m⁻² per month) under (c) ambient CO₂ (aCO₂) and (d) eCO₂ conditions. The grey area represents the standard deviation.

4.4 Soil P pool profiles under ambient CO₂ and elevated CO₂

We explored the distribution of the inorganic and organic soil P and their sorbed fraction within the soil layers and under different CO₂ conditions (Fig. S3). Both the ambient and eCO₂ simulations have a close inorganic soil P distribution at the topsoil layer (0–30 cm) (0.85 vs. 0.9 g P m⁻², respectively) and similar organic soil P distribution (0.85 vs. 0.9 g P m⁻², respectively).

However, the organic soil P and sorbed forms of inorganic and organic soil P profiles do not change significantly between different sets due to the similar parameterization of the processes that control these pools (processes that are related to the physical aspects of soils, and thus do not change under eCO₂ conditions) and the same parameter values being used for both ambient and eCO₂ runs.

Moreover, the soil P within 30 cm soil depth for ambient and eCO₂ conditions is at 14.7 and 14.56 g P m⁻², respectively, and the total ecosystem P for both ambient and eCO₂ conditions is at 35.97 g P m⁻². However, the slightly lower soil P in the eCO₂ conditions is due to the higher plant P demand compared to the ambient condition, which explains the higher allocated P vegetation (10 %) under eCO₂ conditions.

5 Discussion

Studies show the significant role of the tropical forests, and Amazonia in particular, in C uptake and regulating atmospheric CO₂ (Brienen et al., 2015; Phillips et al., 2017). As soil P availability is low in the majority of Amazonia (Quesada et al., 2012), the competition for nutrients by both plant and soil communities is high (Lloyd et al., 2001). The responses of these communities to eCO₂ under P-limited conditions remains uncertain (Fleischer et al., 2019). These responses in P-enabled models are represented in different ways regarding the excess C, which is not used for plant growth due to P limitation. Either growth is directly downregulated by taking the minimum labile plant C, N, and P (Goll et al., 2017), and photosynthesis is downregulated via V_cmax and J_max (Zhu et al., 2016; Yang et al., 2014; Comins and McMurtrie, 1993). Models like JULES-CNP downregulate NPP via respiration of excess carbon that cannot be used for growth due to plant nutrient constraints (Haverd et al., 2018). The estimated CUE depends on the modelling approach. Models that downregulate the photosynthetic capacity and GPP consequently (Zhu et al., 2016; Yang et al., 2014; Comins and McMurtrie, 1993) simulate a positive CUE response to CO₂ fertilization, while models that downregulate the NPP and respire the excess C (Haverd et al., 2018) simulate a negative CUE response (Fleischer et al., 2019), which is in line with field studies showing lower CUE when nutrient availability declines (Vicca et al., 2012b). However, this remains a major uncertainty in understanding the implication of P limitation on terrestrial biogeochemical cycles.

The JULES-CNP structure represents key P processes in both plant and soil pools and can be applied to the Amazon region using existing soil (Quesada et al., 2011) and foliar structural and nutrient (Fyllas et al., 2009) data for parameterization. The model can be applied globally and under future climate projections using global soil P data (Sun et al., 2021) for model initialization and PFT-specific plant (Zechmeister-Boltenstern et al., 2015) and soil stoichiometries (Zechmeister-Boltenstern et al., 2015; Tipping et al. (2016), sorption, and weathering ratios (based on lithological class-specific data from the GliM lithological map, Hartmann and Moosdorf, 2012, and soil shielding from Hartmann et al., 2014).

5.1 Evaluation of model performance

At the study site, JULES-CNP could reproduce the magnitude of soil organic and inorganic P pools and fluxes. The relative distribution of total organic P, total inorganic P and residue P fractions of total P in soils under Brazilian eucalyptus plantations (Costa et al., 2016) shows an inorganic P fraction of 28 % from total soil P, which is close to our estimation of 24 %, and organic P fraction of 30 % from total soil P, which is higher than our estimated fraction of 18 %. Thus, we may need to improve the process representation or the parameters that control the organic P concentration, such as litter flux and decomposition, soil organic P mineralization, and immobilization in the future.

Our estimated maximum P uptake, which represents the actual available P for plant uptake (Goll et al., 2017) for both ambient and eCO₂ is highly correlated with the plant P demand (R²=0.96 and 0.52, respectively). The plant P demand depends on the GPP changes, which are reflected by the WUE (Hatfield and Dold, 2019). Hence, under ambient CO₂, JULES-CNP simulates lower GPP and plant P demand during the dry season than during the wet season. Sufficient P uptake during these periods results in the lowest P limitation and thus the lowest simulated excess C. Nevertheless, under eCO₂ the same pattern is simulated, but a higher availability of soil P is simulated due to the stomatal closure in the dry season. Hence, due to the plant's more efficient water usage, the soil moisture in the dry season is higher (Xu et al., 2016), which impacts our capped P uptake flux (Eq. 27) and increases the uptake capacity.

Overall, JULES-CNP reproduced observed C pools and fluxes in acceptable ranges compared to the measurements. However, when using the JULES default V_cmax estimation method (Eq. 40), the model slightly underestimates the total GPP (2.9 vs. 3–3.5 kg C m⁻² yr⁻¹). Therefore, in this version of the model we used the improved V_cmax estimation method based on N and P (Eq. 46), which resulted in a final estimated GPP closer to the measurements (3.06 kg C m⁻² yr⁻¹).

Our results show an increase in GPP (21 %) in response to eCO₂ that is higher than the average increase of GPP reported in mature eucalyptus forests (11 %) that are also growing under low P soils at the free-air CO₂ enrichment experiment (EucFACE) facility in Australia (Jiang et al., 2020). This can be attributed to the lower decrease of biomass growth response estimated by JULES-CNP (−3 %) compared to the measurements from mature eucalyptus forests (−8 %) (Ellsworth et al., 2017) due to P limitation, which has been shown to impact the above-ground biomass growth response in mature forests (Körner et al., 2005; Ryan, 2013; Klein et al., 2016).

In order to estimate the biomass production (BP), we deducted the excess C fluxes from NPP. Using JULES-C and JULES-CN, the simulated biomass productivity enhancement due to eCO₂ (49 %) is in the middle range of the reported for the different biomes defined by Walker et al. (2021). Moreover, our estimated difference in BP between ambient and eCO₂ conditions (2 %) is close to the estimated difference for mature forests (3 %) (Jiang et al., 2020).

A global estimation for tropical forests using the CASA-CNP model, which includes N and P limitations on terrestrial C cycling, shows that NPP is reduced by 20 % on average due to the insufficient P availability (Wang et al., 2010), which is close to our estimated P limitation of 24 %. This finding is in line with a field study that shows a strong correlation between the total NPP and the soil-available P (Aragão et al., 2009).

The estimated decrease in NPP in response to eCO₂ as a result of P limitation is in line with the findings from the CLM-CNP model at five tropical forests (Yang et al., 2014), which indicates a dependency of CO₂ fertilization on the processes that affect P availability or uptake.

Our estimated CUE (0.31) is close to that of Jiang et al. (2020) for mature eucalyptus forests (0.31 ± 0.03) and to the measurements from our study site (0.34 ± 0.1). There is currently a lack of representation of stand age in JULES-CNP, a variable that can significantly affect CUE (e.g. mature trees are less responsive to the nutrient limitations) (De Lucia et al., 2007; Norby et al., 2016). However, the recent development of the Robust Ecosystem Demography (RED) model in JULES (Argles et al., 2020) and its integration into JULES-CNP in the future can address this issue.

Under low P availability, all available P is considered to be adsorbed or taken by plant and microbes for further consumption, with leaching considered to be minor within the timescales of our study period (Went and Stark, 1968; Bruijnzeel, 1991; Neff et al., 2000). Despite studies that show the possibility of P fixation as a source of available P for plants (Van Langenhove et al., 2020; Gross et al., 2021) due to the strong fixation of P in the soil (Aerts and Chapin, 2000; Goodale et al., 2002), the P deposited is unlikely to be available to plants in the short term (de Vries and Butterbach-Bahl, 2014); for this reason this version of JULES-CNP did not include P deposition. However, both P deposition and leaching are likely to have a very important role in sustaining the productivity of tropical forests in the Amazon over longer timescales (Van Langenhove et al., 2020) and need to be considered in future studies. Moreover, biochemical mineralization is also not included in the current version of JULES-CNP, which only accounts for total mineralization. However, models that include this process show no significant difference between total and biochemical mineralized P, which could be due to complexity of identifying the inclination of mineralization vs. uptake (Martins et al., 2021). Lastly, in order to capture the plant internal nutrient impact on the C storage, future work should focus on implementing recent developments, such as the non-structural carbohydrate (NSC) pools (Jones et al., 2020), in JULES-CNP.

5.1.1 Evaluation of model performance at test sites

Overall, inclusion of P processes in JULES-CNP improved the previously overestimated C fluxes and pools using JULES-C and JULES-CN. Generally, the biomass productivity tends to follow the observed P availability (Fig. 4), where the sites with higher available P for uptake simulated higher productivity, which is in line with observations of P availability in the Amazon (Aragão et al., 2009). Nevertheless, this tendency could be altered if the natural conditions in these forests are perturbated. For instance, in cases of high-mortality events in these P-limited sites (Malhi et al., 2009; Pyle et al., 2009), regrowing forests developing over the highly weathered Oxisols with limited available P (Davidson et al., 2004) results in the limitation shifting from P to N (Herbert et al., 2003). Hence, the controlling processes under N limitation will be N related, and thus processes such as N leaching or outgassing (Yang et al., 2014) will define the productivity. This shift in limitation conditions is not represented by JULES-CNP; therefore, at a few tested sites the model overestimated the P limitation and underestimated the productivity to the point that it was below the measured values. Moreover, the higher (than other sites) BP in JULES-C and JULES-CN at the Gigante Peninsula is related to the higher solar radiation in the forcing data at that site (Fig. S8).

The estimated litterfall and respiration fluxes were considerably lower with JULES-CNP than JULES-C and JULES-CN due to the lower simulated NPP with the former, which is in closer agreement with the observations at all sites. Consequently, the total vegetation and soil C pools have lower values under the P limitation (Malhi et al., 2009), which could not be captured by JULES-C and JULES-CN and successfully represented by JULES-CNP.

As shown in Fig. 5, JULES-CNP is highly sensitive to the five parameters needed to run JULES-CNP in addition to JULES-C and JULES-CN, and these were all prescribed for simulations at test sites. The successful model performance at these sites demonstrates the importance of these parameters in JULES-CNP, with implications for global-scale simulations.

5.2 Inter-model comparison of JULES-C, JULES-CN, and JULES-CNP

The comparison of simulated GPP enhancement across JULES versions for the first year is within the middle range of the first-year CO₂ responses of the C and CN models studied by Fleischer et al. (2019) evaluating simulated eCO₂ effects at a site in Manaus using the same meteorological forcing and methodology used in this study for a range of dynamic global vegetation models. However, comparison for 15 years of eCO₂ data shows that the simulated response with JULES-CNP is on the higher end of Fleischer et al. (2019) study, which is due to the higher estimated biomass growth by JULES-CNP (Table S1 in the Supplement). Similarly, when using JULES-CNP our estimated GPP enhancement is on the higher end of model estimations in Fleischer et al. (2019). Moreover, comparing the GPP responses between different versions of (JULES-C, JULES-CN, and JULES-CNP), JULES-CNP shows a slightly higher response to CO₂ fertilization associated with the higher WUE changes (Xiao et al., 2013) (Fig. S4). This is due to the higher sensitivity of the plant to water availability than P availability in the P-limited system (He and Dijkstra, 2014). Hence, under eCO₂, simulated transpiration is decreased (Sampaio et al., 2021) and photosynthesis is enhanced compared to ambient CO₂ due to the water-saving strategies of the plants and stomatal closure (Medlyn et al., 2016).

To that end, the monthly changes in WUE in JULES-CNP are highly correlated with the GPP, and thus the lowest and highest WUE values are found in the same periods as corresponding GPP values, which is similar to responses captured with the models studied by Fleischer et al. (2019) (Table S1).

Our estimated NPP enhancement using JULES-C and JULES-CN for both the first year and the 15-year period is within the middle range of the models in Fleischer et al. (2019). Nevertheless, the JULES-CNP response of BP is in the lower band of the CNP models in Fleischer et al. (2019) and close to the estimations from the CABLE (Haverd et al., 2018) and ORCHIDEE (Goll et al., 2017) models, which may be due to the similar representation of P processes and limitations between these models. However, our results show a 29 % decrease in NPP using JULES-CNP compared to JULES-C and JULES-CN, which is smaller than the differences between CLM-CNP and CLM-CN (51 % decrease) (Yang et al., 2014). The lower estimated decrease in JULES highlights the need for further study of the responses of corresponding plant C pools and fluxes to the changes in soil and plant P. Therefore, future work should be focused on the improvement of the total P availability and the plant C feedbacks. Moreover, there are other environmental factors (such as temperature) that show a possible impact on the CO₂ elevation and the changes of NPP (Baig et al., 2015) and need further examination in our model.

The CUE estimations of first-year and 15-year responses to CO₂ elevation from JULES-C and JULES-CN are in the middle range of C and CN models in Fleischer et al. (2019). However, the estimated CUE using JULES-CNP for first year and after 15 years are in the low range of CNP models reported by Fleischer et al. (2019), which is due to the same reason discussed for NPP comparison.

Finally, our estimated total biomass enhancement (ΔC_veg) using JULES-C and JULES-CN for the first-year and 15-year data are in the middle range of C and CN models from Fleischer et al. (2019) and are in lower range of CNP models from Fleischer et al. (2019) when using JULES-CNP. Nevertheless, while JULES-CNP includes the trait-based parameters (Harper et al., 2016), other functions such as flexible C allocation and spatial variation of biomass turnover are still missing, and future model improvement should be focused on their inclusion.

6 Conclusions

Land ecosystems are a significant sink of atmospheric CO₂ and thus buffer the anthropogenic increase of this flux. While tropical forests contribute substantially to the global land C sink, observational studies show a stalled increase in carbon gains over the last decade (Hubau et al., 2020; Brienen et al., 2015). However, modelling studies that lack representation of P cycling processes predict an increasing sink (Fleischer et al., 2019; Fernández-Martínez et al., 2019). This is particularly relevant for efforts to mitigate dangerous climate change and assumptions about the future efficacy of the land C sink. Therefore, in this study we presented the full terrestrial P cycling and its feedback on the C cycle within the JULES framework. Our results show that the model is capable of representing plant and soil P pools and fluxes at a site in the central Amazon and across the extended P-limited test sites in the Amazon, Gigante Peninsula (Panama), and Hawaii chronosequences provided using site-level data for model parameterization. Moreover, the model estimated a significant NPP limitation under ambient CO₂ due to the high P deficiency at these sites, which is representative of the central Amazon and tropical P-limited sites, and elevated CO₂ resulted in a further subsequent decrease in the land C sink capacity relative to the model without P limitation. While our study is a step toward the full nutrient cycling representation in ESMs, it can also help the empirical community to test different hypotheses (i.e. dynamic allocation and stoichiometry) and generate targeted experimental measurements (Medlyn et al., 2015).