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

21 Most Land Surface Models (LSMs), the land components of Earth system models (ESMs), include 22 representation of nitrogen (N) limitation on ecosystem productivity. However only few of these models have 23 incorporated phosphorus (P) cycling. In tropical ecosystems, this is likely to be important as N tends to be 24 abundant but the availability of rock-derived elements, such as P, can be very low. Thus, without a 25 representation of P cycling, tropical forest response in areas such as Amazonia to rising atmospheric CO 2 26 conditions remains highly uncertain. In this study, we introduced P dynamics and its interactions with the N and 27 carbon (C) cycles into the Joint UK Land Environment Simulator (JULES). The new model (JULES-CNP) 28 includes the representation of P stocks in vegetation and soil pools, as well as key processes controlling fluxes between these pools. We 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 parameterise, calibrate and evaluate JULES-CNP. Novel soil and plant P pool observations are used for parameterisation and calibration and the model is evaluated against C fluxes and stocks, and for those soil P pools not used for parameterisation/calibration. We then apply the model under elevated CO 2 (600 ppm) at our study site to quantify the impact of P limitation on CO 2 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 2 . We estimate P to limit net primary productivity (NPP) by 24% under current CO 2 and by 46% under elevated CO 2 . Under elevated CO 2 , biomass in simulations accounting for CNP increase by 10% relative to at 39 contemporary CO 2 , although it is 5% lower compared with CN and C-only simulations. Our results highlight the potential for high P limitation and therefore lower CO 2 fertilization capacity in the Amazon forest with low fertility soils.


45
Land ecosystems currently take up about 30% of anthropogenic CO2 emissions (Friedlingstein et al., 2020), thus buffering the anthropogenic increase in atmospheric CO2. 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

118
Here, we describe the development and implementation of the terrestrial P cycle in the Joint UK Land   Description of the plant P pool (Pp) follows Zhu et al., (2016) and is estimated as the difference between the 193 input, plant uptake ! "# (eq.26) and output of this pool, plant litter flux ! "#$ (eq.28), with both fluxes 194 expressed in kg P m -2 yr -1 as follows: 195 196 43 * 40 = 3 5* − 3 "$0 (eq.5)

256
Therefore, actual NPP is for our purposes equal to Biomass Production (BP), and is calculated as potential NPP 257 minus excess C (lost to the plant through autotrophic respiration), with the latter the C that cannot be used to 258 growth new plant tissue due to insufficient plant nutrient supply. Hence, if the system is limited by the 259 availability of N and/or P, NPP will be adjusted to match the growth that can be supported with the limited N or 260 P supply, with any excess carbohydrate lost through excess C.

261
The total excess C term (y 0 ) (kg C m -2 yr -1 ) is calculated as:

