Diagnosing Phosphorus Limitation in Subtropical Forests in China under Climate Warming

Abstract

Phosphorus (P) is a fundamental component of plant organisms. Most of the Earth System Models (ESMs) project increases in future Net Primary Productivity (NPP) due to climate warming. However, large uncertainty exists in projected NPP due to future P limitation. Subtropical China is a region with high vegetation NPP, but its forests are mostly P limited. In this study, we used the simulations of Atmospheric-Vegetation Interaction Model 2 (AVIM2) to diagnose the P limitation in forests in this region, and found that climate warming in the period of 1951–2010 had enhanced P limitation. The P demand during 1981–2010 for Evergreen Broad-leaved Forest (EBF) and Evergreen Needle-leaved Forest (ENF) are 1.67 and 1.8 times than that during 1951–1980, respectively. The observed current Available Soil P (ASP) density in 4 representative forest sites in subtropical China varied between 940 $\mathrm{mg}·{\mathrm{m}}^{-2}$ and 2365 $\mathrm{mg}·{\mathrm{m}}^{-2}$, and the P demands account for 0.86% to 25.5% of the ASP for the period of 1951–2010. Future P demands are estimated to account for 3.2% to 68.3% of the current ASP at the end of this century for RCP8.5. Therefore, forests, especially plantations, in subtropical China are facing high risks of P limitation.

1. Introduction

Phosphorus (P) is a key component of biomacromolecules in all organisms, and it plays a crucial role in terrestrial ecosystems processes [1]. Theoretical and experimental studies have demonstrated that the terrestrial carbon (C) cycle is strongly affected by P availability [2,3], and the stoichiometry of C, nitrogen (N) and P in plant and soil affects assimilation and allocation of C in plant and ultimately affects the functioning of terrestrial ecosystems and the absorption of CO₂ [4,5]. P is derived primarily from rock weathering, and it is not clear if there is additional external P entering the ecosystem due to climate warming. However, the increase of C in ecosystem due to climate warming is certain, which implies that the ecosystem will consume more P to maintain the consistent C:P ratio. Therefore, climate warming may result in the P limitation in ecosystem which will finally influence its C uptakes.

Climate warming has profoundly affected the terrestrial ecosystem and global biogeochemical cycle, as the average global surface temperature has increased almost 1 °C since preindustrial times [6]. Previous studies have shown that climate warming enhances photosynthesis, and thus promotes vegetation growth [7,8]. Net primary productivity (NPP) is the difference of gross photosynthesis minus autotrophic respiration, which represents the function of ecosystem and is considered as a critical indicator for vegetation growth. Large-scale model simulations have shown that both global [9,10,11] and China mean NPP have increased significantly in recent decades [12,13,14,15,16]. The increase of NPP implies that vegetation takes up and stores more C in the ecosystem, which indicates more P consumption in the ecosystem. However, ecosystems begin their existence with a fixed complement of P from which even very small losses cannot readily be replenished [17].

Globally, N and P are the most common nutrients limiting plant growth and soil C storage [18]. Previous studies show that some tropical forests and savanna ecosystems are P limited because of their old soil age [19,20]. Mahowald et al. (2008) predicted that as humans continue to preferentially increase N deposition, ecosystems may shift from N limitation to P limitation [21]. Subtropical China is the region with the highest vegetation gross primary productivity (GPP), and its regional total GPP contributes the largest proportions in China’s total GPP [22]. Previous experiments have shown that some of the subtropical forests in China are P limited, and some are both P and N limited [23,24,25,26,27]. With the intensification of climate warming, it is expected that P will be a crucial nutrient factor affecting the C uptake and storage in subtropical China. However, the spatial and temporal distribution of current P consumptions and future P limitation in subtropical China are still unclear. Therefore, it is of great importance to explore the dynamic changes of P demand in subtropical China, which will be of great help for understanding the potential threats of P deficit in subtropical China and to make effective strategy to cope with climate change.

Modeling and diagnosing are two major approaches to studying the P limitation in ecosystems. In terms of the modelling approach, some process-based models have been developed in the last 20 years for explicitly modeling P cycle. For example, Wang et al. (2010) developed a P model to simulate the global distribution of P limitation on NPP at steady state, and they found that tropical vegetation NPP was reduced by about 20% on average by P limitation [28]. Goll et al. (2012) incorporated a P cycle into a land surface model and found that the average plant uptake of P over 1970–1999 in tropical and temperate forests was in the range of 0.18–0.71 $\mathrm{g}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ [29]. Yang et al. (2014) introduced P dynamics and C-N-P interactions into the community land model version 4 with a prognostic C and N model (CLM4-CN) to simulate the effects of elevated CO₂ concentration on plant P availability and got meaningful results, and they found that tropical forest response to elevated CO₂ can only be predicted if the interactions between C cycle and nutrient dynamics are well understood and represented in models [30]. Therefore, more scientific experiments are needed for understanding the C-P interactions, and the modeling approach for understanding P limitation still remains highly uncertain.

In terms of the diagnostic approach, Peñuelas et al. (2013) used globally uniform stoichiometric ratios for vegetation and soil and calculated N and P demands for given changes in C pools during 2000–2099, and found that global P demand are ranged between −0.9 and 4.3 $\mathrm{Pg}\mathrm{P}$ for the lower projection for Representative Concentration Pathway (RCP) 2.6 and between −1.7 and 6.5$\mathrm{Pg}\mathrm{P}$ for RCP8.5, with negative values indicating a P surplus [31]. However, the size of the global plant-accessible labile P soil pool is only 13.2$\mathrm{Pg}\mathrm{P}$, and, hence, the limited P availability is likely to reduce future carbon storage by natural ecosystems during this century. Wieder, et al. (2015) used the NPP outputs from the Coupled Climate Carbon Cycle Model Intercomparing Projects (CMIP5) models over 1860–2100 and estimated annual new N and P inputs to evaluate how CO₂ fertilization effects could be constrained by nutrient availability [32]. They found that accounting for N and P limitation would lower projected end-of-century estimates of NPP by 19%. The approaches of Peñuelas et al. (2013) [31] and Wieder et al. (2015) [32] were based on C stock and NPP changes, respectively. Sun et al. (2017) used both C stock and NPP-based approaches in a consistent framework to investigate whether available soil P (ASP) can support the P demand with four scenarios of soil P availability, and they found that there are large differences between the NPP-based demand and C stock-based demand, and the diagnostic approach still remains uncertain [33].

Previous diagnostic approaches showed contradictory results in prediction of P limitation in subtropical China. For example, Sun et al. (2017) [33] estimated the values of additional P demand under RCP8.5 in subtropical China were negative based on C stock method, but positive based on NPP method. Therefore, the estimation of P limitation in subtropical China still remains highly uncertain.

The uncertainty in estimates of P limitation mainly comes from three sources: The uncertainty of the C:P ratios for plant organs, the NPP allocation to leaves, roots and wood and the estimates of the ASP. To refine the P limitation study in subtropical China, we obtained biome-specific C:P ratios for plant organs, and the values of NPP allocation to leaves, roots and wood from literatures which reported field measurements in subtropical China. Furthermore, we conducted field investigations in four representative forest ecosystems to obtain the current ASP density. The main objectives of this study are: (1) To diagnose the past P consumption in subtropical forests in China based on field observations; (2) to investigate future P demand under climate warming; and (3) to understand factors controlling the spatial distribution of P demand in subtropical China. The ultimate objective is to find the potential risk of P deficit in subtropical China and to provide scientific suggestions for forest managers to cope with climate change. To achieve these objectives, we used NPP outputs from an atmosphere–vegetation interaction model (AVIM2) in combination with observed stoichiometric relationships to diagnose the past P demand and with the field-observed ASP to diagnose the P limitation under RCP8.5. This approach can improve our understanding of the P limitation in China’s subtropical forest, and the results may provide useful information for forest managers to cope with climate change.

2. Materials and Methods

2.1. Study Area

The subtropical region in mainland China, which consists of 16 provincial administrative regions, with a total area of 2.26 million square kilometers, was taken as the study area (Figure 1). The climate in the study area is warm and humid and suitable for vegetation growth. Forests dominate 42% of the study area, including 13.29% evergreen broad-leaved forest (EBF), 27.39% evergreen needle-leaved forest (ENF) and 1.38% deciduous broad-leaved forest (DBF).

2.2. Data Sources

The meteorological data for driving the AVMI2 model is from the Global Meteorological Forcing Dataset (GMFD) of Princeton University, with a spatial resolution of 1° × 1° and a time interval of 3 h [34]. The vegetation map used in this model is the global land cover dataset (GLC2000), which is available at http://www-gem.jrc.it/glc2000/objectivesGLC2000.htm, and the soil texture data is a combination of three regional maps of similar scales (1:14,000,000, 1:4,000,000 and 1:5,000,000) [35]. All of the model driven datasets are interpolated to the spatial resolution of 0.1° × 0.1°. The future temperature and NPP data are the simulations of BCC-CSM 1.1, which was obtained from CMIP5 multi-model dataset through the Program for Climate Model Diagnosis and Intercomparison (PCMDI) at http://pcmdi9.llnl.gov/esgf-web-fe/.

2.3. Methods

2.3.1. Atmosphere–Vegetation Interaction Model (AVIM2)

The AVIM2 used in this study is a process-based model which was first developed by Ji (1995) [36], and it is the land surface component of the BCC-CSM 1.1. BCC-CSM 1.1 is one to the ESMs, which participate in CMIP5. AVIM2 consists of three modules: Plant-growth module, soil–vegetation–atmosphere transfer module and soil C and N dynamics module. The detailed description of the modules in AVIM2 can be found in Ji (1995) [36] and Huang et al. (2007) [37]. The AVIM2 simulations are widely used in studying the responses of terrestrial ecosystem C cycle to climate change in China [16,38]. In AVIM2, GPP is the aggregate amount of leaf photosynthesis at every model time step, and NPP is the difference between GPP and a similarly aggregated amount of autotrophic respiration (Ra, as the sum of maintenance and growth respiration).

2.3.2. Diagnose of Additional P Demand

The method used to diagnose the additional P demand in forest ecosystems based on NPP changes is mainly in reference to that of Sun et al. (2017) [33]. Here, additional means the demand due to the increases of C storage since 1950. We assume that 100% of the plant uptake of P by NPP is equal to the P supply from all sources at 1950; that is, the additional P demand was zero at 1950. Then, the additional P demand ($\Delta {\mathrm{P}}_{i,\left(t\right)}$) in ecosystems related to NPP change is calculated as the total P demand minus changes in recycled P supply from leaf P resorption ($\Delta P_{RSB,i\left(t\right)}$) and changes in the NPP uptake of mineralized P ($\Delta P_{MR,i\left(t\right)}$):

$$\Delta {\mathrm{P}}_{i,\left(t\right)}={\displaystyle \sum }_{j=1}^3\frac{\Delta NP{P}_{i,j,\left(t\right)}\times d_{i,j}}{r_{i,j}}-\Delta P_{RSB,i\left(t\right)}-\Delta P_{MR,i\left(t\right)}$$