265
where y ( and y # are the excess C fluxes due to growth (g) and spread (s) and are assumed to be rapidly respired 266 by plants.

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

272
The litter production in JULES before limitation is estimated based on the as follows: 288 289 is the inverse of whole plant C:P ratio, is inverse plant C:N ratio,

299
(see Clark et al., (2011) for more detail), Π , is nutrient-unlimited, or potential, NPP (kg C m -2 yr -1 ), y ( is excess 300 C due to either P or N limitation for plant growth (kg C m -2 yr -1 ) and y # is excess C due to either P or N 301 limitation for vegetation spreading (kg C m -2 yr -1 ).

303
Equations 20 and 22 are solved by first setting y ( = 0.0 to find the total plant P (eq. 20) and N demand (eq.22).

304
If the P and N demand for growth are less than the available P and N and fractional coverage ( ) (NPP fraction

311
Similarly, in order to estimate the P and N demand for spreading (eq. 21 and 23), initially the excess C from

318
Plant P uptake ( * 5* ) (arrow a in Fig 1) is estimated based on the P demand for growth and spreading (∅ 0 ) and 319 the root uptake capacity ( +.H ) (kg P kg -1 C yr -1 ), as follows:

323
Description of the plant P uptake ( * 5* ) varies spatially depending on the root uptake capacity ( +.H ) followed 324 by Goll et al., (2017). Therefore, in regions with limited P supply, the plant P uptake is limited to the +.H and 325 consequently impacts the excess C and BP.

333
Description of the litter production of P ( 3 % "$0 ) (arrow b in Fig 1) follows JULES-CN as in Wiltshire et al.,

334
(2021) and is calculated based on the litter flux of C (kg C m -2 yr -1 ) using leaf, root and wood turnovers (yr -1 ),

335
and through the vegetation dynamics due to large-scale disturbance and litter production density, as follows:  408 409

484
where $%0 is the estimated intercept and #" is the slope of the linear regression derived for the Vcmax estimation.

509
Where " is specific leaf density (kg C m -2 per unit LAI), : is balanced (or seasonal maximum) leaf area index 510 (m 2 m -2 ), B" is allometric coefficient (kg C m -2 ) and B" is allometric exponent.

527
The organic and inorganic soil P assumed to be always at equilibrium with the relative sorbed pools (Wang, equilibrium state of sorbed and free-soil P. Moreover, as the magnitude of changes in the occluded and parent 532 material pools are insignificant over a short-term (20 years) simulation period (Vitousek et al., 1997), these two 533 pools were prescribed using observations. Remaining parameters used to describe soil P fluxes (eq.s 27-44)

534
were prescribed using values from the literature (Table 3).

536
We used a combination of data from Study site and the nearby site K34 for model evaluation of C fluxes (GPP,

537
NPP) and C pools (soil and vegetation C, leaf, root and stem C) with no calibration on plant and soil organic and 538 soil inorganic P polls included (Table 3).

558
We evaluate the impact of including a P cycle in JULES using three model configurations (JULES C, CN and 559 CNP). We apply JULES in all three configurations using present day climate under both ambient CO2 and

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

593
To test the sensitivity of the P and C related processes to the model P parameters, six sets of simulations were

631
JULES CNP-CNP could reproduce the plant and soil C ( Figure.2 and Table 5) and N pools and fluxes ( Figure   632 S6 and Table 6) pools and fluxes under ambient CO2. Our results show that simulated GPP, is within the range 633 of measurement (3.02 kg C m -2 yr -1 model vs 3-3.5 kg C m -2 yr -1 observed, respectively, Table 5).

635
Simulated NPP, is close to the measured values (NPP: 1.14 -1.31 observed vs 1.26 modelled kg C m -2 yr -1 ) with 636 autotropic respiration (RESP) also closely following the observations (1.98 observed vs 1.81 modelled kg C m -2 637 yr -1 ). Biomass production is estimated as a difference between NPP and the amount of C which is not fixed by 638 plants due to the insufficient P in the system (excess C) (eq. 27). The excess C flux is highly dependent on the 639 plant P and the overall P availability to satisfy demand (Table 5). Simulated flux of excess C is 0.3 kg C m -2 yr -1 640 under ambient CO2. In JULES-CNP this flux is subtracted from NPP in order to give the BP (eq. 17) ( Table 5).

670
The results indicate that among all the corresponding C and P pools and fluxes, the excess C flux -which 671 demonstrates P limitation to growth -shows the highest sensitivity to changes in C:P ratios, KP and 672 K ]^)bcd , K`a )bcd . A decrease in plant C:P results in a large increase in excess C. This is due to the higher plant 673 P demand as a result of lower plant C:P ratios. An increase in the uptake factor and maximum sorbed organic 674 and inorganic P also results in an increase in excess C. This is due to the higher uptake demand through higher 675 uptake capacity (due to higher KP) and lower available P for uptake due to higher organic and inorganic sorbed 676 P (due to higher K ]^)bcd , K`a )bcd ). Since the total P in the system is lower than the plant demand, the uptake 677 capacity and sorbed P, higher P limitation is placed on growth (decreasing BP) which results in an increase in 678 excess C and decrease in plant C, but also soil C which is a result of lower litter input (Figure 4). Total soil P

679
shows low sensitivity to changes in plant C:P and uptake factor but high sensitivity to maximum inorganic 680 sorbed P. Moreover, sorbed P shows middle to high sensitivity to maximum organic and inorganic sorbed P 681 respectively ( Figure. S5). Nevertheless, organic and inorganic P adsorption coefficients (K \]^_)]^, K \]^_)`a )

682
show no sensitivity to C and P pools and fluxes. This is due to limiting the organic and inorganic P sorption 683 terms controlled only by maximum sorption, hence no effect applied by organic and inorganic adsorption

694
The eCO2 simulation using JULES CNP yields a higher GPP compared to the ambient CO2 (0.83 (kg C m -2 yr -1 ) 695 increase), as a result of CO2 fertilization. Moreover, due to the GPP increase, NPP and RESP follows the same 696 trend and increased compared to ambient CO2 (NPP: 0.49 and RESP:0.3 (kg C m -2 yr -1 ) increase) ( Table 5). The

697
total simulated vegetation C pool increases under eCO2 compared to ambient CO2 (0.41 kg C m -2 ), hence the 698 estimated plant P (estimated as a fraction of C:P ratios) increases as well (+0.45 (g P m -2 )) ( Fig 6, Table 4).

699
Thus, the simulated plant P demand is higher, and as the total available soil P for uptake is limited, the simulated 700 excess C flux increases to 0.51(kg C m -2 yr -1 The simulated organic soil P under eCO2 yields close to the ambient CO2 (1.6 g P m -2 ) ( Table 5). This is due to 708 the same parameterization of the output fluxes from this pool for eCO2 and ambient CO2. The simulated pool of ambient CO2 which is due to the same parameterizing of sorption function (maximum sorption capacity) from 713 the ambient CO2 run as explained in calibration section. Moreover, the modelled occluded and weathered soil P 714 yield similar to those in the ambient CO2 simulation (Table 5) which is due to the same prescribed observational 715 data that was used for this simulation.

740
Corresponding simulated CUE during the 1 st year and 15 years shows an increase of 24% and 20% in response 741 to eCO2 using JULES C/CN respectively. However, using JULES CNP, simulated CUE for the 1 st and after 15 742 years is reduced by 7% and17% in response to eCO2.

744
Simulated total biomass (leaf, fine root and wood C) (DCveg) using JULES C/CN for the 1 st and 15 years of

756
To understand further the CP-cycle dynamics, we studied the monthly averaged plant P demand and the relative 757 (limited) P uptake (eq. 26) under both ambient and elevated CO2 conditions ( Figure. 6).

759
Under ambient CO2 condition the highest GPP is estimated at 3.5±0.19 kg C m -2 month -1 in July and the lowest 760 at 2.06±0.61kg C m -2 month -1 in October ( Figure. 6-a). The estimated WUE and SMCL in October is among the 761 lowest estimated monthly values at 2.3±0.51 kg CO2/kg H2O and 526.2±31 kg m -2 respectively ( Figure. 6-c).

762
The highest P demand is estimated at 0.4±0.02 g P m -2 month -1 in July and the lowest demand at 0.2±0.08 g P m -763 2 month -1 in October. Consequently, the highest and lowest uptake (0.32±0.01 and 0.19±0.07 g P m -2 month -1 , 764 respectively). The excess C for the highest and lowest GPP and demand periods are estimated at 0.4±15 and 765 0.04±0.07 kg C m -2 month -1 , respectively.

767
However, similar to ambient CO2, under eCO2 condition the highest estimated GPP is in July at 4.36±0.21 kg C 768 m -2 month -1 and lowest for October 3.02±0.75 kg C m -2 month -1 ( Figure. 6-b). The estimated WUE and soil 769 moisture content (SMCL) for the lowest GPP period is among the lowest monthly estimated values at 3.5±0.74 770 kg CO2/kg H2O and 552±33 kg m -2 for October respectively ( Figure. 6-d). The highest P demand is estimated 771 for July at 0.51±0.02 g P m -2 month -1 with the uptake flux of 0.31±0.02 g P m -2 month -1 and the lowest demand 772 is estimated for October at 0.32±0.1 g P m -2 month -1 with the estimated uptake flux of 0.26±0.06 g P m -2 month -773 1 . The highest excess C flux is also for July at 1.01±0.17 kg C m -2 month -1 and lowest for October 0.27±0.29 kg 774 C m -2 month -1 , respectively.

776
However, despite the P limitation in both eCO2 and ambient CO2 conditions, the P uptake flux under eCO2 is 777 higher than the ambient CO2 condition. This is due to the higher WUE and increased SMCL (controlling uptake 778 capacity (eq. 27)) under eCO2 condition, hence more water availability during the dry season to maintain 779 productivity and critically transport P to the plant (see eq. 27), compared to ambient CO2 condition ( Figure.  soil P distribution at the topsoil layer (0-30cm) (0.85 vs. 0.9 (g P m -2 ) respectively) as well as similar organic 795 soil P distribution (0.85 vs 0.9 (g P m -2 ) respectively).

797
However, the organic soil P and sorbed forms of inorganic and organic soil P profiles are not changing 798 significantly between different sets due to the similar parameterization of the processes that control these pools (processes which are related to the physical aspects of soils, hence not changing under eCO2 condition) and the

823
However, this remains a major uncertainty in understanding the implication of P limitation on terrestrial 824 biogeochemical cycles.

825
Our new developments include major P processes in both plant and soil pools and can be applied to the Amazon

837
Eucalyptus plantations (Costa et al., 2016) shows inorganic P fraction of 28% from total soil P which is close to 838 our estimation of 24% and organic P fraction of 30% from total soil P which is higher than our estimated 839 fraction of 18%. Thus, we may need to improve the process representation or parameters that control the organic 840 P concentration, such as litter flux and decomposition, soil organic P mineralization, and immobilization in the 841 future.

843
Our estimated maximum P uptake, which represents the actual available P for plant uptake (Goll et

853
Overall, JULES-CNP reproduced the observed C pools and fluxes which are in the acceptable ranges compared 854 to the measurements. However, using the JULES default Vcmax estimation method (eq. 40), the model slightly 855 underestimates the total GPP (2.9 kg C m -2 yr -1 vs. 3-3.5 kg C m -2 yr -1 ). Therefore, in this version of the model,

856
we used the improved Vcmax estimation method based on N and P (eq. 46) which resulted a final estimated GPP 857 closer to the measurements (3.06 kg C m -2 yr -1 ).

859
Our results show an increase in GPP (21%) in response to eCO2 which is higher than the average increase of 860 GPP reported in mature eucalyptus forests (11%), also growing under low P soils at the free air CO2 enrichment

866
In order to estimate the biomass production (BP), we deducted the excess C fluxes from the NPP. Using JULES 867 C/CN models our estimated biomass productivity enhancement due to eCO2 (49%) is in the middle range of the 868 reported various studies from different biomes by Walker et al., (2021). Moreover, our estimated difference of 869 BP between ambient and eCO2 conditions (2%) is close to the estimated difference for mature forests (3%) 870 (Jiang et al., 2020).

871
A global estimation for tropical forests using CASACNP model which includes N and P limitations on 872 terrestrial C cycling, shows that NPP is reduced by 20% on average due to the insufficient P availability (Wang,