(1)

where i and t represent forest type and year, respectively, and the numbers 1, 2 and 3 indexed by j represent the C pool for root, wood and leaf, respectively. The $r_{i,j}$ is the C:P ratios in the corresponding C pool, and $d_{i,j}$ is the fraction of NPP allocated to the C pool of leaf, wood and root, respectively.

${\Delta \mathit{NPP}}_{i,j,\left(t\right)}$ is the difference of NPP between two years:

$${\Delta \mathit{NPP}}_{i,j,\left(t\right)}=NP{P}_{i,j,\left(t\right)}-NP{P}_{i,j,\left(t-1\right)}$$

(2)

The changes of leaf P resorption available for the increment of NPP of the next year was assumed to be only associated with the leaf NPP change, which is given by:

$$\Delta P_{RSB,i,\left(t\right)}=\frac{{{\displaystyle \sum }}_{j=1}^3\Delta NP{P}_{i,j,\left(t\right)}\times d_{i,3}\times RS{B}_i}{r_{i,3}}$$

(3)

where $RS{B}_i$ is the adsorption coefficient.

The mineralization of P from SOM and litter is another source of inorganic P for NPP uptake, which is the difference between microbial immobilization and gross soil organic matter and litter mineralization. This flux ($\Delta P_{MR,i\left(t\right)}$) was inferred using two idealized scenarios. Scenario 1 (${\mathrm{S}}_1$) assumes low mineralization, that is, only P in fallen leaves can be mineralized and support NPP in the following year, while P in wood and root components of litter is unavailable for vegetation to recycle. This is given by:

$${\Delta P}_{\mathit{MR},i,S_1\left(t\right)}=\frac{\Delta NP{P}_{i,\left(t\right)}\times d_{i,3}\times \left(1-RS{B}_i\right)}{r_{i,3}}$$

(4)

Scenario 2 (${\mathrm{S}}_2$) assumes high mineralization, that is, 90% of P in the wood and root litter is recycled in addition to all P in fallen leaves each year. This is given by:

$${\Delta P}_{\mathit{MR},i,S_{2}t}=\frac{\Delta NP{P}_{i,\left(t\right)}\times d_{i,3}\times \left(1-RS{B}_i\right)}{r_{i,3}}+fr\mathrm{a}c\times \Delta NP{P}_{i,\left(t\right)}\times \left(\frac{d_{i,2}}{r_{i,2}}+\frac{d_{i,1}}{r_{i,1}}\right)$$

(5)

where frac is the fraction of P in litter reused in the ecosystems, which is set as 0.9.

The accumulated additional P demand (AAPD) for a given period is given by:

$$AAPD={\displaystyle \sum }_{t1}^{t2}\Delta {\mathrm{P}}_{i,\left(t\right)}$$

(6)

where $t_1$ and $t_2$ represent the start and end year, respectively.

The parameters used for diagnosing P limitation are listed in

Table 1. Parameters and sources for diagnosing phosphorus limitation in subtropical forests in China.

Pool		${\mathit{r}}_{\mathit{i},\mathit{j}}$	Ref	${\mathit{d}}_{\mathit{i},\mathit{j}}$	Ref	RSB	Ref
EBF	root	433.62 (177.02–585.82)	[39,40]	0.22	[22]	0.58	[41]
	wood	479.62 (269.02–828.16)	[39,40]	0.72
	leaf	431.41 (208.11–779.13)	[39,42,43,44,45,46]	0.06
ENF	root	386.92 (224.84–404.80)	[39,40]	0.19		0.58
	wood	726.67 (554.01-1025.48)	[39,40]	0.73
	leaf	568.55 (333.12–843.06)	[39,42,43,44,45,46]	0.08
DBF	root	513 (410–615)	[39,40]	0.24		0.46
	wood	2187 (1750–2625)	[39]	0.70
	leaf	622 (192–1052)	[39,42,43,44,45,46]	0.06

EBF, ENF and DBF are the abbreviations for evergreen broad-leaved forest, evergreen needle-leaved forest and deciduous broad-leaved forest, respectively.

, in which $r_{i,j}$ is obtained from literatures that reported the field measurements in subtropical forests in China. We used the fractions of C pools (leaves, wood and roots) to total vegetation C to represent NPP allocated to leaves, wood and roots, and $d_{i,j}$ were obtained by averaging over 200 field observations in subtropical forests [22].

2.3.3. Theil-Sen Regression

The Theil-Sen regression analysis is used to study the trends of NPP and P demand in this study. Theil-Sen regression is a method for robust linear regression that chooses the median slope among all lines through pairs of two-dimensional sample points. It has been called the most popular nonparametric technique for estimating a linear trend and is widely used in geoscience and remote sensing studies [47,48].

For a set of two-dimensional points ($x_i$, $y_i$), the slope ($Se{n}_{i,j}$) is estimated as follows:

$$Se{n}_{i,j}=\mathrm{Median}[\frac{y_i-y_j}{x_i-x_j}:x_i\ne x_j,i<j=1,2,\dots ,n]$$

(7)

The Mann-Kendall trend test is used to test the significance of the trends [49,50].

2.3.4. Field Sampling

The field investigations were conducted in July 2018 in Shennongjia, Qian Yanzhou, Mt. Gongga and Dinghushang. The locations of the forest sites, including Shennongjia (109°59′–ll0°58′ E, 3l°15′–3l°57′ N), Mt. Gongga (29°20′–30°20′ N, 101°30′–102°15′ E), Dinghushang (23°09′–23°11′ N, 112°31′–112°34′ E) and Qian Yanzhou (26°12′–26°34′ N, 115°10′–115°23′ E), are shown in Figure 1.

Shennongjia is located in the broad-leaved forest ecological zone of Qinba Mountain, which is a key area for biodiversity conservation in China. The forest in Shennongjia represents a natural mountainous forest in the northern subtropical China. The annual mean temperature is 10.6 °C, and the annual total precipitation is 1219 mm in Shennongjia [51]. The soil samples were obtained at the altitudes of 760 m for EBF and 2080 m for ENF. The soils of both EBF and ENF are yellow-brown clay.

Mt. Gongga is located in the transitional area between China’s eastern monsoon subtropics and the frigid area of the Tibetan Plateau. The forest in Mt. Gongga represents a natural mountainous forest in western subtropical China. The annual mean temperature is 4.4 °C and the annual total precipitation is 1910 mm in Mt. Gongga. The soil samples were obtained at the altitudes of 2265 m for EBF and 3082 m for ENF. The soil for EBF is brown silt soil and that for ENF is dark-brown silt soil.

Dinghushan is located in southern tropical China, and its forest represents a southern tropical natural forest in China. The annual total precipitation is 1680 mm and the annual mean temperature is 22.3 °C in Dinghushan [52]. The sampling plots were located at the altitudes of 89 m for EBF and 210 m for ENF. The soils for both EBF and ENF are red loam.

Qian Yanzhou is located at the central subtropical region in China, with an annual total precipitation of 1485 mm and an annual mean temperature of 17.9 °C. The plantation in Qian Yanzhou represents a man-made forest in the central subtropical region. The soil samples were obtained at the altitudes of 55 m for EBF and 103 m for ENF. The soils are red sandy loam soils for both EBF and ENF.

Four soil quadrats (10 × 10 m each) were set under each EBF and ENF in each forest site. Five soil samples were randomly taken at the top 50 cm of soils in each quadrat, and they were evenly mixed together. For each forest type in each site, four soil samples were taken back for analysis, and all soil samples were air-dried and had small stones and any visible plant and root residues removed, and were then ground and passed through a 2 mm sieve for analysis.

2.3.5. Laboratory Analysis

The ASP is the main source of P that plants can absorb and utilize, which represents the actual P supply capacity in soil. Different approaches for calculating ASP have resulted in a large discrepancy in the estimates, which is the major source of uncertainty in P limitation study [28,29,30,33]. To reduce this uncertainty, we conducted field observations in representative forest sites in the subtropical region to estimate ASP. ASP is considered as the sum of resin-extractable P, bicarbonate-extractable inorganic P and organic P. As our soil samples were obtained from the moderate weathering rocks, and the soil pH values of the soil samples are between 3.5 and 5.2, the Bray-I method was used for analysis of ASP [53].

3. Results

3.1. Spatial Pattern of NPP and NPP Trends

Figure 2a shows the spatial distribution of forest NPP averaged over 1950 to 2010 in the study area. The average NPP varies between 400 $\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 900$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, and the higher NPP values are mainly distributed in the southeastern and western forest areas, including forests in the Fujian, Jiangsu, Sichuan and Yunnan provinces. All forest NPP show positive trends during 1950–2010 (Figure 2b), and the Mann-Kendall test confirms that the increasing trends are all significant at 95% confidence level. The trends in NPP are in the range of 1–3$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$ in the study area, with the higher values (greater than 2$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$) in central subtropical China, including southwestern Sichuan, Guizhou, Yunnan and most of the Jiangxi province.

3.2. Changes of NPP for Different Forest Types

The changes of NPP time series for EBF, ENF and DBF are compared in Figure 3. The mean NPP averaged over 1950–2010 for EBF is 805 $\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, which is the highest among the three forest types, followed by 573$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ for DBF and 452$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ for ENF. The slope of EBF is 1.93 $\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$, which is the highest among those of the three forest types, while the slopes of ENF and DBF are 1.$33\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$ and 1.12 $\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$, respectively. The average NPP slope for the three forest types is 1.46$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$. Piao et al. (2011) reported the NPP slopes in Southern China were greater than 3$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ during 1980–2002 [54], and Gu et al. (2017) estimated an NPP slope of 2$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$ in southern China during 1960–2010 [15]. The study period for Gu et al. (2017) [15] is close to that of this study, and their estimates are close to our estimates for EBF, but slightly greater than our estimates for ENF and DBF.

3.3. Spatial Pattern of Additional Accumulated Phosphorus Demand

The spatial distribution of AAPD for the period of 1951–2010 for ${\mathrm{S}}_1$ is similar as that for ${\mathrm{S}}_2$, but differs from ${\mathrm{S}}_2$ in values. The AAPD are between 0.064–0.3$\mathrm{gC}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_1$ and between 0.005–0.064$\mathrm{gC}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_2$ (Figure 4). The higher values mainly distributed in central subtropical China, including northeastern and southwestern Sichuan, and northwestern Yunnan, with values between 0.25 and 0.30$\mathrm{gC}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_1$ and 0.016 and 0.064$\mathrm{gC}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_2$.

Figure 5 compares the AAPD for the first (1951–1980) and the latter (1981–2010) 30 years for ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$. For both ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$, the AAPD in the latter 30 years was mostly higher than that in the first 30 years and the higher AAPD values (greater than 0.128$\mathrm{gC}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_1$ and greater than 0.016$\mathrm{gC}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_2$) sporadically distributed in the western subtropical China in the first 30 years, but occupied large central and eastern subtropical regions in the latter 30 years.

We take ${\mathrm{S}}_1$ as an example to illustrate the changes of the spatial pattern of AAPD. For the first 30 years, the AAPD was below 0.032 $\mathrm{gC}·{\mathrm{m}}^{-2}$ in 16% of the total forest area, and was in the range of 0.032–0.064 $\mathrm{gC}·{\mathrm{m}}^{-2}$ and 0.064–0.256 $\mathrm{gC}·{\mathrm{m}}^{-2}$ in 47% and 37% of the total forest area, respectively. However, for the latter 30 years, only 1% of the total forest area was covered by AAPD below 0.032 $\mathrm{gC}·{\mathrm{m}}^{-2}$, while 18% and 80% of the total forest areas were covered by AAPD in the range of 0.032–0.064 $\mathrm{gC}·{\mathrm{m}}^{-2}$ and 0.064–0.256 $\mathrm{gC}·{\mathrm{m}}^{-2}$, respectively.

3.4. Changes of Additional Phosphorus Demand for Different Forest Types

The additional P demand averaged over forest types increased greatly during 1951–2010 (Figure 6). Linear regression analyses show that the slopes of additional P demand for EBF, ENF and DBF are 3.52, 1.85 and 0.81 $\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ for ${\mathrm{S}}_1$, and are 0.35, 0.18 and 0.08 $\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ for ${\mathrm{S}}_2$, respectively. The positive trends in the additional P demand for EBF, ENF and DBF are all significant, at 99% confidence level.

The AAPD values in the period of 1981–2010 are much greater than those in 1951–1980 for both EBF and ENF. Take ${\mathrm{S}}_1$ for EBF as an example; the AAPD values are 85.7 $\mathrm{mg}·{\mathrm{m}}^{-2}$ and 142.77$\mathrm{mg}·{\mathrm{m}}^{-2}$ for 1951–1980 and 1981–2010, respectively, and the latter is 1.67 times more than the former. For ENF, the AAPD for 1981–2010 is 1.8 times more than that for 1951–1980. The exceptional case occurred in DBF. The AAPD for DBF is 64$\mathrm{mg}·{\mathrm{m}}^{-2}$ during 1951–1980 and is 26$\mathrm{mg}·{\mathrm{m}}^{-2}$ during 1981–2010, which is only 0.41 times the former, implying a decreased growth rate in the latter 30 years. However, DBF occupies very small areas (1.38%) of subtropical China, and the decreased growth rate in AAPD is not common in this region.

3.5. The Latitudinal Distribution of Additional Phosphorus Demand in Subtropical China

We define the zones in the north of 30° N, between 25° N and 30° N, and in the south of 25° N in the study area as northern, central and southern subtropical China, respectively. The NPP trends averaged over northern, central and southern subtropical China for EBF are 1.63$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, 2.02 $\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 1.79$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, respectively, and those for ENF are 1.23$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, 1.38$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 1.27$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, respectively. For both EBF and ENF, the highest NPP trends occurred in central subtropical China, and the lowest in northern subtropical China, and the NPP trend for EBF was higher than that for ENF in the whole subtropical China.

The AAPD trends over 1951–2010 for EBF in northern, central and southern subtropical China are 3.44$\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, 4.25$\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 3.91$\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, respectively, and those for ENF are 2.00$\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, 2.14$\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 1.91$\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, respectively. The latitudinal distribution of AAPD trends is consistent with that of NPP trends. For both EBF and ENF, the highest AAPD trend occurs in central subtropical China, and the lowest in northern subtropical China. While the AAPD trends for EBF are higher than those for ENF in the whole subtropical China. Therefore, the central subtropical forests, especially the EBF, have the highest increasing trends of NPP and P demands over 1951–2010 due to climate warming.

3.6. Observed Soil Available Phosphorus in Forest Sites

Based on the field investigations conducted in the summer of 2018, the observed ASP in Shennongjia, Qian Yanzhou, Mt.Gongga and Dinghushan are 1.24$\mathrm{mg}·{\mathrm{kg}}^{-1}$, 1.43$\mathrm{mg}·{\mathrm{kg}}^{-1}$, 5.24$\mathrm{mg}·{\mathrm{kg}}^{-1}$ and 3.15$\mathrm{mg}·{\mathrm{kg}}^{-1}$ for EBF, and 1.63$\mathrm{mg}·{\mathrm{kg}}^{-1}$, 1.45$\mathrm{mg}·{\mathrm{kg}}^{-1}$, 1.98$\mathrm{mg}·{\mathrm{kg}}^{-1}$ and 1.53$\mathrm{mg}·{\mathrm{kg}}^{-1}$ for ENF, respectively. The spatial variations of ASP range from 1.24 $\mathrm{mg}·{\mathrm{kg}}^{-1}$ to 5.24 $\mathrm{mg}·{\mathrm{kg}}^{-1}$ among forest sites for EBF, which are greater than that for ENF. The latter range from 1.45 $\mathrm{mg}·{\mathrm{kg}}^{-1}$ to 1.98 $\mathrm{mg}·{\mathrm{kg}}^{-1}$ among forest sites. Based on the soil bulk densities observed in corresponding plots, ASP could be changed to ASP density (

Table 2. Observed available soil phosphorus (ASP) ($\mathrm{mg}·{\mathrm{m}}^{-2}$) in 2018 and the accumulated additional phosphorus demand (AAPD) ($\mathrm{mg}·{\mathrm{m}}^{-2}$ ) for the period of 1951–2010 for ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$.

Forest Site	ASP		$$\mathbf{AAPD}(\mathbf{AAPD}/\mathbf{ASP})\mathbf{for}{\mathit{S}}_1$$		$$\mathbf{AAPD}(\mathbf{AAPD}/\mathbf{ASP})\mathbf{for}{\mathit{S}}_2$$
	EBF	ENF	EBF	ENF	EBF	ENF
Shennongjia	812	1068	173 (21.4)	99 (9.3)	17 (2.1)	10 (0.9)
Qian Yanzhou	937	950	239 (25.5)	120 (12.7)	24 (2.6)	12 (1.3)
Mt. Gongga	3432	1297	231 (6.7)	113 (0.87)	23 (0.7)	11 (0.1)
Dinghushan	2063	1001	180 (8.7)	86 (0.86)	18 (0.9)	9 (0.1)

The numbers in the brackets are the values of AAPD/ASP (%). EBF and ENF are the abbreviations for evergreen broad leaf forest and evergreen needle leaf forest, respectively. ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$ are the scenarios for the low and high mineralization rates, respectively.

). The ASP densities in Shennongjia, Qian Yanzhou, Mt.Gongga and Dinghushan are 812 $\mathrm{mg}·{\mathrm{m}}^{-2}$, 937 $\mathrm{mg}·{\mathrm{m}}^{-2}$, 3432 $\mathrm{mg}·{\mathrm{m}}^{-2}$ and 2063 $\mathrm{mg}·{\mathrm{m}}^{-2}$ for EBF and 1068$\mathrm{mg}·{\mathrm{m}}^{-2}$, 950$\mathrm{mg}·{\mathrm{m}}^{-2}$, 1297$\mathrm{mg}·{\mathrm{m}}^{-2}$ and 1001$\mathrm{mg}·{\mathrm{m}}^{-2}$ for ENF, respectively. The variations of ASP density among forest sites for EBF are also larger than those for ENF. In addition, the differences in soil bulk densities in each forest plot also contribute to the ASP density differences among forest sites and forest types.

3.7. Past Phosphorus Consumptions in Forest Sites

To estimate the past P consumptions in the four representative forest sites, we calculated the AAPD for the period of 1951–2010 and their proportions to ASP and AAPD/ASP (Table 2). The AAPD are greater in EBF than in ENF in all sites. For ${\mathrm{S}}_1$, the AAPD and AAPD/ASP in Qian Yanzhou are 239$\mathrm{mg}·{\mathrm{m}}^{-2}$ and 25.5% for EBF, and 120$\mathrm{mg}·{\mathrm{m}}^{-2}$ and 12.7% for ENF, which are the highest among the forest sites. Followed by Shennongjia, with an AAPD and AAPD/ASP of 173$\mathrm{mg}·{\mathrm{m}}^{-2}$ and 21.4% for EBF, and 99$\mathrm{mg}·{\mathrm{m}}^{-2}$ and 9.3% for ENF. The AAPD/ASP in Mt. Gongga and Dinghushan are relatively smaller than that in Qian Yanzhou and Shennongjia. The AAPD and AAPD/ASP for ${\mathrm{S}}_1$ are nine to ten times those for ${\mathrm{S}}_2$.

3.8. Prediction of Future Phosphorus Limitation in Forest Sites

To reveal if the past increasing trend in additional P demand will continue in subtropical China, we used the outputs of BCC-CSM 1.1 to estimate the future AAPD in the four forest sites. The temperature will be increased by 4 °C, 3.4 °C, 2.9 °C and 2.3 °C, while the NPP will be increased by 327$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, 368$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$, 371$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 381$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ in Shennongjia, Qian Yanzhou, Mt. Gongga and Dinghushan, respectively, by 2099 under RCP8.5, compared with that in 2018 (

Table 3. Predicted temperature and NPP difference as well as P limitation in forest sites in subtropical China under RCP8.5.

Forest Site	ΔT (°C)	ΔNPP $(\mathbf{g}\mathbf{C}·{\mathbf{m}}^{-2}·{\mathbf{a}}^{-1})$	$\mathbf{AAPD}(\mathbf{m}\mathbf{g}·{\mathbf{m}}^{-2})$
			${\mathbf{S}}_1$	${\mathbf{S}}_2$
Shennongjia	3.99	327	568	57
Qian Yanzhou	3.41	368	644	65
Mt. Gongga	2.91	371	651	65
Dinghushan	2.32	381	668	67

$\Delta$T and $\Delta$NPP represent temperature and NPP difference between 2099 and 2018. The accumulated additional phosphorus demand (AAPD) is for the period of 2018–2099, and ASP is the available soil phosphorus observed in 2018. ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$ are the scenarios for the low and high mineralization rates, respectively.

). Future climate warming will stimulate the increase of NPP at an average rate of 4.2$\mathrm{gC}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-2}$ in the four sites, which are more than double that during 1951–2010. For the four sites, future AAPD are estimated in the range of 568–668$\mathrm{mg}·{\mathrm{m}}^{-2}$ and 57–67$\mathrm{mg}·{\mathrm{m}}^{-2}$.

The AAPD/ASP indicates the proportions of P consumption to P supply. The higher the values, the higher the degree of potential P limitation. By the end of the 21st century, for RCP8.5, the AAPD/ASP values for the four forest sites are predicted to be in the range of 32.2% to 68.4% under ${\mathrm{S}}_1$, and in the range of 3.3% to 6.3% under ${\mathrm{S}}_2$. For RCP4.5, the AAPD/ASP values for the four forest sites are predicted in the range of 12.8% to 32.1% under ${\mathrm{S}}_1$, and in the range of 1.3% to 3.2% under ${\mathrm{S}}_2$ (

Table 4. Comparison of AADP/ASP (%) for RCP8.5 and RCP4.5 among forest sites for the low and high mineralization rate scenarios.

Scenario	Shennongjia	Qian Yanzhou	Mt. Gongga	Dinghushan
RCP4.5 (${\mathrm{S}}_1$)	32.1	24.7	12.8	27
RCP8.5 (${\mathrm{S}}_1$)	62.2	68.4	32.2	47
RCP4.5 (${\mathrm{S}}_2$)	3.2	2.5	1.3	2.7
RCP8.5 (${\mathrm{S}}_2$)	6.3	6.9	3.3	4.7

${\mathrm{S}}_1$ and ${\mathrm{S}}_2$ are the scenarios for the low and high mineralization rates, respectively. The accumulated additional phosphorus demand (AAPD) is for the period of 2018–2099, and ASP is the available soil phosphorus observed in 2018.

).

Compared with the predicted AAPD/ASP in RCP4.5, the predicted AAPD/ASP in RCP8.5 is greater in each forest site. For both S₁ and S₂, from RCP8.5 to RCP4.5, the AAPD/ASP decreased by about 48.3%, 63.8%, 60.2% and 42.6% in Shennongjia, Qian Yanzhou, Mt. Gongga and Dinghushan, respectively.

4. Discussion

In this study, we analyzed the spatial and temporal variations of model-simulated NPP during 1951–2010 in subtropical China. Based on the variations of NPP, we estimated the AAPD over 1951–2010. We predicted future P limitation based on field observations of current ASP and model simulated future NPP. Our results demonstrated that climate warming has caused the increase of NPP and P demands in subtropical forests in China, and future climate change will likely continue to increase NPP and enhance P limitations in subtropical China.

4.1. Factors Affecting Phosphorus Limitation in Subtropical China

Though the additional P demand over 1951–2010 varied among forest sites and forest types, the P limitations are largely dependent on current ASP. For example, the AAPD for EBF in Mt. Gongga is higher (231$\mathrm{mg}·{\mathrm{m}}^{-2}$), but the AAPD/ASP there is the lowest (6.7%) among the EBFs in the four sites (Table 2), which indicates the importance of the ASP in determining the intensity of P limitation, and hence, it is necessary to investigate the spatial distribution of ASP in subtropical China.

Our estimates of ASP density show a large variety among sites and forest types. The ASP density is controlled by many factors, such as the parent material, soil weathering regime, biota, and the relative importance of these controls during soil development [51]. Therefore, the large variety in observed ASP density is understandable. Cao et al. (2014) reported that the observed ASP densities in the central subtropical region (109°45′ E, 26°50′ N) for ENF are in the range of 620–840 $\mathrm{mg}·{\mathrm{m}}^{-2}$ [55], which are lower than our estimate of the ASP density for ENF in Qian Yanzhou. Zhang et al. (2005) demonstrated that the average ASP in tropical and subtropical China is 4800$\mathrm{mg}·{\mathrm{m}}^{-3}$ in the top 50 cm soil [51], which could be changed to a density of 2400 $\mathrm{g}·{\mathrm{m}}^{-2}$ for comparison in this study, and is close to our estimates of ASP density in Mt. Gongga.

The AAPD/ASP for both EBF and ENF in Qian Yanzhou are the highest among those in the forest sites, and are due to higher AAPD and lower ASP. The forests in Qian Yanzhou are plantations, which implies human activities may have decreased the ASP there. According to the eighth national inventory of forest resources, there are 41.06 million hectares of plantations in subtropical China, which account for 19.8% of China’s total forest area. Cao et al. (2014) reported an even lower (less than 840 $\mathrm{mg}·{\mathrm{m}}^{-2}$) ASP in plantations in the central subtropical region near Qian Yanzhou [55], which implies lower ASP may widely exist in the subtropical plantations.

4.2. Comparison of Past and Future Phosphorus Limitations in Subtropical China

By the end of this century, the model-predicted NPP in the four forest sites will continue to increase under RCP8.5. To compare the intensity of P limitation in 1951–2010 with that in 2018–2099, we calculated the AAPD/$\Delta t$ in the two periods, where ∆t is the length of the time period. The AAPD/$\Delta t$ averaged over the four sites for ${\mathrm{S}}_1\mathrm{and}{\mathrm{S}}_2$ are, respectively, 2.58 $\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 0.26 $\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ for 1951–2010, and 7.62 $\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ and 0.77 $\mathrm{mg}·{\mathrm{m}}^{-2}·{\mathrm{a}}^{-1}$ for 2018–2099. The AAPD/$\Delta t$ for 2018–2099 is about 2.95 times more than that for 1951–2010, which indicates P limitations will be greatly enhanced in subtropical forests in China under RCP8.5.

The intensity of P limitation in Qian Yanzhou is the highest among the forest sites for the period of 1951–2010, and the future AAPD/∆t in Qian Yanzhou is also the highest, which implies the plantations in central subtropical China are likely to face the highest risks of P limitation in the future.

4.3. Comparison of Phosphorus Limitations between Low and High Mineralization Rate Scenarios

Mineralization processing is an important source of P supply. The mechanism affecting the mineralization process of organic P is complex. In this paper, the ${\mathrm{S}}_1$ scenario assumes that only the litters of leaf could be mineralized, and the ${\mathrm{S}}_2$ scenario assumes that all litters, including root, stem and leaf, could be mineralized. The estimated AAPD over 1951–2010 for ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$ are quite different. Take ENF as an example, the AAPD averaged over ENF were estimated as 160 $\mathrm{m}gP·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_1$, and 16 $\mathrm{m}gP·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_2$. From ${\mathrm{S}}_1$ to ${\mathrm{S}}_2$, the AAPD decreased by 90%. Therefore, the intensity of the mineralization rate is extremely important for determining the degree of P limitation in subtropical China. The human-managed subtropical plantations, with wood harvest, may correspond to ${\mathrm{S}}_1$, and the natural subtropical forests without human disturbance are more likely to correspond to ${\mathrm{S}}_2$. The observed ASP in Qian Yanzhou plantations is the lowest among forest sites, which confirms that human activities may decrease forest ASP.

4.4. Comparison of Phosphorus Limitation with Other Studies

Previous studies on P limitation revealed that the major sources of uncertainty are related to the estimates of ASP, the stoichiometry of plant organs and the proportion of NPP allocation to plant organs [33]. In this study, we obtained biome-specific C:P ratios for leaves, wood and roots, and NPP allocation to C pools (leaves, wood and roots) from literatures which reported field measurements in subtropical China. Furthermore, we conducted field observations in representative forest sites to estimate ASP, which could reduce the uncertainty of P limitation studies in subtropical China.

If we define the P deficit at the end of this century as the difference between ASP in 2018 and the AAPD for the period of 2018–2099, then the P deficit will be 372, 300, 1714 and 865 $\mathrm{mg}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_1$, and 883, 879, 2299 and 1465 $\mathrm{mg}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_2$ for RCP8.5 in Shennongjia, Qian Yanzhou, Mt. Gongga and Dinghushan, respectively. Sun et al. (2017) estimated the P deficit in these areas were in the range of −25,000 $\mathrm{mg}·{\mathrm{m}}^{-2}$ to 15,000 $\mathrm{mg}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_1$, and 5000 $\mathrm{mg}·{\mathrm{m}}^{-2}$ to 25,000 $\mathrm{mg}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_2$ [33]. Our study refines this estimation to the range between 300 and 1714 $\mathrm{mg}·{\mathrm{m}}^{-2}$ and 879 and 2299 $\mathrm{mg}·{\mathrm{m}}^{-2}$ for ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$, respectively. In this study, we assumed 1950 as the year of steady status, which may lower our estimates of P deficit when compared with other studies which assumed 1860 as the year of steady status.

4.5. Limitations and Further Studies

While great efforts have been made to reduce the uncertainties in the diagnosing of the P limitations in subtropical China, there are still some sources of uncertainty that were not taken into account in this study. First, the synergistic effect of N and P is neglected, which may result in higher estimation of AAPD in N-limited areas. Second, ecosystem-level manipulation experiments have shown that climate warming will gradually change the C:P ratios in ecosystems [1,56,57], while human activities, such as fertilization and pollution, also changes C:P ratios in ecosystems [1,56,57]. However, the C:P ratios for plant organs are assumed to be constant during the study periods. Third, the allocations of vegetation carbon to roots, stem and leaves are assumed to be constant, but are changed with nutrient availability [58].

Further studies are needed to better understand and quantify the effects of N-P interactions, dynamic C:P ratios and vegetation carbon allocations on the P limitation. A comprehensive sensitivity analysis would help us to improve our understanding of the uncertainties, and will be an important part of our future research. In addition, the distribution of current ASP in subtropical China has large heterogeneity [59], and is a major factor affecting the degrees of P limitation, and, hence, further studies need to conduct more field investigations to understand the spatial distribution of ASP and P limitation in subtropical China.

5. Conclusions

We diagnose P limitation in subtropical forests in China under climate warming based on model-simulated NPP. The forest NPP in subtropical China increased significantly in recent decades, and the P limitation has been strengthened due to climate warming. The AAPD in 1981–2010 are 1.67 and 1.8 times that in the former 30 years for EBF and ENF, respectively. The higher AAPD sporadically distributed in the western subtropical China for the period of 1951–1980, but covered large areas of central and eastern subtropical China for the period of 1981–2010.

The observed current ASP density in forest sites varied from 812 to 3432 $\mathrm{mg}·{\mathrm{m}}^{-2}$ in EBF and 950 to 1297 $\mathrm{mg}·{\mathrm{m}}^{-2}$ in ENF, and the AAPD/ASP for EBF were greater than those for ENF. The AAPD for ${\mathrm{S}}_1$ is about 10 times that for ${\mathrm{S}}_2$. The AAPD/ASP in the four sites varied between 6.7% and 25.5% for EBF, and between 0.86% and 12.7% for ENF for the period of 1951–2010 for ${\mathrm{S}}_1$.

Future climate change will enhance P limitation in the four sites. The increasing trends of forest NPP will be doubled, and the AAPD/ASP are estimated to be in the range of 32.2–68.3% and 3.2–6.9% for ${\mathrm{S}}_1$ and ${\mathrm{S}}_2$, respectively, under RCP8.5 at the end of this century.

Though climate warming is one of the major factors in increasing the intensity of P limitation in subtropical forest in China, human disturbance is also another important factor. The ${\mathrm{S}}_1$ corresponds to low recycling of P by plants, which is more likely to correspond to forests with human disturbance. Therefore, future P limitation in subtropical China forest to some degree depends on human activities, and intense human disturbance will strengthen P limitation. In addition, the amount of current ASP is another factor that affects the intensity of P limitation in subtropical China, and understanding the current ASP distribution is essential for understanding the intensity of P limitation in subtropical China.