Soil phosphorus crisis in the Tibetan alpine permafrost region

Domain: earth_science
Published: 2025-07-05
Editorial outcome: minor_revision

Abstract
Phosphorus (P) is an essential nutrient for living systems and is critical to the functioning of ecosystems. Permafrost areas have a huge reservoir of soil P that is currently not used very much; however, the direction and magnitude of changes in soil P stocks across the Tibetan alpine permafrost region over recent decades remain unclear and the P budget has not been well assessed. Here we use a unique combination of a soil resampling method and a modified process-balanced model to assess the historical dynamics of soil P pools (0–30 cm depth) and the key flows of P in ecosystems across Tibetan alpine permafrost region. Compared with the 1980s, the soil P stock decreases dramatically by 36.1% in the 2020s, decreasing from 346.5 to 221.4 Tg P (1 Tg = 10^12 g) during the last three decades. Water erosion accounts for 82.3% of the total soil P outflow. Our projections suggest that the soil P stock will only be 20.3% of the 1980s stock by the end of this century, leading to an unprecedented crisis of P limitation in permafrost regions.

Introduction
Phosphorus (P) is a key element in organic molecules (e.g., nucleic acids, ATP and phospholipids) and has an important role in regulating the growth, development and reproduction of living organisms1–4. Permafrost regions occupy about 24% of the land surface in the northern hemisphere5, but their soil P reserves account for >30% of the northern hemisphere total (soil depth of 0–30 cm, calculated according to ref. 6). Permafrost soils contain a huge reservoir of P as a result of the rejuvenation of soils, weak leaching processes and the low uptake of nutrients by plants, and can therefore be regarded as important reserves of P in terrestrial ecosystems6.

The recent rapid changes in the Earth’s climate are simultaneously causing the widespread degradation of permafrost7 and altering the hydrology of high-latitude/high-altitude regions8, exposing mineral nutrients previously protected in permafrost to transport and processing in lakes, rivers, and estuaries9,10. The unidirectional flow of P from terrestrial ecosystems to freshwater bodies and, eventually, the oceans is considered to be the greatest sustainability challenge in the P cycle11. The redistribution of P on spatiotemporal scales could regulate the composition of local, downstream, and even coastal plant communities, as well as the net primary productivity and health of the environment12,13. In addition, the availability of P alters the trajectory of the permafrost carbon (C) cycle by influencing gross primary productivity and net ecosystem productivity14,15, and the associated fluxes of C in aquatic sediments12 by influencing soil microbial processes and greenhouse gas emissions. This can have a positive or negative feedback role in global climate change. Therefore, quantifying the temporal changes of the P stored in permafrost is crucial and urgently needed in guiding strategies to tackle the grand challenge of future P limitation11,16, as well as the environmental crisis caused by huge losses of P into water bodies that pushes the Earth’s environment to its boundary, beyond which come abrupt environmental changes with catastrophic impacts17,18.

Several approaches have been used to explore the historical dynamics of soil P, including soil sample analysis and P budget models11,19,20. The resampling method can accurately estimate historical changes in the soil pool of elements at large spatial scales21–23. Process-balanced modeling approaches have been used to detect flows and their drivers in cropland soils over time1,2, demonstrating that water erosion11, absorption by plants19, food production2,20 and/or the recycling of animal manure24 might have a major impact on the soil P balance. Despite advances in our understanding of the P budget over recent decades, substantial uncertainties remain about how much P is lost from soils1,11,19,25,26. The main cause of these discrepancies is that our understanding of current processes and historical soil nutrient availability data are lagging behind the needs of mass balance models27.

As the largest extent of alpine permafrost in the world28, the Tibetan Plateau is an important pool of mineral nutrient elements globally, with 25% of the material flux reaching the ocean coming from this area, despite the fact that it represents only 5% of the Earth’s land surface6,29. A few studies have used observations from non-resampling methods to explore the historical changes in total P30 and available P (AP)31 in alpine soils. However, the quantity of SPS that have been lost across Tibetan alpine permafrost regions over recent decades remains unclear, and the P budget of permafrost has not been well assessed. Moreover, substantial uncertainties exist in the projected direction and magnitude of soil P changes, originating from the scarcity of observational broad-scale soil resampling data.

In this study, we used a unique combination of soil resampling and model simulations to assess the temporal changes of soil P pools and the key flows of P in the Tibetan alpine permafrost region (Fig. 1). First, we performed a broad-scale repeat soil survey to obtain data for the soil P density (SPD) in the 1980s and the 2020s. We then built a process-based model of the P flow in ecosystems (MPFE) to simulate changes in the soil P budget from the 1980s to the 2020s according to the mass balance principle. A notable difference between this work and previous studies is that the soil P data were repeatedly measured in situ for the two time periods to verify and improve the accuracy of the MPEF. We predicted the soil P budget up to the year 2100 using the MPEF according to the Shared Socioeconomic Pathway (SSP) SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 future climate scenarios, providing scientific support for predicting future P balances in permafrost regions.

Results
Rapid decreases in soil P since the 1980s
We built a unique dataset composed of two parts separated by three decades of time to evaluate the spatiotemporal changes in the density and stock of P in topsoil (0–30 cm depth) across the alpine permafrost regions of the Tibetan Plateau (Supplementary Figs. 1, 2, Supplementary Method 1 and Supplementary Table 1). There was a significant positive correlation between the SPD in the 1980s and 2020s (R2 = 0.32, P < 0.001), indicating that the spatial pattern of the SPD had not changed significantly in the past three decades (Fig. 2a, b). However, the SPD showed a significant decrease across the alpine permafrost regions from the 1980s to the 2020s, with the mean SPD decreasing from 1583.8 (±67.6 SE) to 1046.4 (±32.0 SE) kg P ha−1 (P < 0.001) (Fig. 2c). Further analyses revealed that the SPD tended to decrease among 190 of the total 255 resampling plots and the decline in SPD remained significant until >240 resampling plots were removed from the initial analysis with 255 resampling plots (Supplementary Figs. 3 and 4), demonstrating the robustness of SPD dynamics observed in this study. This decreasing pattern held true in alpine grassland, forest and cropland ecosystems (P < 0.01), with the net loss of SPD estimated to be 533.5 kg P ha−1, 414.8 kg P ha−1, and 715.0 kg P ha−1, respectively (Fig. 2c, Supplementary Figs. 3 and 4). The rapid decrease in soil P (14.9 kg ha−1 yr−1) across the Tibetan Plateau was within the range of the global soil P depletion rate (4–19 kg ha−1 yr−1)11.

In terms of the spatial pattern of historical changes in the SPD, the largest decreases in SPD occurred in the central and northeastern Tibetan Plateau, whereas the smallest decrease in SPD was in the western Tibetan Plateau (Supplementary Fig. 3). This may be related to the spatial differences in the quantity of water-induced P loss (Supplementary Fig. 5).

The SPS of the Tibetan permafrost regions decreased dramatically by 36.1% in the 2020s compared with the 1980s, with a gross regional-scale change in the surface SPS of 125.1 Tg P (1 Tg = 10^12 g) over the last three decades (P < 0.01) (Fig. 2d). The average net soil P loss of the three ecosystems was approximately 3.5 Tg P yr−1, which is roughly four-fifths of the annual application of chemical P fertilizer to farmland in China (4.2 Tg P yr−1) and one-fifth of the annual application of chemical P fertilizer globally (16.4 Tg P yr−1) (Supplementary Fig. 6). Such a rapid depletion in SPS over a short period of time has not previously been detected and poses a potential crisis in the limitation of P in terrestrial alpine permafrost ecosystems.

Historically negative budget for soil P
The database of the MPFE model developed to simulate the historical trajectories of the P budget consists of the main key P inflow and outflow subsystems (e.g., weathering, atmospheric deposition, fertilizer, absorption by plants, the deposition of livestock excretions, the return of plant litter, livestock feed, livestock, forest and grain products, the burning of firewood and soil erosion) (Supplementary Method 2 and Supplementary Table 2). Similar to conventional modeling approaches11,19,20, we compared our calculations of the main inflow and outflow subsystems with those of other studies. Most of our modeled P flows were well within the range of published P flows (Supplementary Tables 3–5). To further test the accuracy of our model, we compared the observations and simulations of the change in the SPD between final and first years of observations (Fig. 3a, b). We found three important results: (1) there was fairly good agreement between the simulated and observed temporal changes in SPD at the plot scale (R2 = 0.87, P < 0.001) and county scale (R2 = 0.90, P < 0.001); (2) the results of the slope test indicated that the simulated and observational data either increased or decreased in equal proportions (P > 0.05); and (3) the mixed-effects model showed no significant difference between the geometric mean of the observed and modeled SPD (P > 0.05). For a few plots, however, the model had difficulty in capturing the observed temporal changes in SPS. This could have been related to the thin soil thickness used to calculate the SPS preventing us from accounting for vertical migration, or for the uptake of P by plant roots32.

The inflow density of P increased slightly from 4.5, 4.4, 4.7 and 16.2 kg P ha−1 yr−1 to 4.8, 4.6, 4.8, and 39.3 kg P ha−1 yr−1 for the overall ecosystem and the grassland, forest and cropland ecosystems, respectively, from the 1980s to the 2020s (P < 0.001) (Fig. 4a–d). Weathering made the highest contribution to the soil P inflow (70.4%) under the effects of the freeze–thaw cycle33 and glacial denudation34 (Supplementary Discussion 1). The P outflow density decreased significantly from 25.6, 26.4, and 20.6 kg P ha−1 yr−1 to 18.9, 19.4, and 15.4 kg P ha−1 yr−1 for the overall ecosystem and the grassland and forest ecosystems, respectively, from the 1980s to the 2020s (P < 0.01) (Fig. 4a–c). The P outflow density of the cropland ecosystem increased significantly from the 1980s to the 2020s (P < 0.01) (Fig. 4d). However, the soil P outflow far exceeded the soil P inflow, and therefore the total soil P budget in the alpine permafrost region of the Tibetan Plateau was negative during the last three decades.

Over the same time period, the absorption of P by plants and the dissolved P in runoff accounted for only a small proportion of the total P outflow (13.7%) (Fig. 4e–h, Supplementary Figs. 7 and 8). Among the sub-terms of the absorption of P by plants, the mean annual consumption of P by humans was only 0.3, 0.1, 0.3 and 38.6 kg P ha−1 yr−1 for the overall ecosystem and the grassland, forest and cropland ecosystems, respectively, between the 1980s and 2020s (Supplementary Fig. 9a–d). Since 2006, the implementation of energy substitution projects35 has led to a significant decrease in the consumption of P by burning animal dung and firewood (Supplementary Fig. 9 and Supplementary Discussion 1). The interannual change in dissolved P was not significant due to the spatial heterogeneity of precipitation across the alpine permafrost region of the Tibetan Plateau36, whereas the increased absorption of P by plants was caused by the recent favorable climatic conditions for plant growth21 (Supplementary Discussion 1).

Erosion caused the largest outflow of soil P across the Tibetan Plateau (86.3%) (Fig. 4e–h and Supplementary Fig. 10). Our estimates were further verified by other measurements, such as model simulation and isotope tracing methods (Supplementary Table 4), as well as the seven main regional riverine exports of P (Supplementary Table 5). Among them, the outflow of P caused by water erosion (82.3%) was far greater than that caused by wind erosion (4.0%) (Supplementary Fig. 10). Our area-related calculations resulted in an average outflow of P from the alpine permafrost region caused by water erosion of about 16.1 kg P ha−1 yr−1 (Supplementary Fig. 10a), which is slightly higher than that of previous calculations worldwide (Supplementary Table 6). The annual precipitation is concentrated in the summer months, and the complex conditions of differences in elevation, low vegetation coverage, and the weak resistance to erosion all lead to strong water erosion in this alpine permafrost region37.

The recent trend of warming and humidification has led to the degradation of glaciers and permafrost soils. The increased duration of soil thawing and the advance of the start date of soil thawing could intensify the erodibility of alpine soils and hence increase the supply of sediment38,39. Previous research on northern rivers has established marked increases in discharge and the oceanward flux of alkalinity and associated ions across the pan-Arctic region since the early to mid-20th century, which has been attributed to intensification of the hydrological cycle40. However, the increase in vegetation coverage caused by the warming and humidification on the Tibetan Plateau may reduce the intensity of erosion. The insignificant historical change in water erosion across the Tibetan Plateau may have arisen from the counteracting effects of rainfall erosivity factor and vegetation factor (Supplementary Fig. 11). Even if the intensity of soil erosion shows marginal historical changes, the outflow of soil P caused by erosion gradually decreases due to the continuous reduction in soil P concentrations. In addition, the widely distributed floating dunes and desertification in the northwest Tibetan Plateau, on the banks of the Yarlung Zangbo River and in the Qaidam Basin, provide abundant sources of material for wind erosion41. Dust aerosols can easily enter the middle and upper troposphere and be transported for long distances over East Asia and the Pacific Ocean due to a combination of the high terrain and intense atmospheric activities over the Tibetan Plateau42.

Potential scarcity of soil P by 2100
Our results highlight that the outflow of soil P will be higher than the inflow of soil P throughout the 21st century, although the loss of soil P gradually decreases under the SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5 scenarios (Fig. 5a, b and Supplementary Fig. 12). In general, the future outflow of soil P in the alpine permafrost region is still dominated by erosion and absorption by plants (Fig. 5c), whereas the amounts of P outflow from human consumption and dissolved P in runoff are relatively low (Supplementary Fig. 13). The accumulation of biomass increases the P demand of plants under improved growth conditions for alpine plants (caused by warming, humidification and CO2 enrichment) and extension of the length of the growing season21,43. During the next few decades, the amount of soil P absorbed by plants may become gradually closer to the outflow of soil P caused by erosion and become a non-negligible cause of soil P deficiency in permafrost areas.

The results of the change-point tests showed that soil P will continue to decrease rapidly from the 2020s to the 2050s, with an average rate of 10.0 kg P ha−1 yr−1 (Fig. 5d). The projected net loss of soil P will then gradually slow down, with the mean rate decreasing to 6.8 kg P ha−1 yr−1 by the end of the 21st century (Fig. 5d), mainly due to a decrease in the amount of P outflow via erosion. The predicted mean cumulative loss of the SPS for the four scenarios by the end of this century is 152.3 Tg P (Supplementary Fig. 14), which is just 20.3% of the 1980s stock, and thus highlights a serious challenge for future nutrient management and sustainable development in permafrost regions.

Methods
Study area
The Tibetan Plateau, with an average altitude higher than 4000 m, covers the whole of the Tibet Autonomous Region and Qinghai Province, some of Gansu, Sichuan, and Yunnan provinces, and the Xinjiang Uygur Autonomous Region. The Tibetan Plateau covers the area (25.99–39.82°N, 73.50–104.67°E) (Supplementary Fig. 1). The Tibetan permafrost accounts for 75% of the area of alpine permafrost in the northern hemisphere66. The Tibetan Plateau has experienced significant increases in the mean annual temperature, by 0.4°C decade−1 67, and precipitation has also shown a general upward trend during the last three decades43. An increase in warming, moisture-induced thawing of the permafrost, and increased erosion have all been reported in recent decades38,39. The two dominant natural ecosystems on the Tibetan Plateau are alpine grassland and forest68. The main herbaceous species in the alpine grassland are Kobresia pygmaea, Kobresia littledalei, Stipa purpurea, and Stipa glareosa69. The main woody species include Abies georgei var, Picea likiangensi, Picea crassifolia, Abies faxoniana, Picea brachytyla and Betula albo-sinensis70. The croplands are mainly located in the gullies and valleys of Huangshui, Yarlung Tsangpo River, Lancang River, Niyang River, Nianchu River, Yalong River, and Shiquan River71. As important ecological protection areas and key agricultural production areas across the Tibetan Plateau, they provide important guarantees for local grain self-sufficiency72. The main grain crops are the two most common crops on the plateau—namely, highland barley and wheat. Following the World Reference Base for Soil Resources classification, the major soil types in alpine grasslands are humic cambisols and haplic xerosols, while forests are largely characterized by humic cambisols, orthic luvisol, and dystric podzoluvisols, and croplands by haplic kastanozems and gelic leptosols73.

Urbanization on the Tibetan Plateau has been relatively slow during the last three decades (Supplementary Figs. 16 and 17). The Tibetan Plateau is still dominated by natural ecosystems, although the impervious area has increased from 5454 to 15,997 ha during the last three decades. In 2020, the impervious area of the Tibetan Plateau accounted for <0.007% of the plateau (Supplementary Fig. 17). Historically, farmers and herders have had their own unique lifestyles, consuming mainly beef, mutton, milk, fur, and grain. The dung from livestock (e.g., yaks, yellow cattle, Pien Niu, goats, and sheep) is burned as fuel for cooking and heating in the grassland and cropland areas, whereas firewood is burned for cooking and heating in the forested areas74,75. There are few areas of intensive pasture and most farmers are self-sufficient and follow traditional grazing and farming patterns.

Field resampling and laboratory measurements
We compiled a broad-scale repeat-sampling dataset based on two multi-year field campaigns in the 1980s and the 2020s to detect any changes in the SPD and stock across the alpine permafrost regions of the Tibetan Plateau. Historical soil data were obtained from the Second National Soil Inventory of China and the First Tibetan Plateau Scientific Expedition Program. The data from the Second National Soil Inventory of China were compiled by the Agriculture and Animal Husbandry Bureau or the Agricultural Bureau, whereas the data from the First Tibetan Plateau Scientific Expedition Program were compiled by the Tibetan Plateau Research Team, Chinese Academy of Sciences, during the 1980s (Supplementary Table 1).

Historical data for the physical and chemical properties of soils are reported in individual monographs for each province, city, and county. These monographs contain soil data of the typical profiles of different soil series and record detailed information, including the geographic allocation, records of landmarks, parent material, soil type, topographic features, land use type, vegetation type, the dominant species, and the physicochemical properties of each soil genetic horizon. The original soil P dataset comprised 605 plots in the 1980s. As a result of practical constraints (e.g., natural hazards, road rebuilding, and human disturbance; see Supplementary Method 1), a total of 255 perfectly matched resampling plots were obtained (Supplementary Fig. 1). Frequency distribution analysis of SPD revealed similar means between the original and resampled plots. A mixed-effects model further confirmed no significant difference between the two datasets (P > 0.05). The distribution and representativeness of the plots are depicted in Supplementary Fig. 2. Historical data in the 1980s were recorded from 131 counties across the Tibetan Plateau, and resampling plots were obtained from 93 counties in the 2020s, accounting for 71% of the total 131 counties (Supplementary Fig. 1b). To further evaluate the accuracy of our plots, we compared all the plots with other independent references. We made secondary use of different maps, including vegetation-type maps, soil-type maps, topographic maps (elevation, slope, and terrain), land-use maps (residential areas, rivers, and lakes), road maps, and administrative-boundary maps that were developed before 1990s (more details see Supplementary Method 1). In addition, we also consulted several staff who participated in the 1980s campaign to determine the exact locations of the sampling plots. Plots in close proximity to electric power lines, roads, railways, or drainage ditches were avoided. These plots were set throughout the geographical extent of alpine forests, grasslands and cropland, covering a wide range of climates (mean annual temperature: −12.6–14.3°C; mean annual precipitation: 22.4–1299.5 mm), vegetation characteristic (net primary productivity: 6.1–1173.7 g C m−2 yr−1; normalized difference vegetation index (NDVI): 0–0.8), soil properties (soil P concentration: 0.2–1.6 g kg−1; soil AP concentration: 0.4–97.6 mg kg−1) across the plateau. The topography, parent material, soil type, vegetation type, dominant species, and management strategies of the 255 resampling plots did not change during the intra-sampling period.

We randomly set up three (20×20 m) quadrats at intervals of 10 m and collected soil samples from the center of each quadrat for each plot (Supplementary Fig. 1c). At each sampling quadrat, a square pit covering an area of (1×1 m) was dug and a rule was placed inside the pit along the profile (repeated three times). We took three soil cores in each plot to generate descriptions of their vertical properties (e.g., layering, color, and texture for each layer) for comparisons with those from the 1980s. The consistency between each layer of the profile in the 2020s and the historical record was the key piece of evidence used to determine whether the soil sample was collected correctly. Undisturbed soil samples were collected from different soil depth intervals according to the soil genetic horizon. A 2 kg sample of each soil layer was placed into a Ziploc bag and labels were affixed inside and outside the bag, recording the profile number, location, depth, and other information. To test the soil bulk density, 100 cm3 steel cylinders were pushed into the soil profile with the center of the cylinder aimed at the central depth of each soil genetic horizon. If the target soil genetic horizon was thinner than the diameter of the steel cylinders, the adjacent layers shared the soil bulk density data of the thicker soil genetic horizon. Where large stones and tree roots were present, we moved the square pit horizontally (for all soil genetic horizons) until we could successfully collect soil samples from along the whole profile.

The soil samples collected in the 2020s were air-dried to a constant weight in a ventilated room (~15–20°C). The soil samples were then ground and passed through 2 and 0.149 mm sieves after any stones, plant litter, and roots had been removed. The P concentration of each soil sample was determined by the molybdenum antimony resistance to colorimetric method76, and the soil AP concentration was determined by the NaHCO3 extraction/Mo-Sb colorimetric method77. The two analysis methods were identical to the methods used in the 1980s. The samples used to determine the bulk density were dried to constant weight in an oven at 105°C. The mean P concentration and bulk density determined for three repeat samples from each soil layer were used in all further calculations. We selected the soil horizon corresponding to a thickness of 30 cm, which includes most of the plant roots and is the most active soil layer in the biogeochemical cycle of P78. Although there was a significant positive correlation between the P concentration in the two time periods at the same plot, the soil P concentration decreased significantly from the 1980s to the 2020s (Supplementary Fig. 18).

The SPD in the top 30 cm was calculated for each soil profile using the following equation:
SPD = Σ i=1 n Ti × BDi × SPCi / H (1)
where SPD is the SPD (kg P ha−1), T is the layer thickness (cm), BD is the bulk density (g cm−3), SPC is the soil P concentration (g P kg−1), i is a particular soil layer, and H is the number of profile layers within 30 cm.

For each resampled plot, the change in the SPD was calculated as
ΔSPD = SPD2020s − SPD1980s (2)
where ΔSPD represents the change in the SPD of the 0–30 cm depth soil layer during the last three decades (kg P ha−1), SPD2020s is the SPD of the 0–30 cm depth soil layer obtained in the 2020s, and SPD1980s is the SPD of the 0–30 cm depth soil layer obtained in 1980s. The calculation method for the soil AP density is consistent with that used for the SPD.

The soil P stock (SPS) was calculated according to the distribution area of different soil types and their corresponding SPD in grassland, forest, and cropland ecosystems because the SPD varies significantly in different soil types and even in the same soil type in different ecosystems. We then added the SPS of the three ecosystems to obtain the total SPS across the Tibetan Plateau. The soil type map was obtained from China’s soil atlas (scale 1:50,000)79, SoilGrid (250 m; https://www.soilgrids.org/), and the FAO (50 km; https://www.fao.org/). An ensemble calculation of the SPS was adopted to reduce the uncertainty caused by the three different resolutions of these soil maps. The vegetation maps of the 1980s68 and 2020s80 were obtained from China’s vegetation atlas (scale 1:1,000,000). The grassland ecosystem includes alpine meadow, steppe, desert, and cushion vegetation.

We compared the results of changes in the SPS during the two time periods. The change in SPS was calculated as
ΔSPS = SPS2020s − SPS1980s (3)
where ΔSPS represents the change in SPS of the 0–30 cm depth soil layer during the last three decades (kg P), SPS2020s is the SPS of the 0–30 cm depth soil obtained in the 2020s, and SPS1980s is the SPS of the 0–30 cm depth soil layer obtained in the 1980s.

Simulations of P flow in ecosystems
Our general conceptual framework for the soil P budget in grassland, forest and cropland ecosystems is based on the principle of mass balance (Fig. 1). The MPFE consists of the main key P inflow and outflow subsystems, including weathering, atmospheric deposition, fertilizer, absorption by plants, deposition via livestock excretions, litter return, livestock feed, livestock, forest and grain products, the burning of livestock excretion products (animal dung) for fuel, and soil erosion (Supplementary Method 2 and Supplementary Table 2). This framework has previously been widely applied to P flows on national and global scales1,11,32,81,82. Similar to conventional modeling approaches, we compared our calculations of the main P inflow and outflow subsystems with those of other studies. A notable difference between this work and previous studies is that the soil P data were repeatedly measured in situ for the two time periods to verify and improve the accuracy of the MPEF. In particular, we compared the observations and simulations of the change in the SPD between final and first years of observations (see Statistical analyses, below). The net change in soil P is the P inflow minus the P outflow:
Soil P budget = Soil P inflow − Soil P outflow (4)
Soil P inflow = P released from weathering + Atmospheric P deposition + P deposition via excretion + P return from plant litter + P fertilizer (5)
Soil P outflow = P absorption by plants + P loss via erosion + Loss of dissolved P (6)

Briefly, the inputs of P from weathering, atmospheric deposition and fertilizer are the P inflow subsystems, whereas the P loss from soil erosion (particulate P) and runoff (dissolved P) are the P outflow subsystems. However, the calculations for the plant-related P flows are more complex and include the net accumulation by plants, the return of plant litter to the soil, feeding by livestock, the deposition of livestock excretions, and human consumption (e.g., livestock, forest and grain products, and the burning of livestock excretion products for fuel)19,82 (Supplementary Figs. 19–21, Supplementary Data 1, 2). All the inflows and outflows of P were calculated at a resolution of an 8×8 km grid and are given in Supplementary Method 2.

Supplementary Table 3 compares the estimates of the main plant-related P flows on the Tibetan Plateau with those reported by other researchers. Although it is not easy to compare all the results obtained in this study with those reported previously due to differences in the system boundaries, approaches, and scope, it is possible to compare some of the flows calculated here with a few other studies conducted on the Tibetan Plateau and in China and Asia. Most of our modeled P flows were of the same order of magnitude as previous P flows (Supplementary Table 3). However, previous studies found larger amounts of P stored in livestock products, primarily due to the large number of yak and sheep in the alpine meadow ecosystem, whereas our study covered all types of alpine grasslands. The sparse population, low vegetation coverage and small number of livestock in the central and western areas of the plateau (i.e., the alpine steppe and alpine desert) resulted in a low P density for livestock products per unit area. The carcass weight, which is the weight of an animal after slaughter and the removal of the fur, head, and genitourinary organs, is the common indicator of livestock meat production. In the present study, however, the total amount of P in livestock products was the sum of the amount of P in meat, milk, and livestock byproducts (blood, adipose tissue, skeleton, viscera, digestive content, and skin hair). There are large differences in the P concentration in different fractions, such as the skull, gut, and blood1. Such differences in calculation methods may lead to differences in the P outflow of livestock products.

The Revised Universal Soil Loss Equation (RUSLE) and the Revised Wind Erosion Equation (RWEQ) models were used to calculate the erosion of soil by water and wind, respectively (Supplementary Method 2). To verify the accuracy of the estimation of the RUSLE and RWEQ models, our results were compared with other measurements (e.g., model simulation, isotopic tracing, and sand collection from the field). Our estimates of the P flows caused by erosion were of the same order of magnitude as previous P flows (Supplementary Table 4). For water erosion, the calculated potential soil P loss in this study was well within the range of published on-site riverine P exports for the Yangtze, Yellow, Lancang, Nujiang, Yarlung Tsangpo, Heihe, and Shule rivers (Supplementary Table 5). Although the riverine P loads of a few hydrological stations may not completely cover the P soil erosion status of the whole basin, they can provide an important reference for verifying the transport of P by sediments. Although the P transported from soils might be redeposited within water catchments (i.e., wetlands, lakes, ponds, riparian zones or buffer strips), the potential threats of P loss from soils are not reduced11.

Future scenario for the time period 2021–2100
Future climate data projected by 27 global climate models (GCMs) of the Coupled Model Intercomparison Project Phase 6 (CMIP6) from 2021 to 2100 were downloaded from the World Climate Research Program (WCRP) dataset (https://esgf-node.llnl.gov/search/cmip6/). The raw GCMs from the WCRP dataset were at low spatial and temporal resolutions and therefore downscaling and bias-correction methods were used to generate the site-scale daily datasets83. Four SSP scenarios (including the SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5 future climate scenarios of the CMIP6 models) were selected for projecting future soil P budgets (Supplementary Table 7). The SSP2-4.5 and SSP5-8.5 scenarios were applied to predict future P flows on the Tibetan Plateau from the 27 GCMs of CMIP6—namely, ACCESS-CM2, ACCESS-ESM1-5, BCC-CSM2-MR, CanESM5, CanESM5-CanOE, CNRM-CM, CNRM-CM6-1-HR, CNRM-ESM, EC-Earth3, EC-Earth3-Veg, FGOALS-g3, GFDL-ESM4, GISS-E2-1-G, INM-CM4-8, INM-CM5-0, IPSL-CM, MIROC6, MIROC-ES2L, MPI-ESM1-2-HR, MPI-ESM1-2-LR, MRI-ESM, UKESM1-0-LL, CMCC-CM2-SR5, CIESM, GFDL-CM4, MRI-ESM2-0 and NESM3. The SSP1-2.6 and SSP3-7.0 scenarios do not include the CMCC-CM2-SR5, CIESM, GFDL-CM4, and NESM3 models. The data were integrated using the multi-model ensemble method to reduce the uncertainty caused by an individual GCM. The average result from multiple GCMs is widely used in climate projections based on the assumption that the errors tend to cancel each other out if the GCMs are independent84–86.

Future population data were derived from gridded population and economy datasets under the four SSP scenarios87. To verify the accuracy of future predictions under each scenario, we standardized the meteorological (1985–2020) and population (2010–2020) data during the common period. Based on the linear relationship between the data series before and after this treatment, we converted the modeled result to a historical result and then scaled the future data accordingly88.

Statistical analyses
Descriptive statistics included the frequency distributions, means, and 95% confidence intervals of the means of the SPD at soil depth of 0–30 cm for both time periods. The relationship between the SPD in the 1980s and 2020s was analyzed using the ordinary least-squares (OLS) regression model. The temporal changes in SPDs were examined with linear mixed-effects models using the lme489 and the lmerTest90 packages in R version 4.2.1, in which the fixed effect was the sampling period and the random effect was the sampling plot. The effects of the sampling period on the SPD were considered significant if P < 0.05. The homoscedasticity and normality of the residuals were tested before the linear mixed modeling analyses and in-normal transformation was conducted where necessary. The historical trend of the outflow and inflow of the SPD from the 1980s to the 2020s was analyzed by the OLS regression model.

To test the accuracy of the dynamic model for P, a simulation vs. observation comparison of the change in the SPD between the observations in the 1980s and the 2020s was analyzed at both the sampling plot and county scale. First, the OLS fit of the linear equation was applied to calculate the relationship between the simulated vs. observational SPD data for the two sampling periods. The coefficient of determination can reflect the correlation of the SPD between the simulated and actual sampling data. Second, the slope of the fitted equation was compared with the 1:1 line using the tidyverse package to determine whether the simulated and observed data increased or decreased in equal proportions91. A mixed-effects model was then conducted using the lme489 and lmerTest90 packages in R version 4.2.1 to test for possible random effects due to the measurements being made by different methods across sampling plots and counties. The effects of different methods on the SPD were considered significant if P < 0.05. Through these three verification methods, the model was shown to roughly match the in situ decrease in soil P on the Tibetan Plateau over the last three decades.

Predicted trends in the soil P budget and their associated change points for the time period 2021–2100 were analyzed using the Mann–Kendall and change-point tests using the trend92 and cpm93 packages, respectively, in R version 4.2.1. The Kolmogorov–Smirnov test was used to extract the change points because the data did not conform to a Gaussian distribution.

Data availability
All data are available in the article and its Supplementary Information. Source data are provided with this paper. The daily rainfall, daily temperature and daily wind speed are available at https://data.cma.cn/; the Daily snow depth is available at https://data.tpdc.ac.cn/home; the NDVI is available at http://www.csdata.org/; the N deposition is available at https://doi.org/10.3334/ORNLDAAC/1220; the ANPP is available at http://www.glass.umd.edu/Download.html; the DEM is available at https://www.jpl.nasa.gov/missions/shuttle-radar-topography-mission-srtm; the LUCC (land-use and land-cover change) is available at https://www.earth-system-science-data.net/; the soil property data (e.g., sand, silt, clay and SOC) are available at https://www.soilgrids.org/; the CaCO3 data are available at https://gaez.fao.org/pages/hwsd; the CNRDv1.0 dataset is available at https://data.tpdc.ac.cn/zh-hans/data/8b6a12c7-c8f9-465a-b449-852fbff51853/; and the future population data are available at https://geography.nuist.edu.cn/2019/1113/c1954a147560/pagem.psp. The SPD data and the geographical/environmental characteristics for the resampling plots are available in Supplementary Data 3.

Code availability
Figures were created in Microsoft PowerPoint software, R 4.2.1, MATLAB R2022a, Adobe Acrobat DC, and ArcGIS 10.6 (https://www.esri.com/). The data analysis codes have been deposited in the figshare database at https://figshare.com/articles/code/29134574.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "With the deepening of the Second Tibetan Plateau Scientific Expedition and Research Program, research on the ecological problems of the Tibetan Plateau has become a hotspot in China, and thousands of researchers have been involved in this massive endeavor. The reviewer has also had the honor of participating in this work and has gained an understanding of the ecological problems of the Tibetan Plateau. This manuscript collects a huge amount of data and analyzes the soil phosphorus crisis on the Tibetan Plateau. It also predicts the evolution of soil phosphorus under different climate scenarios, which is of practical significance for a comprehensive understanding of soil phosphorus changes on the Plateau. For the following reasons, I have taken a negative view of the publication of the manuscript and recommend that it be resubmitted to other Journals specializing in Soil Science. First, the ecological problems facing the Tibetan Plateau in the context of global change are many and varied, including glacier retreat, permafrost thinning, lake expansion, reduced water resources, grassland degradation, declining biodiversity, and increased risk of natural disasters (e.g., glacial lake outbursts); even soil nutrient depletion is considered to be a serious problem, with soil nitrogen and organic matter depletion being the most pressing, and soil phosphorus not yet high on the list. Second, it is believed that for soil phosphorus stock, the study of soil active phosphorus is more meaningful and important, yet the authors seem to neglect this. Third, the omission of farmland land is the shortcoming of this study: 1) it is precisely because of the complexity of farmland P that it needs to be studied; 2) it should not be ignored just because of its small size; farmland is generally located in gullies and riverways, precisely the areas where soil phosphorus converges, and carrying out an in-depth study of it would be more meaningful and valuable. Some specific comments are as follows:\nLine 31-33: Soil phosphorus in this manuscript was calculated at a depth of 0-30 cm, and the data cited here are from 0-100 cm. For consistency of data, the authors should have calculated the proportion of phosphorus accounted for in the 0-30 cm soil layer, rather than 0-100 cm.\nLine 41-43: Here should be a clearer explanation of how soil phosphorus affects soil carbon.\nLine 91-98: Reliability of data? Comparison with other studies is recommended to verify the reliability of the results. There are more reports on changes in soil nutrients (including nitrogen, phosphorus and potassium) on the Tibetan Plateau, and a 33 percent decline from 1980s to 2020s is a huge value.\nLine 197-199: The reason for this should be explained here.\nLine 229-233: These control measures were not derived from this study, but are generally accepted as common sense. In addition, most areas of the plateau are not suitable for afforestation, particularly the north-western region.", "reviewer_2": "The authors employed modified process-balanced model to evaluate the budget of phosphorus in Tibetan Plateau during 1980 and 2020 and predicted the potential P changes by the end of this century. Overall, this study is interesting and contributes to the changes in soil P regarding climate change. However, from the perspective of modeling, there were several issues should be clarified further.\nMajor comments:\n1. The long-term repeated sampling of soils is the basis for this study, and that long-term repeated sampling is a magnificent achievement and the source of much new understanding of soil dynamics. Unfortunately, there are hidden uncertainties behind the apparent precision with which results are stated. To the best of my knowledge, the accurate latitudinal and longitudinal information of sampling sites in the 1980s was lacked. How to account for the uncertainty involved in the analysis of paired samples? Is the size of collected soil dataset enough to quantify the spatio-temporal changes in soil P (n=211)?\n2. One of the most important conclusions is the role of soil erosion in soil P dynamics (81%). Regarding the changing magnitude of soil P, the hypothesis of consistent changing rates of soil P concentrations should be reject (Line 410-413 in the Supplementary Material). The changes in soil P might be nonlinear, and the changing rates are closely correlated with the initial P pool.\n3. I am not sure whether it is credible to estimate the soil P stock based on a legacy soil type map (1:1,000,000). How to quantify the prediction uncertainty of soil P stock estimation? In addition, the soil types and vegetation types might be changed during 40 years. The advanced methods, such as machine learning, might yield robust performance.\n4. The sensitivity analysis of the process-based models should be conducted to quantify the uncertainty contributions of the various parameters to the models.\n5. If soil erosion explains most of the P depletion, it must be pointed out that during the last several decades overall soil erosion become less severe as vegetation covergae as a main control of soil erosion has become significantly better.\n6. The introduction section should summarize the current understanding of soil P changes in the Tibetan Plateau. It will help to identify the key knowledge gaps that are addressed in this study.\nMinor comments:\nL22: One decimal place is enough.\nL51-53: I don’t think it’s a drawback of soil sampling. We cannot demand other irrelevant information from soil measurements.\nL78: Four decades?\nL88-90: It is better to briefly explain the reason for this spatial pattern.\nL118: “The proven MPFE model suggested that the…” - replaced by “The…”\nL163-165: Is it possible to model or incorporate the information on eating preference and population in the future? The livestock products act as important P outflow."}, "author_response": "Response to Reviewers\nFirst of all, we would like to thank the reviewers for their constructive comments and suggestions, which we believe have substantially strengthened our manuscript. We have considered all the suggestions and comments very seriously and have adjusted our manuscript accordingly. Below, we provide a detailed point-by-point response to the reviews. To clarify how and where we addressed the suggestions, we have included our responses directly after the comments and suggestions. The actual changes (in red) are indicated in the revised manuscript and revised supplementary information, with their corresponding page and line numbers.\n[Comment 1] With the deepening of the Second Tibetan Plateau Scientific Expedition and Research Program, research on the ecological problems of the Tibetan Plateau has become a hotspot in China, and thousands of researchers have been involved in this massive endeavor. The reviewer has also had the honor of participating in this work and has gained an understanding of the ecological problems of the Tibetan Plateau. This manuscript collects a huge amount of data and analyzes the soil phosphorus crisis on the Tibetan Plateau. It also predicts the evolution of soil phosphorus under different climate scenarios, which is of practical significance for a comprehensive understanding of soil phosphorus changes on the Plateau. For the following reasons, I have taken a negative view of the publication of the manuscript and recommend that it be resubmitted to other Journals specializing in Soil Science.\n[Response 1] We thank the reviewer for the positive evaluation of the manuscript and the constructive criticism. The Tibetan Plateau, mainly located in southwestern China, is referred to as the world's Third Pole. It covers an area of about 2.5 million km2, with an average altitude of 4,000 m (Zhang, et al., 2019). The Tibetan Plateau is an important ecological safety barrier and strategic resource reserve base (Wang, et al., 2017). The Second Tibetan Plateau Scientific Expedition and Research was launched and implemented in 2017. The leader of the expedition team was Yao Tandong, an academician of the Chinese Academy of Sciences. It involved more than 28,000 participants in conducting scientific research covering the Tibetan Plateau (Times, 2024). The author is also honored to participate in this research. Through extensive literature data collection and field sampling, we obtained a long-term soil P dataset, which is of great significance to support the sustainable use and future management of P in permafrost regions.\nWe revised the three main issues of your concern to the best of our ability, including the importance of P in the permafrost region, soil AP, and soil P in the cropland ecosystem. For details, please see the following Response 2-Response 4.\n1）The importance of P in permafrost region\nWe have rephrased the relevant sentence and added more information on the\n2\nimportance of soil P in permafrost regions in the Introduction. Please see [Response 2].\n2）Soil AP data\nWe have made additional determinations of soil AP and added the supplementary results and discussion of soil AP in the revised version. Please see [Response 3].\n3）Soil P in cropland ecosystem\nWe have supplemented the estimation of soil P change in the cropland ecosystem. Firstly, we already have some resampled soil samples from 11 cropland plots, and the soil P has been determined. Secondly, we have conducted supplementary field sampling during August and September of 2024 in the hope that the estimation of soil P change on the Tibetan Plateau could be more accurate. Thirdly, P inflow and outflow in the cropland ecosystem were incorporated into the MPEE model. Please see [Response 4].\nBy addressing your and Reviewer 2’s comments, we feel that the revised Manuscript (MS) has been greatly improved. Our results are of great importance for the sustainable use of nutrients and the implementing nutrient management practices in permafrost areas. Hence, they will be of great interest to the readership of Nature Communications from multiple perspectives:\n1) Research topic:\na. As noted by the reviewer, this study focuses on P budgets across the permafrost regions, instead of the carbon (C) and nitrogen (N) budgets. While some papers on C and N budgets across the permafrost regions were published (Ding, et al., 2017; Kou, et al., 2020; Chen, et al., 2022; Huang, 2022), no research has yet conducted the analysis for P.\nb. Quantifying the temporal changes of the P stored in permafrost is crucial and urgently needed in guiding strategies to tackle the grand challenge of future P limitation (Alewell, et al., 2020; Barbieri, et al., 2022), as well as the environmental crisis caused by huge losses of P into water bodies that pushes the Earth's environment to its boundary, beyond which come abrupt environmental changes with catastrophic impacts (Chen and Graedel, 2016; Brownlie, et al., 2022).\n2) Methods and dataset:\na. We used a unique resampling method to estimate the dynamics of soil P pools in\n3\npermafrost regions over the last three decades.\nb. This is the first use of the dynamic background of soil P data for the two time periods to verify and improve the accuracy of the P flow in ecosystems. It is of great significance to reduce the uncertainty of estimating the past biogeochemical cycling of P as well as future predictions of nutrient resource management.\nc. As the reviewer has mentioned, this manuscript collects a huge amount of data and analyzes the soil P crisis on the Tibetan Plateau. It also predicts the evolution of soil P under different climate scenarios, which is of practical significance for a comprehensive understanding of soil P changes on the plateau.\n3) Findings:\na. An astonishing depletion in soil P stocks over a short period of time has not previously been detected. The average net soil P loss of the Tibetan alpine permafrost region was about 3.5 Tg P y−1, which is roughly four-fifths of the annual application of chemical P fertilizer to cropland in China and one-fifth of the annual application of chemical P fertilizer globally.\nb. We propose that the rapid loss of soil P would increase the risk of P limitation in permafrost regions in the future. Large amounts of soil P are transferred from the permafrost to the adjacent water body, contributing to the potential high risk of decreasing biodiversity and the production of algal blooms, threatening the health of the environment as well as affecting the trajectory of ecosystem C fluxes.\nc. We predict that soil P will continue to decrease rapidly from the 2020s to the 2050s. The tipping point for the P budget is around 2050. Soil and water conservation needs to be considered as soon as possible to maintain soil health and the nutrient status of soils as well as decrease the ecologically negative impacts on interior freshwater lakes and ocean waters.\n[Comment 2] First, the ecological problems facing the Tibetan Plateau in the context of global change are many and varied, including glacier retreat, permafrost thinning, lake expansion, reduced water resources, grassland degradation, declining biodiversity, and increased risk of natural disasters (e.g., glacial lake outbursts); even soil nutrient depletion is considered to be a serious problem, with soil nitrogen and organic matter depletion being the most pressing, and soil phosphorus not yet high on the list.\n[Response 2] Sorry for the confusion. Over the past 40 years, the climate and environment in this region have experienced dramatic changes, including accelerated glacier retreat and imbalance of the Asian water tower (Yao, et al., 2022), permafrost thinning (Yang, et al., 2010), lake expansion (Xu, et al., 2024), grassland degradation (Wang, et al., 2022a), increased risk of natural disasters (e.g., glacial lake outbursts) (Nie, et al., 2021), and C and N dynamics (Kou, et al., 2020; Chen, et al., 2022). These transformations have affected the structure and function of the regional ecosystem, with significant implications for human survival and development.\nQuantifying the temporal changes of the P stored in permafrost is crucial and urgently needed in guiding strategies to tackle the grand challenge of future P limitation (Alewell, et al., 2020; Barbieri, et al., 2022), as well as the environmental crisis caused by huge losses of P into water bodies that pushes the Earth's environment to its boundary, beyond which come abrupt environmental changes with catastrophic impacts (Chen and Graedel, 2016; Brownlie, et al., 2022).\nSpecifically, P is a key element in organic molecules (e.g., nucleic acids, ATP and phospholipids) and has an important role in regulating the growth, development and reproduction of living organisms (Sattari, et al., 2012; Sattari, et al., 2016; Brownlie, et al., 2021). Permafrost regions occupy about 24% of the land surface in the northern hemisphere (Zhang, et al., 2008), but their soil P reserves account for >30% of the northern hemisphere total (He, et al., 2021). Permafrost soils contain a huge reservoir of P as a result of the rejuvenation of soils, weak leaching processes, and the low uptake of nutrients by plants, and can therefore be regarded as important reserves of P in terrestrial ecosystems (He, et al., 2021).\nThe recent rapid changes in the Earth’s climate are simultaneously causing the widespread degradation of permafrost (Slater and Lawrence, 2013) and altering the hydrology of high-latitude/high-altitude regions (Peterson, et al., 2006), exposing mineral nutrients previously protected in permafrost to transport and processing in lakes, rivers and estuaries (Abbott, et al., 2015). The unidirectional flow of P from terrestrial ecosystems to freshwater bodies and, eventually, the oceans is considered to be the greatest sustainability challenge in the P cycle. The redistribution of P on spatiotemporal scales could regulate the composition of local, downstream and even coastal plant communities, the net primary productivity and the health of the environment (Brownlie, et al., 2021). In addition, the availability of P affects the trajectory of the permafrost C cycle by influencing gross primary productivity and net ecosystem productivity (Yang, et al., 2021), and the associated fluxes of C in aquatic sediments (Brownlie, et al., 2021) by influencing soil microbial processes and greenhouse gas emissions. This can have a positive or negative feedback role in global climate change.\nThe P present in soils is primarily derived from the slow processes of rock weathering, but can be rapidly replenished by fertilizers in human-managed soil systems (Bouwman, et al., 2013). The unidirectional flow of P from terrestrial ecosystems (e.g., soils, mineral reserves, and agricultural products) to freshwaters and eventually the oceans is considered to be the greatest sustainability challenge in the P cycle (Alewell, et al., 2020). The increasing consumption of food produced by human agriculture and the longer growing seasons caused by the current global climate crisis have increased the P requirements of terrestrial ecosystems (Sattari, et al., 2012; Liu, et al., 2016). It has been suggested that almost half of the global AP will be depleted by the end of this century (Van Vuuren, et al., 2010). In the Introduction, we have added more information on the importance of soil P in permafrost regions. Please see Revised Manuscript page 2, line 40-42; page 3, line 48-52, line 54-57, line 66-71.\n[Comment 3] Second, it is believed that for soil phosphorus stock, the study of soil active phosphorus is more meaningful and important, yet the authors seem to neglect this.\n[Response 3] Very good comments. Fortunately, the soil samples are well preserved. According to your comment, we have made additional determination of soil AP, which could provide more information to the reader. We have added the result and discussion of soil AP in the revised version. Please see Revised Manuscript page 8-9, line 208-210; Supplementary Text 4 (Revised Supplementary Information page 27, line 668-686); and Supplementary Figs. 22 and 23.\nThe supplementary result and discussion are as follows:\nDiscussion\n“Arctic and alpine permafrost ecosystems are generally considered to be nitrogen (N) limited, and the rapid loss of soil P could increase the risk of P limitation in permafrost regions. Particularly, soil AP in the grassland ecosystem (occupying the largest area of the Tibetan Plateau) also decreased significantly, indicating that the soil P supply became insufficient (Supplementary Text 4). During the last few decades, an enhanced supply of available N in soils derived from the release of originally frozen N or accelerated N mineralization has been observed with climate warming or the thawing of permafrost.”\nSupplementary text 4\n“In general, there was no significant change in soil AP density across the Tibetan Plateau during the two sampling periods (Supplementary Fig. 22). However, soil AP density showed different variation trends for different ecosystems. Grassland showed an overall decreasing pattern of AP dynamics, while the AP density increased significantly in forest and cropland (P < 0.05; Supplementary Fig. 22). The soil AP dynamics were linked to the differences in climate, vegetation biomass, soil texture, and agronomic management practices. Soil AP dynamics over the entire Tibetan grassland were positively related to mean annual temperature (MAT), normalized difference vegetation index (NDVI), nitrogen (N) deposition, and sand content, and negatively correlated with soil clay content (P < 0.05; Supplementary Fig. 23); whereas, no correlation was established between changes in AP and mean annual precipitation (MAP). Firstly, the high temperatures and increased biomass/litter may exist in response to the elevated AP densities associated with accelerated plant growth and organic matter decomposition rates. Secondly, a higher N deposition rate can stimulate phosphatase activities and thus increase the inorganic P supply, which could indirectly promote the bioavailability of soil P. Thirdly, clay particles generally have the smallest diameter and largest specific surface area. However, the clay particles in alpine soil are negatively charged, which will repel the also negatively charged AP. Lastly, a certain amount of chemical fertilizer and organic fertilizer could directly enhance the soil AP densities in the cropland ecosystem.”.\nUnlike the soil total P, the soil AP varies greatly in different months and seasons in permafrost areas due to the strong freeze-thaw effect. The observed investigation shows that the range of variation of soil AP in permafrost areas can reach more than 200 times within a year (Zhou, et al., 2020; Huang, et al., 2021). This phenomenon might be caused by continuous freeze-thaw cycles affecting the complex adsorption and desorption processes of soil P. For example, previous studies have proven that freeze-thaw cycles significantly decrease the soil adsorption capacity of P and increase the soil AP (Fan, et al., 2014), which could also cause P depletion associated with runoff and erosion (Yakutina, 2011). The rapid changes in AP of permafrost soil make it difficult for MPEF to incorporate it because it requires extensive field sampling and laboratory analysis to reveal the dynamic changes of soil AP in different months. Therefore, we have added the Future improvements in the revised version. Please see Revised Manuscript page 11, line 269-272.\nThe Future improvements is as follows: “The analysis presented in this paper is a systematic attempt to assess the historical changes in soil P pools and the key flows of P in permafrost regions using a unique combination of a soil resampling method and a modified process-balanced model, but both the methodology and spatial resolution of the MPFE model need additional work. For instance, the soil AP varies greatly in different months and seasons in permafrost areas due to the strong freeze–thaw effect. Therefore, high-temporal-resolution soil AP data need to be incorporated into the MPFE model in the future”. Thank you for your understanding!\n[Comment 4] Third, the omission of farmland land is the shortcoming of this study: 1) it is precisely because of the complexity of farmland P that it needs to be studied; 2) it should not be ignored just because of its small size; farmland is generally located in and riverways, precisely the areas where soil phosphorus converges, and carrying out an in-depth study of it would be more meaningful and valuable.\n[Response 4] Very good comment. Regarding the soil P change in farmland ecosystems, we are sorry that our original manuscript did not consider it comprehensively. Actually, P is an important nutrient for maintaining crop yields in agroecosystems. Soil P dynamics across agricultural systems may reflect contrasting agricultural practices and broader agricultural policies. The export of P due to harvesting tends to increase soil P limitation. Although the area of arable soil is relatively small across the Tibetan Plateau, the failure to include cropland in the original manuscript could lead to biases in the calculation of soil P stocks.\nCropland is mainly located in gullies and valleys of Huangshui, Yarlung Tsangpo River, Lancang River, Niyang River, Nianchu River, Yalong River and Shiquan River (Duan, et al., 2024). As the important ecological protection areas and the key agricultural production areas across the Tibetan Plateau, they provide important guarantees for local grain self-sufficiency (Wang, et al., 2024). Following your suggestions, we have supplemented the estimation of soil P change in arable soils. Firstly, we have already resampled soil samples from 11 cropland plots, and the soil P has been determined. Secondly, we have conducted supplementary field sampling during August and September of 2024 in the hope that the estimation of soil P change in the Tibetan Plateau would be more accurate. An additional 33 cropland plots were collected, bringing the total of 44 cropland plots included in the revised version (Fig. R1). Cropland management practices, including fertilizer application and crop yield, were also investigated. The field sampling method and sample analysis method were the same as our previous methods. We comprehensively compiled and summarized soil P change and inflows/outflows in the cropland ecosystem. Thirdly, we analyzed the historical changes of P concentration and stock in cropland soils, and then incorporated P data from cropland ecosystems into the MPFE model.\nPlease see the Figs. 1-5, Results and Discussion (Revised Manuscript page 4, line 94; page 6, line 130, 136 and 144), Methods (Revised Manuscript page 12, line 295-300, line 303; page 13, line 331; page 14, line 377; page 15, line 393) and Supplementary Text 2 (Revised Supplementary Information page 5, line 115-120; page 6, line 131-134; page 8, line 208-210; page 9-10, line 234-248; page 11, line 271-277; page 12-13, line 316-323; page 16, line 421-423; page 17, line 438-443).\nFig. R1 Spatial distribution of supplementary resampling plots for cropland on the Tibetan Plateau\nSome specific comments are as follows:\n[Comment 5] Line 31-33: Soil phosphorus in this manuscript was calculated at a depth of 0-30 cm, and the data cited here are from 0-100 cm. For consistency of data, the authors should have calculated the proportion of phosphorus accounted for in the 0-30 cm soil layer, rather than 0-100 cm.\n[Response 5] This is a good catch. We have updated the data on soil P stock for the depth of 0-30 cm. Please see Revised Manuscript page 2, line 32-33.\n[Comment 6] Line 41-43: Here should be a clearer explanation of how soil phosphorus affects soil carbon.\n[Response 6] Thank you for the careful review. We have rewritten this sentence and added clear information here. Please see Revised Manuscript page 2-3, line 44-46.\n[Comment 7] Line 91-98: Reliability of data? Comparison with other studies is recommended to verify the reliability of the results. There are more reports on changes in soil nutrients (including nitrogen, phosphorus and potassium) on the Tibetan Plateau, and a 33 percent decline from 1980s to 2020s is a huge value.\n[Response 7] Thank you very much for the careful review. In order to verify the reliability of the results, we compared the direction and magnitude of soil P changes with previous studies. We found an overall significant decrease in soil P across the Tibetan Plateau, which concurred with national and regional soil investigations (non-resampling method) and machine learning techniques reporting a significant decline in topsoil P (Yang, et al., 2014) and AP (Tian, et al., 2019) stocks from the 1980s to the 2010s. Due to the lack of initial soil P stock data using the resampling method, it is not possible to compare the amount of soil P loss in the present study with those of other studies. However, the rapid decrease in soil P (14.9 kg ha−1 yr−1) across the Tibetan Plateau was within the range of the global soil P depletion rate (4–19 kg ha−1 yr−1) (Alewell, et al., 2020). Therefore, in order to enhance the credibility of the results, we strengthened the comparative analysis between the results of this study and those of previous studies. Please see Revised Manuscript page 4, line 95-96.\nSince the largest part of soil P loss is caused by water erosion, we also compared the data of P loss caused by water erosion on the Tibetan Plateau with other studies. Our area-related calculations resulted in an average outflow of P from the alpine permafrost region caused by water erosion of about 16.1 kg P ha−1 yr−1, which is slightly higher than that of previous calculations worldwide (Supplementary Table 6). The main reason could be that the annual precipitation is concentrated in the summer months and the complex conditions, including differences in elevation, low vegetation coverage, and weak resistance to erosion, all led to strong water erosion in this alpine permafrost region. Please see Revised Manuscript page 7, line 157-160 and Revised Supplementary Information Supplementary Table 6.\n[Comment 8] Line 197-199: The reason for this should be explained here.\n[Response 8] Thanks for this point. We have added the reason here. Please see Revised Manuscript page 9, line 217-218.\n[Comment 9] Line 229-233: These control measures were not derived from this study, but are generally accepted as common sense. In addition, most areas of the plateau are not suitable for afforestation, particularly the north-western region.\n[Response 9] Thanks for this suggestion. After careful consideration, we have deleted this sentence. Please see Revised Manuscript page 10, line 250.\n[Comment 1] The long-term repeated sampling of soils is the basis for this study, and that long-term repeated sampling is a magnificent achievement and the source of much new understanding of soil dynamics.\n[Response 1] Thanks for the reviewer’s positive comment.\n[Comment 2] Unfortunately, there are hidden uncertainties behind the apparent precision with which results are stated. To the best of my knowledge, the accurate latitudinal and longitudinal information of sampling plots in the 1980s was lacked. How to account for the uncertainty involved in the analysis of paired samples? Is the size of collected soil dataset enough to quantify the spatio-temporal changes in soil P (n=211)?\n[Response 2] Very good comment! Following your comments, we have made the following two major revisions:\n⚫ Indeed, as you said, how to accurately locate sampling points is an important issue in this study. To improve the precision of positioning the resampling plots, we added more detailed information for the accurate identification of resampling plots. Please see Supplementary Text 1 (Revised Supplementary Information page 2-3, line 23-55).\nFirstly, the resampling at the plot scale was determined based on geographical information for the historical sampling and digitalized historical maps, including administrative divisions, types of vegetation, records of landmarks, the dominant plant species, the soil type, and the parent material. A few counties, towns, and villages have changed their place names. To further verify the exact locations of the sampling plots, we drove to the target provinces, counties, towns, and villages where the plots were located and carried out on-plot verification. We contacted local government staff to help us determine the previous and present place names and to guide the way to the sampling plots.\nSecondly, the soil genetic horizon was an important reference. We took three soil cores in each plot to generate descriptions of their vertical properties (e.g., layering, color, and texture for each layer) for comparisons with those from the 1980s. In general, the soil genetic horizon of the Tibetan Plateau is relatively stable owing to the low amounts of human disturbance, and thus the consistency between each layer of the profile in the 2020s and the historical record was the key piece of evidence used to determine whether the plot was selected correctly.\nFigure R2. Layers of soil profile\nThirdly, to further assess the accuracy of our plots, we compared all the locations with other independent references, such as secondary use in different maps, vegetation type maps, soil type maps, topographic maps (elevation, slope, and terrain), land use maps (residential areas, rivers, and lakes), road maps, and administrative boundary maps that were developed before 1990s (scaled at 1:50,000–1:1000,000). We overlaid different maps to obtain a range that was likely to include the sampling locations from the 1980s campaign. The distance between 91% of the relocated plots (233 plots) and of those in the 1980s was <100 m (see details in Supplementary Fig. 19). For 22 plots, the low-resolution independent references failed to capture geographic information in sufficient detail because those plots were distributed in vegetation transition zones or intrazonal soil areas. We further checked photos of the vegetation type, topographic information, soil type, and genetic horizon of the 22 resampling plots from the 2020s campaign, which were consistent with historical written records.\nFourthly, we also consulted several staff who participated in the 1980s campaign to determine the exact locations of the sampling plots (Figure R3). We apologize for missing this important information in the previous version. We have added this additional information on the exact location of the resample plots in the revised versio\nn.\n[Redacted]\nFigure R3. Communication with sampling staff who participated the 1980s campaign\nLastly, in a few repeated sampling plots, a relatively high change in soil P concentration between the two periods most likely arises from the lack of sufficiently accurate latitudinal and longitudinal geographic information on legacy pits, together with the high spatial variability in soil P. To avoid confusion, we replace the “sampling sites/pits” with “sampling plots” in the revised version.\n⚫ Regarding the number of sampling plots, we do agree that the more sampling plots, the more accurate the calculation of soil P stock would be. Following the reviewer 1’s suggestions, we have added the soil P data from cropland in the revised version. We supplemented soil P data from 44 cropland plots. Although the number of plots may still be insufficient, we've done our best to get as many resampling plots as possible for such a high-altitude and harsh natural environment. To quantify soil P changes on regional scales, more long-term soil sampling is needed, but it requires significant human resources and financial support.\nThe Tibetan Plateau covers the area (25.99–39.82°N, 73.50–104.67°E), and has a broad environmental gradient. The Tibetan permafrost accounts for 75% of the area of alpine permafrost in the northern hemisphere (Yang, et al., 2010). The Tibetan Plateau has experienced significant increases in the mean annual temperature of 0.4 °C decade–1 (Pepin, et al., 2015), and precipitation has also shown a general upward trend during the last three decades (Wang, et al., 2022b). This area serves as an ideal platform for exploring the dynamics of soil P over a broad geographic scale. There is great spatial variability in climatic conditions. The dominant ecosystems on the plateau are alpine grassland and forest, which shift from alpine forest in the southeastern area to alpine grassland in the northwestern area. Cropland is mainly located in gullies and valleys of Huangshui, Yarlung Zangbo River, Lancang River, Niyang River, Nianchu River, Yalong River, and Shiquan River (Duan, et al., 2024). Those 255 sample plots cover the major ecosystem types (forest, grassland, and cropland) and soil orders (humic cambisols, haplic xerosols, orthic luvisol and dystric podzoluvisols, and ect.) in this region according to the World Reference Base for Soil Resources (Lu, et al., 2004), which can provide valuable basic data for the calculation of soil P change. These plots were set throughout the geographical extent of alpine forests, grasslands and cropland, covering a wide range of climates (mean annual temperature (MAT): -12.6–14.3 °C; mean annual precipitation (MAP): 22.4–1299.5 mm); vegetation characteristic (net primary productivity (NPP): 6.1–1173.7 g C m-2 yr-1; normalized difference vegetation index (NDVI): 0–0.8), soil properties (soil total phosphorus (P): 0.2–1.6 g kg-1; soil available phosphorus (AP): 0.4–97.6 mg kg-1) across the plateau (Fig. R4). Thus, the sampling plots are typical and representative, and we have added this background information in the revised version. Please see Revised Manuscript page 13, line 330-336.\n[Comment 3] One of the most important conclusions is the role of soil erosion in soil P dynamics (81%). Regarding the changing magnitude of soil P, the hypothesis of consistent changing rates of soil P concentrations should be reject (Line 410-413 in the Supplementary Material). The changes in soil P might be nonlinear, and the changing rates are closely correlated with the initial P pool.\n[Response 3] We thank the reviewer for the insightful comment that the change rate of soil P concentration might be nonlinear. If the annual soil P concentration for each plot is available, it can provide precise data for estimating soil P stock changes and improving the accuracy of the MPFE model. To be honest, annual survey sampling for such a large area is very difficult; it requires significant human resources and financial support. Since we do not have continuous soil P data for each plot, how to interpolate the data of the missing years is a major problem in the current model simulation.\nInspired by your comment, we reanalyzed the data and found that the change rate of soil P was not consistent during the past three decades. As you guessed, the changing rates of soil P are closely correlated with the initial soil P concentration, and show an exponential pattern of change (Supplementary Fig. 21). Specifically, soil with a high initial P concentration tends to lose P, and the higher the initial value, the higher the loss rate. However, soil with a low initial P concentration tends to accumulate P, and the lower the initial value, the higher the accumulation rate. Therefore, we used a space-for-time approach to fit the functional relationship between the change rate of soil P and the initial P concentration for the increasing and decreasing sites, respectively. We have replaced the consistent changing rates with exponential changing rates for interpolating the soil P concentration, and then reanalyzed relevant P flows in the revised version. Please see Revised Supplementary Information page 19, line 495-501; page 20, line 532-535.\nSpace-for-time substitution is a common method to investigate long-term phenomena in soil, geography, and ecology (Pickett, 1989; Huang, et al., 2019a). It has been proved that this approach is successful in predicting soil nutrients with successional dynamics (Hartley, et al., 2021; Wang, et al., 2022a). Although we use a space-for-time approach to make the predicted soil P concentration consistent with the actual trend of soil P concentration as much as possible, it may have a certain gap with the actual soil P concentration for a given year. Further improvements in the accuracy of this study can be achieved through the acquisition of more precise soil P concentration data. Due to this point, together with your comments, we added the potential uncertainty of changes in soil type on soil P stock changes in the Limitation section of the revised version. Please see Revised Manuscript Page 10, line 260-265.\nThe Limitation is as follows: “The annual soil P concentration can provide precise data for estimating the soil P stock and improving the accuracy of the MPFE model. However, annual survey sampling for such a large area is very difficult because of the huge challenges in terms of labor and experimental costs. We used a space-for-time substitution method to interpolate the soil P concentration of missing years from the 1980s to 2020s, which may translate to discrepancies with the actual soil P concentration for a given year.”.\n[Comment 4] I am not sure whether it is credible to estimate the soil P stock based on a legacy soil type map (1:1,000,000). How to quantify the prediction uncertainty of soil P stock estimation? In addition, the soil types and vegetation types might be changed during 40 years. The advanced methods, such as machine learning, might yield robust performance.\n[Response 4] Very good comment. We have summarized your comments into two main questions. ①Can we use advanced methods (i.e., machine learning) to calculate soil P stocks? ②If soil type and vegetation type maps are used to calculate soil P stocks, how can we reduce the prediction uncertainty, including the resolution (1:1,000,000), and the potential historical changes in vegetation type and soil type?\nFirstly, following your suggestions, we have tried to use machine learning for estimating the soil P stock to reduce the uncertainty in the revised version. Widely used digital soil mapping techniques were adopted to model the spatiotemporal pattern of elements. It was assumed that soil nutrients can be modeled as a function of soil-forming factors, such as landform variables, climate, lithology, human activity, vegetation characteristics, and physical and chemical properties of soil (Table R1).\nCategory Name Type Unit Data source Spatial Time average Reference\nBulk density Static g cm-3 SoilGrid 1000m current https://www.isric.org/explore/soilgrids\npH (H O) Static - SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nOrganic carbon density Static hg m-3 SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nTotal nitrogen Static % SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nSand content Static % SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nSilt content Static % SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nSoil\nClay content Static % SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nCoarse fragments\nStatic % SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nvolumetric total\nCation exchange\nStatic mmol kg-1 SoilGrid 1000m current https://www.isric.org/explore/soilgrids\ncapacity at ph7\nVolumetric Water\nStatic ‰ SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nContent at 10kPa\nStatic ‰ SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nContent at 10kPa\nVolumetric Water\nStatic ‰ SoilGrid 1000m current https://www.isric.org/explore/soilgrids\nContent at 1500kPa\nSoil Respiration Static g m-2 yr-1 Earth Data 1000m 1961-2016 https://daac.ornl.gov/\nSoil autotrophic\nStatic g m-2 yr-1 Earth Data 1000m 1961-2016 https://daac.ornl.gov/\nrespiration\nSoil heterotrophic\nStatic g m-2 yr-1 Earth Data 1000m 1961-2016 https://daac.ornl.gov/\nrespiration\nhttps://figshare.com/articles/dataset/Glob\nal_Aridity_Index_and_Potential_Evapotr\nAridity Index Static - ai_et0 1000m 1970-2000\nanspiration_ET0_Climate_Database_v2/\n7504448/3\nClimate\nhttps://figshare.com/articles/dataset/Glob\nal_Aridity_Index_and_Potential_Evapotr\nEvapotranspiration Static - ai_et0 1000m 1970-2000\nanspiration_ET0_Climate_Database_v2/\n7504448/3\nMinimum temperature Static ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nMaximum temperature Static ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nPrecipitation Static mm WorldClim 1000m 1970-2000 https://worldclim.org/\nSolar radiation Static kJ m-2 day-1 WorldClim 1000m 1970-2000 https://worldclim.org/\nWind speed Static m s-1 WorldClim 1000m 1970-2000 https://worldclim.org/\nWater vapor pressure Static kPa WorldClim 1000m 1970-2000 https://worldclim.org/\nMean Diurnal Range\n(Mean of monthly (max Static ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\ntemp - min temp))\nIsothermality\n(TempDiurnalRange/Te\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nmp AnnualRange)\n(×100)\nTemperature Seasonality\n(standard deviation Static ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\n×100)\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nWarmest Month\nMin Temperature of\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nColdest Month\nTemperature Annual\nRange\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\n(TempMaxWarmestM-\nTempMinColdestM)\nMean Temperature of\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nWettest Quarter\nMean Temperature of\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nDriest Quarter\nMean Temperature of\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nWarmest Quarter\nMean Temperature of\nStatic ℃ WorldClim 1000m 1970-2000 https://worldclim.org/\nColdest Quarter\nPrecipitation of Wettest Static mm WorldClim 1000m 1970-2000 https://worldclim.org/\nPrecipitation of Driest\nStatic mm WorldClim 1000m 1970-2000 https://worldclim.org/\nMonth\nPrecipitation\nSeasonality (Coefficient Static mm WorldClim 1000m 1970-2000 https://worldclim.org/\nof Variation)\nPrecipitation of Wettest\nStatic mm WorldClim 1000m 1970-2000 https://worldclim.org/\nQuarter\nPrecipitation of Driest\nStatic mm WorldClim 1000m 1970-2000 https://worldclim.org/\nQuarter\nPrecipitation of\nStatic mm WorldClim 1000m 1970-2000 https://worldclim.org/\nWarmest Quarter\nPrecipitation of Coldest\nStatic mm WorldClim 1000m 1970-2000 https://worldclim.org/\nQuarter\nChange of Precipitation Dynamic mm WorldClim 1000m 1985-2020 https://worldclim.org/\nChange of Average\nDynamic ℃ WorldClim 1000m 1985-2020 https://worldclim.org/\ntemperature\nDynamic ℃ WorldClim 1000m 1985-2020 https://worldclim.org/\ntemperature\nChange of Maximum\nDynamic ℃ WorldClim 1000m 1985-2020 https://worldclim.org/\ntemperature\nNitrogen deposition nhx Static g m-2 yr-1 Earth Data 0.5° 1970-2020 https://daac.ornl.gov/\nNitrogen deposition noy Static g m-2 yr-1 Earth Data 0.5° 1970-2020 https://daac.ornl.gov/\nhttps://www.jpl.nasa.gov/missions/shuttl\nElevation Static m SRTM 1000m Current\ne-radar-topography-mission-srtm\nBased on Elevation data, it is calculated\nSlope Static degrees SAGA 1000m Current\nby SAGA software\nBased on Elevation data, it is calculated\nAspect Static - SAGA 1000m Current\nLandform by SAGA software\nBased on Elevation data, it is calculated\nFlow Accumulation Static - SAGA 1000m Current\nby SAGA software\nBased on Elevation data, it is calculated\nLS Factor Static - SAGA 1000m Current\nby SAGA software\nMelton Ruggedness Static - SAGA 1000m Current Based on Elevation data, it is calculated\nMultiresolution index of Based on Elevation data, it is calculated\nStatic - SAGA 1000m Current\nthe ridge top flatness by SAGA software\nMultiresolution index of Based on Elevation data, it is calculated\nStatic - SAGA 1000m Current\nvalley bottom flatness by SAGA software\nBased on Elevation data, it is calculated\nStream Power Index Static - SAGA 1000m Current\nby SAGA software\nTerrain Ruggedness Based on Elevation data, it is calculated\nStatic - SAGA 1000m Current\nIndex by SAGA software\nTerrain surface Based on Elevation data, it is calculated\nStatic - SAGA 1000m Current\nconvexity by SAGA software\nBased on Elevation data, it is calculated\nTerrain surface texture Static - SAGA 1000m Current\nby SAGA software\nTopograohic Wetness Based on Elevation data, it is calculated\nStatic - SAGA 1000m Current\nIndex by SAGA software\nhttps://zenodo.org/record/1464846#.YoX\nLithology Lithology classes Static - USGS 250m Current\nUVOhBw2\nGrazing intensity Dynamic number - 8000m 1985-2020 This study\nactivity\nNDVI Static - CSDATA 8000m 1985-2020 http://www.csdata.org/\nhttp://www.glass.umd.edu/Download.ht\nNPP Static g C m-2 yr-1 Glass 0.05° 1985-2018\nml;\n1980s and\nVegetation Vegetation type Dynamic - TPDC 5000m https://data.tpdc.ac.cn/\n2020s\nChange of NDVI Dynamic - CSDATA 8000m 1985-2020 http://www.csdata.org/\nhttp://www.glass.umd.edu/Download.ht\nChange of NPP Dynamic g C m-2 yr-1 Glass 0.05° 1985-2018\nml;\nWe constructed an indicator system comprising 68 feature variables covering six aspects: soil, climate, landform, lithology, human activity, and vegetation variables. These feature variables were categorized into static and dynamic types to further clarify the impact of background factors and changes on soil P density. Assuming that changes in soil P density can be modeled as functions of factors such as soil, climate, landform, lithology, human activity, and vegetation variables, we adopted a combined method to select independent variables. (1) The correlations between all feature variables and the target variables (P density in the 1980s, P density in the 2020s, and changes in P density from the 1980s to the 2020s) were analyzed; then feature variables with significant correlations were selected (P<0.05). (2) We added the significantly correlated feature variables into Recursive Feature Elimination (RFE) to further reduce the dimensionality of the feature set, thereby enhancing model accuracy and efficiency. (3) We utilized the aforementioned models to perform 30 repeated permutation importance analyses, thereby comprehensively understanding each feature variable’s contribution to model predictions.\nWe employed GridSearchCV to perform hyperparameter tuning for all models, aiming to identify the optimal parameters and ensure reproducibility. However, the optimal tuning processes varied across different models:\nRandom Forest (RF) Model: Tuned parameters included n_estimators (50 to 500), max_depth (0 to 100), min_samples_split (2 to 10), min_samples_leaf (1 to 8), and bootstrap (True and False).\nDecision Tree (DT) Model: Tuned parameters were consistent with the RF model, including max_depth, min_samples_split, and min_samples_leaf. Additionally, max_features values included 1.0, log2, and None; criteria considered were squared error, Friedman MSE, and absolute error.\nSupport Vector Regression (SVR) Model: Features were standardized using StandardScaler to normalize their distributions, which was crucial for SVR performance. Hyperparameters included kernel types (linear, radial basis function, polynomial), regularization parameter C (0.1, 1, 10, 100), epsilon values (0.01, 0.1, 0.2, 0.5), gamma settings (scale, auto), and polynomial degrees (2, 3, 4) for the polynomial kernel.\nArtificial Neural Network (ANN) Model: We used Multi-Layer Perceptron Regressor (MLPRegressor) with a predefined parameter grid, including hidden layer sizes of 50 and 100, activation functions (relu and tanh), solvers (adam and lbfgs), regularization strengths (alpha values of 0.0001, 0.001, 0.01), learning rates (constant and adaptive), and initial learning rates (0.001, 0.0005, 0.0001), to establish a predictive model. Then we configured for a maximum of 10,000 iterations to ensure sufficient training, enabled early stopping by monitoring 10% of validation scores, and set a patience parameter of 2000 iterations to further prevent overfitting.\nBack Propagation Neural Network (BP) Model: We employed a Multi-Layer Perceptron (MLP) Regressor configured for up to 10,000 iterations to allow adequate training, enabled early stopping with 10% of the validation scores, and set the tolerance to 1e-3 to determine convergence criteria. Then we explored various configurations of hidden layer sizes of 50 and 100, activation functions [relu, tanh], regularization parameters alpha values [0.0001, 0.001], and learning rates [constant, adaptive].\neXtreme Gradient Boosting (XGBoost) Model: Tuned parameters included the number of estimators (100, 200, 300, 500), maximum tree depth (3, 5, 7, 10), learning rate (0.01, 0.05, 0.1, 0.2), subsample ratio (0.6, 0.8, 1.0), column sampling by tree (0.6, 0.8, 1.0), gamma values (0, 0.1, 0.2), and minimum child weight (1, 3, 5).\nAll models were tuned using three-fold cross-validation with the objective of minimizing Mean Squared Error (MSE) to identify the optimal hyperparameter sets. After determining the best models, their performance was evaluated using five-fold cross-validation, and the average Root Mean Squared Error (RMSE) was calculated to indicate predictive accuracy. Finally, the models were assessed on an independent test set, obtaining RMSE, pearson correlation coefficients, and R² values to validate their effectiveness. Table R2 showed that the accuracy of machine learning of seven algorithms is relatively low (R2≤0.1).\nTable R2. Prediction accuracy for the seven algorithms. Prediction accuracy was evaluated in terms of the pearson correlation coefficients (R), the coefficient of determination (R2) and root mean square error (RMSE). Multiple Linear Regression (MLR), Random Forest Regression (RF), Decision Tree Regression (DT), Support Vector Regression (SVR), Artificial Neural Network (ANN), Back Propagation Neural Network (BP), and eXtreme Gradient Boosting (XGBoost).\nCross-validate the Post-optimization\nMethods P density\naverage RMSE R R2 RMSE\n1980s 1055.94 -0.12 -0.01 1126.1\n2020s 512.42 -0.14 -0.01 521.68\nMLR\n1980s-\n867.73 -0.07 -0.004 986.70\n2020s\n1980s 975.74 0.20 0.01 1114.82\n2020s 475.94 0.33 0.10 491.57\nRF\n1980s-\n792.52 0.25 0.05 957.60\n2020s\n1980s 1123.77 0.09 -0.30 1279.54\n2020s 510.12 0.32 0.02 513.02\nDT\n1980s-\n923.42 0.07 -0.20 1076.77\n2020s\n1980s 1083.49 -0.12 -0.08 1165.32\n2020s 526.82 -0.14 -0.13 550.25\nSVR\n1980s-\n890.27 -0.07 -0.05 1010.15\n2020s\n1980s 1051.55 -0.12 -0.01 1127.17\n2020s 512.72 -0.14 -0.01 522.32\nANN\n1980s-\n858.37 -0.07 -0.004 986.88\n2020s\n1980s 1894.24 -0.12 -2.03 1952.54\n2020s 1155.86 0.14 -4.26 1189.36\nBP\n1980s-\n1005.62 0.07 -0.29 1116.48\n2020s\n1980s 976.88 0.21 0.03 1103.80\nXGBoost\n2020s 469.99 0.36 0.09 493.89\n30\n1980s-\n837.64 0.27 0.07 947.97\n2020s\nSecondly, since the accuracy of machine learning is relatively low, we still use soil type and vegetation type maps to calculate soil P stock. It is also a common method to calculate soil nutrient stock by multiplying the area of soil types by the corresponding nutrient density both at regional scale (Zhang, et al., 2005; Chapman, et al., 2013; Ding, et al., 2017; Jiang, et al., 2019) and global scale (Yang, et al., 2013). Thus, soil type map is an important basic data to calculate soil storage of C, N, and P. There may be some uncertainties in the estimation accuracy of soil P stock, such as resolution of map, and potential changes in vegetation and soil types over the past decades. We estimated the soil P stock using the following improved methods to reduce the prediction uncertainty in the revised version.\nResolution of soil type map. Firstly, the high-resolution soil type maps were used to calculate soil P stock in the revised version. We used the 1:50,000 (30 m) soil type map to replace 1:100,0000 soil type map for calculating soil P stock. So far, this database is the most detailed and recognized soil type map in China. More than 80,000 scientific researchers participated in the study over 16 years (1979-1995) (Zhang, et al., 2002). Secondly, we also used additional soil type maps with resolutions of 250 m of SoilGrid and 50 km of FAO for recalculation. Thirdly, to further estimate the uncertainty, we used the Monte Carlo simulation to quantitatively test the data source uncertainty and variability (soil types) into the final results. The results showed that the relative uncertainty of soil P stock was estimated to be 1.22% and 2.60%, respectively, between the 1980s and 2020s (Fig. R5). The relative uncertainty of the change in soil P stock from the 1980s to the 2020s was calculated to be 1.89% (Fig. R5). The data indicated a relatively low level of variability in soil P stock between the three soil type selection criteria. An ensemble calculation of the three-resolution soil maps (30 m, 250 m and 50 km) was adopted to reduce the uncertainty of the soil types. Please see Revised Manuscript Page 15, line 380-383 and Fig. R5 below.\nDynamics of vegetation types. There may have been some subtle changes in ecosystem types over the past three decades at extremely high elevations that are sensitive to climate change, and in areas where human activity is intense. For example, urban area expansion, afforestation, returning cropland to forests or grassland, and the rise of forest line and grass line could result in interconversion of ecosystems on the Tibetan Plateau. According to your suggestion, we refer to the vegetation maps of the 1980s and 2020s to calculate the soil P stocks in the revised version. Please see Revised Manuscript Page 15, line 383-384.\nDynamics of soil type. We appreciate the insights about the limitations of the dynamics of the soil type. Soil development is closely related to local climate, vegetation, soil parent material and land use change. Under the influence of climate change and human activities, soil types on the Tibetan Plateau are likely to have changed in the past few decades. Specifically, in the lake and reservoir inundation areas, glacier retreat areas where cropland has been returned to forest or grassland, there may be transformations in the ecosystem, leading to changes in spatial distribution for specific soil types. Thus, we have tried to synthesize long-term dataset or observations for soil types from the existing studies but failed to find recent data. This kind of long-term dataset is very scarce. Frankly, it is still a big challenge to obtain high-precision spatial data of soil types for the two periods over a broad geographic scale, probably because monitoring soil type changes requires extensive field investigation. For example, as we described above, the production of soil type maps in the 1980s took 16 years and involved more than 80,000 researchers. Due to this point, together with your comments, we added the potential uncertainty of changes in soil type on soil\nThe Limitation is as follows: “There may have been some subtle changes in soil types over the past three decades. For example, afforestation, conversion of cropland to forest and grassland, and an upward shift in forest lines could result in interconversion of ecosystems and soil types on the Tibetan Plateau. The lack of continuous, high-precision data on spatial soil types may have led to some bias in our estimation of historical soil P stock changes.”\n[Comment 5] The sensitivity analysis of the process-based models should be conducted to quantify the uncertainty contributions of the various parameters to the models.\n[Response 5] Very good comment. Following your insightful comment, we conducted the sensitivity analysis of various parameters in the process-based models. We employed Monte Carlo simulations to assess the sensitivity of MPFE model outputs to activity data (including statistical and remote sensing data) and parameter variations. To achieve this, we provided the ranges for all activity data and parameters instead of using single-point estimates. We generally assumed that input parameters are independent, allowing each parameter to be sampled independently from its distribution. Assuming a normal distribution for probability distributions, the variability range of certain parameters was based on multiple traceable reference sources, while parameters from a single reference were quantitatively set according to differences in their data's temporal scales. The simulation model conducted 10,000 trials by randomly selecting values from the input distributions, generating sensitivity range results that include the mean, median, standard deviation (SD), coefficient of variation (CV), and the 5th, 25th, 75th, and 95th percentiles. The results indicated that the sensitivity affecting net P flow primarily originated from parameters rather than activity data, and that activity data and parameters related to soil erosion had a significant impact on net P flow (Supplementary Table 8). Please see Supplementary Text 2 (Revised Supplementary Information page 21-22, line 561-575) and\n[Comment 6] If soil erosion explains most of the P depletion, it must be pointed out that during the last several decades overall soil erosion become less severe as vegetation coverage as a main control of soil erosion has become significantly better.\n[Response 6] Sorry for the confusion. Water erosion is the main factor leading to the decrease of soil P stocks across the Tibetan Plateau. There are many factors affecting water erosion, including the rainfall erosivity factor, soil erodibility factor, slope length, steepness factor, the cover and management factor and ect. In the past few decades, the Tibetan Plateau has undergone a process of climate warming and wetting (Kuang and Jiao, 2016), and has faced an intensified threat of water erosion due to the increasing rainfall erosivity factor (Tian, et al., 2020). However, the warming and humidification of climate can also lead to an increase in vegetation coverage and reduce water erosion (Ying, et al., 2024). Thus, climate warming and wetting can result in mixed and unclear responses of water erosion.\nThe changes in climate and vegetation conditions are suggested as the major contributors to change in water erosion (Yetemen, et al., 2019). In the present study, we found an insignificant change in inter-annual erosion over years (P=0.30), while significantly declining (P<0.01) and rising (P<0.05) changes were detected in the contemporaneous C (combined vegetation) and R (rainfall erosivity) factors, respectively (Supplementary Fig. 9). The insignificant change in the amount of water erosion across the Tibetan Plateau may come from the counteracting effect of C and R factors. Our results are consistent with previous studies (Sun, et al., 2019; Ying, et al., 2024). Therefore, the outflow of soil P caused by water erosion decreased mainly due to the continuous reduction in soil P concentrations rather than the change in the amount of water erosion. To avoid ambiguity, we have rephrased the related sentence and discussed the climate warming and wetting, and vegetation coverage on soil P loss in the Results and Discussion section of the revised version. Please see Revised Manuscript page 7, line 169-174; Revised Supplementary Information\nSupplementary Fig. 9.\n[Comment 7] The introduction section should summarize the current understanding of soil P changes in the Tibetan Plateau. It will help to identify the key knowledge gaps that are addressed in this study.\n[Response 7] Thanks for the valuable comment. Following your insightful suggestion, we summarized the current understanding of soil P changes on the Tibetan Plateau to accurately identify the key knowledge gaps in the Introduction section. Please see Revised Manuscript page 3, line 66-71.\nThe Introduction is as follows: “A few studies have used observations from non-resampling methods to explore the historical changes in total P and available P (AP) in alpine soils. However, the quantity of soil P stocks that has been lost across Tibetan alpine permafrost regions over recent decades remains unclear, and the P budget of permafrost has not been well assessed. Moreover, substantial uncertainties exist in the projected direction and magnitude of soil P changes, originating from the scarcity of observational broad-scale soil resampling data.”.\nMinor comments:\n[Comment 8] L22: One decimal place is enough.\n[Response 8] Thank you for the careful review. Corrected accordingly across the revised version.\n[Comment 9] L51-53: I don’t think it’s a drawback of soil sampling. We cannot demand other irrelevant information from soil measurements.\n35\n[Response 9] Thank you for this suggestion. We have revised this sentence, please see Revised Manuscript pages 3, lines 54-55.\n[Comment 10] L78: Four decades?\n[Response 10] Historical soil data were obtained from the Second National Soil Inventory of China and the First Tibetan Plateau Scientific Expedition Program during the 1980s, with an average year of 1985 (Supplementary Table 1). The repeat sampling campaign was conducted in the 2020s based on historical soil sampling plots from the 1980s. The resampled soil P data were less than 40 years, and thus we used three decades in the present study. Thank you for your understanding!\n[Comment 11] L88-90: It is better to briefly explain the reason for this spatial pattern.\n[Response 11] Thanks for the suggestion. We updated the reason here, please see Revised Manuscript page 4, line 99-100.\n[Comment 12] L118: “The proven MPFE model suggested that the…” - replaced by “The…”\n[Response 12] Thanks for the suggestion. Revised accordingly, please see Revised Manuscript page 6, line 129.\n[Comment 13] L163-165: Is it possible to model or incorporate the information on eating preference and population in the future? The livestock products act as important P outflow.\n[Response 13] Thanks for the helpful suggestion. From the perspective of food consumption structure, the rural residents on the Tibetan Plateau still maintain the traditional eating habits and self-sufficient food consumption (Wang, et al., 2021). The coarse grain consumption is mainly highland barley, and the fine grain is mainly wheat. The consumption of meat is high, while the consumption of other foods such as fruits and vegetables are insufficient. Over the past several decades, major changes in human\n36\ndiets have evolved into per capita increases in total caloric intake and in meat consumption (Cui, 2019). Due to the large differences stemming from the inefficient conversion from feed to meat (less than 10%, compared to 30% for crop production), an increasingly animal-based diet implies a significant demand for P consumption (Chen and Graedel, 2016). These changes in human diet and population may affect the calculation of the P budget. Although four SSP scenarios (including the SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5 future climate scenarios of the CMIP6 models) include food and meat consumption and demographic data (Table R3), very detailed data on the Tibetan Plateau are still lacking. Therefore, we have added these Future improvements in the revised version. Please see Revised Manuscript pages 11, lines 272-276.\nThe Future improvements is as follows: “Due to the large differences stemming from the inefficient conversion from feed to meat (less than 10%, compared to 30% for the crop production), changes in dietary structure across the Tibetan Plateau may lead to changes in P demand. More accurate data on eating preferences and populations need to be incorporated into the MPFE model in the future”.\nTable R3 Assumption on key influencing factors of diet and population projection for SSP scenarios\nClimate\nDiet (Popp, et al., 2017) Population (Huang, et al., 2019b)\nscenarios\nLow growth in food consumption, Fertility:low; Mortality:low; Migration:\nSSP1\nlow meat diets medium; Education: high\nMaterial-intensive consumption, Fertility: medium; Mortality: medium;\nSSP2\nmedium meat consumption Migration: medium; Education: medium\nElites:high consumption lifestyles; Fertility: low; Mortality: medium;\nSSP4\nrest; low consumption Migration: medium; Education: low\nMaterial-intensive consumption; Fertility:low; Mortality:low; Migration:\nSSP5\nmeat-rich diets high; Education: high\nOverall, we appreciate the reviewers’ insightful comments. These comments enabled us to have a comprehensive overview of published literature in P dynamic in permafrost regions, and a deeper consideration of the P dynamics in different ecosystems, field sample methods, results interpretation, sensitivity analysis, uncertainty estimation and future directions, and guided us to thoroughly revise the manuscript. By addressing these comments, we feel that the revised manuscript has been greatly improved. Thank you!\n(https://www.globaltimes.cn/page/202408/1318238.shtml, 2024).\n40\nand demand balance in the Qinghai-Tibet Plateau in recent 30 years. Acta Ecol. Sin. 44, 1-16 (2024).\n42"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I would like to express my appreciation for the outstanding revisions the authors have made to their manuscript. The revisions are truly remarkable and have significantly enhanced the quality of the manuscript. The revision notes are exceptionally clear and well-organized, which makes it easy to understand the improvements and the rationale behind them. The modifications have addressed the key issues raised in the previous review process, and the clarity, coherence, and overall presentation of their manuscript have been greatly improved. The arguments are now more compelling, and the data are presented in a more convincing manner. Now it has resulted in a much stronger and more impactful piece of work.\nI am very pleased with the outcome of the revisions and believe that the manuscript is now ready for publication.", "reviewer_3": "General comments\nThis manuscript presents a significant contribution to the understanding of soil phosphorus dynamics in the Tibetan alpine permafrost regions over recent decades. By utilizing a combination of soil resampling methods and a modified process-balanced model, the authors quantify historical changes in soil P stocks and flows from the 1980s to the 2020s and provide projections for the year 2100. The findings highlight a dramatic decline in soil P stock and raise concerns regarding future phosphorus scarcity in these ecosystems. The manuscript is quite well written and well presented.\nHowever, more detail is needed concerning the selection process for the resampling plots. Addressing potential biases introduced by environmental changes or human activities during the three-decade gap is needed to enhance the robustness of the conclusions. It would be helpful to discuss potential uncertainties associated with these scenarios or conditions that may affect the projections, including climate feedback mechanisms.\nThe authors should also consider that as more water flows through the system, there is increased phosphorus loss, which correlates with greater soil erosion. If vegetation increases in presence but soil remains, phosphorus can become increasingly limiting. However, if soils are completely lost and new soil formation begins from bedrock, phosphorus availability may shift as it comes from rocks, while nitrogen and water may become more limiting factors. Therefore, the authors should clearly state the hypothesis regarding the relative importance of phosphorus, nitrogen, and water across different soil degradation scenarios.\nAdditionally, it would be beneficial to provide a more detailed description of the study sites, including specific geographical and environmental characteristics, as well as a thorough account of the sampling methods and data analysis techniques used. Addressing concerns about site comparability and the limitations imposed by the number of data points would strengthen the manuscript.\nMoreover, the implications of phosphorus scarcity focus primarily on terrestrial and aquatic ecosystems but it would be good to further explore ramifications on regional agricultural practices. Including a discussion on potential adaptive strategies for local communities may enhance applicability. The authors should expand on the discussion surrounding agricultural practices and their potential impact on soil recovery, integrating how these practices could influence the dynamics of phosphorus, water and nitrogen in the context of soil degradation and recovery.\nDetailed comments\nRegarding your statements on impacts on water bodies and algal blooms: to make your statements robust and convincing , you must know and describe what level of P concentrations reach the body waters, if not is very speculative.\nClarify which tipping point when you state: “The tipping point for the P budget is around 2050. Soil….”\nWhen you state limitations, please take into account that the limiting role of phosphorus in the future will depend on the remnants of soil that remain after permafrost completely disappears, as these will serve as the initial substrate for plant establishment and the onset of the edaphogenic process. If future plant colonization occurs on this residual soil that lacks frost but is depleted in phosphorus, phosphorus may become a significant limiting factor. However, if this legacy soil is lost to erosion and the edaphogenic process begins with rock weathering and leaching, other resources, such as nitrogen or water, are likely to become more limiting during the initial thousands of years of soil development as commented in previous General comments section of my review above. As comented there, all your statements depend on the capacity of retaining soil.\nIn your statement on “a higher N deposition rate can stimulatephosphatase activities and thus increase the inorganic P supply, which could indirectly promote the bioavailability of soil P….” take into account that this also increases the potential to be leached.\nFig R4: It should be noted that in forested areas, the available phosphorus seems to be very high, while in pasture areas it is very low. This should be stated and discussed, likely due to the more favorable climate and soil protection, as well as having a larger reservoir in the vegetation, which has resulted in less loss of phosphorus.\nRegarding the response letter itself:\nThe response letter demonstrates a good effort to address the reviewers' comments. The authors have made significant revisions to the manuscript, which is commendable. However, some responses could be more detailed and specific. For example:\nReviewer 1:\n[Response 2]: This response is a bit vague. What specific sentence did you rephrase, and what additional information did you add?\n[Response 7]: What specific studies did you compare your results to? What specific details did you compare (e.g., magnitude of soil P loss, spatial patterns)?\nRecommendation:\nAccept after considering the suggested revisions.\nProf. Josep Penuelas"}, "author_response": "Reviewer #3 (Remarks to the Author):\nGeneral comments\n[Comment 1] This manuscript presents a significant contribution to the understanding of soil phosphorus dynamics in the Tibetan alpine permafrost regions over recent decades. By utilizing a combination of soil resampling methods and a modified process-balanced model, the authors quantify historical changes in soil P stocks and flows from the 1980s to the 2020s and provide projections for the year 2100. The findings highlight a dramatic decline in soil P stock and raise concerns regarding future phosphorus scarcity in these ecosystems. The manuscript is quite well written and well presented.\n[Response 1] Thank you for your positive comment.\n[Comment 2] However, more detail is needed concerning the selection process for the resampling plots.\n[Response 2] Very good comment. To address your comment about the sampling strategy, we provided more details about the selection process for the resampling plots.\nFirstly, we conducted statistics on the counties where the historical plots were located and the counties where the resampling plots were selected. In addition, we also analyzed the distribution of the climate conditions of the counties where the resampling plots were selected. Please see Revised Manuscript page 14-15, line 377-379; Supplementary Fig. 1 (b).\nThe added sentence in the revised MS is as follows:\n“…Historical data in the 1980s were recorded from 131 counties across the Tibetan Plateau, and resampling plots were obtained from 93 counties in the 2020s, accounting for 71% of the total 131 counties (Supplementary Fig. 1, b)…”\nThe added figure in the Revised Supplementary Information (Supplementary Fig. 1 b) is as follows:\nSecondly, we added more detailed information on selection process for the resampling plots, such as comparison with other independent literature, consultation with several staff who participated in the past campaign, and avoiding areas disturbed by intense human activities in the Method section. Please see Revised Manuscript page 15, line 379-386.\nThe added sentences in the revised MS are as follows:\n“…To further evaluate the accuracy of our plots, we compared all the plots with other independent references. We made secondary use of different maps, including vegetation-type maps, soil-type maps, topographic maps (elevation, slope, and terrain), land-use maps (residential areas, rivers, and lakes), road maps, and administrative-boundary maps that were developed before 1990s (more details see Supplementary text 1). In addition, we also consulted several staff who participated in the 1980s campaign to determine the exact locations of the sampling plots. Plots in close proximity to electric power lines, roads, railways, or drainage ditches were avoided…”\nThirdly, we added representativeness of selected soil samples according to the soil P density. A histogram of soil P density of the 255 resampling plots (red bars) in comparison to the historical dataset (605 plots, black bars) was provided. Please see Revised Manuscript page 14, line 374-377; Supplementary Fig. 2.\nThe added sentences in the revised MS are as follows:\n“…Frequency distribution analysis of SPD revealed similar means between the original and resampled plots. A mixed-effects model further confirmed no significant difference between the two datasets (P > 0.05). The distribution and representativeness of the plots is depicted in Supplementary Fig. 2.”\nThe added figure in the Revised Supplementary Information (Supplementary Fig. 2) is as follows:\nFig. R2 Representativeness of selected soil plots. (a) Histogram of soil P density of the 255 sub-plots (blue bars) in comparison to the whole 605 plots (red bars). (b) The comparisons of soil P density between recorded plots and resampled plots. ns, non-significant differences.\n[Comment 3] Addressing potential biases introduced by environmental changes or human activities during the three-decade gap is needed to enhance the robustness of the conclusions. It would be helpful to discuss potential uncertainties associated with these scenarios or conditions that may affect the projections, including climate feedback mechanisms.\n3\n[Response 3] Good suggestions. Following your insightful suggestions, we analyzed the influences of climate changes and human activity on soil P changes using the linear mixed-effects model. At least over the past three decades, changes in MAT ( MAT) and human activity ( Human activity) had marginal effects on soil P changes, w△hile MAP, changes in MAP△ ( MAP) and MAT had significant effects on soil P changes. High precipitation and i△ncreased precipitation were expected to have negative impact on soil P density by accelerating water erosion (Fig. R3). However, high temperature could enhance mineral weathering and may potentially promote the soil P availability (Fig. R3). Although the population has increased significantly during the last three decades, the density of the population of farmers and herders is still relatively low on the Tibetan Plateau, resulting in a low consumption of P (Fig. R3).\nBased on the results of linear mixed-effects modeling, We have added the information concerning the effects of climate changes and human activities on soil P changes in the revised MS. Please see Revised Manuscript Pages 11-12, lines 286-300; Supplementary Fig. 15.\nThe supplementary discussion in the revised MS are as follows:\n“Our analysis presents a potential soil P crisis in the Tibetan alpine permafrost region, but the P cycle process is extremely complex, with both climatic factors (i.e., warming and humidification) and human activities possibly causing potential prediction uncertainties. For instance, the warming climate of the Tibetan Plateau is expected to enhance mineral weathering and may potentially promote the availability of soil P, while increasing precipitation (in both amount and magnitude) is expected to have negative impacts on SPD by accelerating water erosion (Supplementary Fig. 15). Some studies have found that plants can alleviate the limitation of P by promoting mycorrhizal symbiosis (Franke, et al., 2018), activating mineral-bound P and uplifting P (biological pumping) from the deep soil layer (Lucas, 2001). It is very difficult for our model to capture these complex climate feedback mechanisms or ecological processes. Although the overall lower intensity of human activities has a relatively small impact on soil P changes (Supplementary Fig. 15), the local construction of highways/railways and mining may increase surface disturbances (Chen, et al., 2022),\n4\nthereby accelerating erosion and the migration rate of P. For those variables, regional-scale data are not yet available for the Tibetan Plateau, which may also cause potential prediction uncertainties.”\nThe added figure in the Revised Supplementary Information (Supplementary Fig. 15) is as follows:\nFig. R3 Dots and lines represent coefficients and 95% confidence intervals in the linear mixed effects model. MAT, mean annual temperature; MAP, mean annual precipitation; ΔMAT, changes in MAT in the past three decades; ΔMAP, changes in MAP over the last three decades. Human activity refers to product and biomass energy consumption. ΔHuman activity, changes in product and biomass energy consumption over the last three decades. ns, non-significant; *, P < 0.05; **, P <.0.01; ***, P < 0.001.\n[Comment 4] The authors should also consider that as more water flows through the system, there is increased phosphorus loss, which correlates with greater soil erosion. If vegetation increases in presence but soil remains, phosphorus can become increasingly limiting. However, if soils are completely lost and new soil formation begins from bedrock, phosphorus availability may shift as it comes from rocks, while nitrogen and water may become more limiting factors. Therefore, the authors should clearly state the hypothesis regarding the relative importance of phosphorus, nitrogen, and water across different soil degradation scenarios.\n5\n[Response 4] Very good comment. As you said, we previously overlooked the relative importance of P, N and water at different stages of soil degradation. Following your insightful comment, we have added supplementary discussion in the revised version. Please see Revised Manuscript page 12, line 301-308.\nThe added discussion in the revised MS is as follows:\n“It should be pointed out that the relative importance of P and other resource limitations (e.g., N and water) varies under different soil degradation scenarios. For example, substantial soil P was translocated to plant biomass (grassland root and forest biomass), with the soil remaining well-preserved, and P limitation may become intense in the presence of a continuous decline in soil P. If soils are completely lost and new soil formation begins from bedrock, however, P availability may increase due to the weakened soil shielding effect and increased rock weathering caused by the frequent freeze–thaw effect (Hartmann, et al., 2014), thereby leading to N and water possibly becoming more limiting factors.”\n[Comment 5] Additionally, it would be beneficial to provide a more detailed description of the study sites, including specific geographical and environmental characteristics, as well as a thorough account of the sampling methods and data analysis techniques used. Addressing concerns about site comparability and the limitations imposed by the number of data points would strengthen the manuscript.\n[Response 5] Thank you for your valuable suggestions. According to your valuable comments, we have added more detailed description of the study plots, sampling methods, data analysis techniques, and site comparability and the limitations imposed by the number of data points.\n1. Description of the study plots.\nFollowing your suggestions, we provided a detailed description of the study plots, including soil P data, specific geographical and environmental characteristics in Supplementary Dataset 3. Please see Revised Manuscript page 22, line 576-578; Supplementary Dataset 3.\nThe summary of the environmental information of all the plots is as follows:\n“These plots were set throughout the geographical extent of alpine forests, grasslands and cropland, covering a wide range of climates (mean annual temperature: −12.6 to 14.3 °C; mean annual precipitation: 22.4–1299.5 mm), vegetation characteristic (net primary productivity: 6.1–1173.7 g C m−2 yr−1; normalized difference vegetation index (NDVI): 0–0.8), soil properties (soil P concentration: 0.2–1.6 g kg−1; soil AP concentration: 0.4–97.6 mg kg−1) across the plateau. The topography, parent material, soil type, vegetation type, dominant species, and management strategies of the 255 resampling plots did not change during the intra-sampling period.”\n2. Description of sampling methods.\nFirstly, we added a schematic diagram of sampling design in the Supplementary Figure 1 (c).\nThe added Supplementary Fig. 1 (c) is as follows:\nFig. R4 Schematic diagram of sampling design. Three square pits were dug in each plot.\nSecondly, we supplemented more information on sampling process in the Method section. Please see Revised Manuscript page 15, line 397-401.\nThe supplementary information in the revised MS is as follows:\n7\n“We took three soil cores in each plot to generate descriptions of their vertical properties (e.g., layering, color, and texture for each layer) for comparisons with those from the 1980s. The consistency between each layer of the profile in the 2020s and the historical record was the key piece of evidence used to determine whether the soil sample was collected correctly.”\nThirdly, we provided more additional instructions for field resampling processes and plots selection processes in Supplementary information. Please see Supplementary Text 1.\n3. Description of data analysis techniques.\nFollowing your suggestions, we provided more description of data analysis techniques in the revised version, including frequency distribution analysis (Revised Manuscript page 14, line 374-375; Supplementary Fig. 2 a), mixed-effects model (Revised Manuscript page 14, line 375-376; Supplementary Fig. 2 b), bootstrap resampling approach (Supplementary Fig. 4), and linear mixed-effects modeling (Supplementary Fig. 15).\nIn addition, we have added code associated with the data analysis techniques in the revised MS. The code of data analysis have been deposited in the figshare database under accession code https://figshare.com/articles/code/29134574. Please see Revised Manuscript page 22, lines 581-582.\n4. Site comparability and the limitations imposed by the number of data points\nAccording to your comments, we added representativeness of selected soil samples according to the soil P density. A histogram of soil P density of the 255 resampling plots (red bars) in comparison to the historical dataset (605 plots, black bars) was provided. Please see Revised Manuscript page 14, line 371, 374-377; Supplementary Fig. 2.\nThe added sentences in the revised MS are as follows:\n8\n“The original soil P dataset comprised 605 plots in the 1980s.”\n“Frequency distribution analysis of SPD revealed similar means between the original and resampled plots. A mixed-effects model further confirmed no significant difference between the two datasets (P > 0.05). The distribution and representativeness of the plots is depicted in Supplementary Fig. 2.”\nTo determine the robustness of observed changes in soil P density, we applied a bootstrap resampling approach following Crowther, et al. (2016). This involved iteratively (n = 10,000 repetitions) removing randomly selected plots from the dataset and recalculating the slope of the linear mixed-effects model (LMM) for soil P density trends using the remaining plots. For each iteration, we recorded the slope estimates and derived their interquartile ranges (IQRs) across different levels of plot removal. Statistical significance was assessed by determining whether the IQR of the bootstrapped slopes overlapped with zero (i.e., the null line), with non-significance indicating that the observed trend could not be reliably distinguished from no change (Crowther, et al., 2016). Our results demonstrated that the statistical significance of the average slope persisted until more than 240 plots were removed from the initial dataset (n=255). This robustness decreasing pattern held true in alpine grassland, forest and cropland ecosystems, respectively. These additional analyses demonstrated the robustness of our conclusion about the decreased soil P on the Tibetan Plateau over the detection period. Please see Revised Manuscript page 4, line 93-99; Supplementary Fig. 4.\nThe added result in the revised MS is as follows:\n“Further analyses revealed that the SPD tended to decrease among 190 of the total 255 resampling plots and the decline in SPD remained significant until >240 resampling plots were removed from the initial analysis with 255 resampling plots (Supplementary Figs. 3 and 4), demonstrating the robustness of SPD dynamics observed in this study. This decreasing pattern held true in alpine grassland, forest and cropland ecosystems (P < 0.01), with the net loss of SPD estimated to be 533.5, 414.8 and 715.0 kg P ha−1, respectively (Fig. 2c, Supplementary Figs. 3 and 4).”\nFig. R5. Bootstrapped estimates of model slope based on different number of resampling plots, of which the model refers to the linear mixed-effects model used to statistically analyze changes in soil P density on the Tibetan Plateau over the period from the 1980s to the 2020s. (a) Overall, (b) Grassland, (c) Forest, (d) Cropland. This analysis estimated the threshold at which soil P density changes lost statistical significance when plots were randomly removed plots from the original dataset. In the analyses, the bootstrap was iterated for 10,000 times, in each of which the resampling\n10\nplots were removed randomly and the slope of the linear mixed-effects model used to test changes in soil P density was estimated based on the remaining plots (Crowther, et al., 2016). The estimated slopes were then used to calculate interquartile range for each number of resampling plots that were removed from the dataset (Crowther, et al., 2016). The average slope was considered non-significant if the interquartile range of the slope overlapped with the horizontal zero line (Crowther, et al., 2016). The dashed line represents the first interquartile range which overlapped with the horizontal zero line. Each point in the diagram represents a bootstrapped estimate of model slope, and each box-plot presents the interquartile range and median slope for each number of resampling plots removed from the analysis.\n[Comment 6] Moreover, the implications of phosphorus scarcity focus primarily on terrestrial and aquatic ecosystems but it would be good to further explore ramifications on regional agricultural practices. Including a discussion on potential adaptive strategies for local communities may enhance applicability. The authors should expand on the discussion surrounding agricultural practices and their potential impact on soil recovery, integrating how these practices could influence the dynamics of phosphorus, water and nitrogen in the context of soil degradation and recovery.\n[Response 6] Very good comment. Regarding the agricultural practices in cropland ecosystems, we are sorry that our original manuscript did not consider it comprehensively. Adding a discussion on potential adaptive strategies for local communities could enhance applicability in cropland ecosystems. Following your insightful suggestions, we have supplemented discussion on potential adaptive strategies in the revised version. Please see Revised Manuscript page 10-11, line 258-270.\nThe supplementary discussion in the revised MS is as follows:\n“Agricultural practices may be able to reduce P loss by adopting conservation agriculture, sustainable farming practices, increasing the utilization of animal manure, or using management systems with increased use of residue inputs (MacDonald, et al., 2011). However, in the long-term, a reduction in P losses through mitigating soil erosion (e.g., mulching, increased vegetation cover, intercropping and conservation agriculture), will be the most effective approach, and it will simultaneously enhance soil health, the general nutritional status of the soil (e.g., P and N), and water retention capacity (Alewell, et al., 2020). Although leguminous plants will mainly improve the N rather than P status of the soil, a recent review found that leguminous crops can utilize less accessible forms of P under P-deficient conditions via several mechanisms, including shaping a more refined root structure, releasing root exudate, and strengthening the connection with mycorrhizae (Franke, et al., 2018). Alfalfa (Medicago sativa L. Alfking), a leguminous plant, is a high-quality forage type widely cultivated on the Tibetan Plateau (Cui, et al., 2017), and thus growing/intercropping with alfalfa may be one of the effective ways to alleviate P limitation in cropland. ”\nDetailed comments\n[Comment 7] Regarding your statements on impacts on water bodies and algal blooms: to make your statements robust and convincing , you must know and describe what level of P concentrations reach the body waters, if not is very speculative.\n[Response 7] Very good comment. Following your suggestions, we have added the descriptions of total P concentrations and microorganisms changes in some permafrost lakes. Please see Revised Manuscript page 9, line 233-236.\nThe supplementary discussion in the revised MS is as follows:\n“…For instance, previous studies found that total P concentrations in 7 lakes have exceeded the eutrophication threshold set by the OECD (>35 mg L−1)(Hart, 1984; Wu, et al., 2021), as well as exhibited a significant increasing trend for both prokaryotic and eukaryotic microorganisms over the past several decades (Ouyang, et al., 2023)….”\n[Comment 8] Clarify which tipping point when you state: “The tipping point for the P budget is around 2050. Soil….”\n[Response 8] Sorry for the confusion. The tipping point for the P budget is around 2050. Soil P will continue to decrease rapidly from the 2020s to the 2050s, with an average\n12\nrate of 10.0 kg P ha−1 yr−1. The projected net loss of soil P will then gradually slow down, with the mean rate decreasing to 6.8 kg P ha−1 yr−1 by the end of the 21st century. Please see Revised Manuscript page 8, line 198-202.\nTherefore, we have clarified these sentences in the revised MS:\n“…The results of the change-point tests showed that soil P will continue to decrease rapidly from the 2020s to the 2050s, with an average rate of 10.0 kg P ha−1 yr−1 (Fig. 5d). The projected net loss of soil P will then gradually slow down, with the mean rate decreasing to 6.8 kg P ha−1 yr−1 by the end of the 21st century (Fig. 5d), mainly due to a decrease in the amount of P outflow via erosion…”. Thank you for your understanding!\n[Comment 9] When you state limitations, please take into account that the limiting role of phosphorus in the future will depend on the remnants of soil that remain after permafrost completely disappears, as these will serve as the initial substrate for plant establishment and the onset of the edaphogenic process. If future plant colonization occurs on this residual soil that lacks frost but is depleted in phosphorus, phosphorus may become a significant limiting factor. However, if this legacy soil is lost to erosion and the edaphogenic process begins with rock weathering and leaching, other resources, such as nitrogen or water, are likely to become more limiting during the initial thousands of years of soil development as commented in previous General comments section of my review above. As commented there, all your statements depend on the capacity of retaining soil.\n[Response 9] Thanks for the helpful suggestion. According to your previous General comments 4, we have added this limitation in the revised version. Please see Revised Manuscript page 12, line 301-308.\nThe added limitation in the revised MS is as follows:\n“It should be pointed out that the relative importance of P and other resource limitations (e.g., N and water) varies under different soil degradation scenarios. For example, substantial soil P was translocated to plant biomass (grassland root and forest biomass), with the soil remaining well-preserved, and P limitation may become intense in the presence of a continuous decline in soil P. If soils are completely lost and new\n13\nsoil formation begins from bedrock, however, P availability may increase due to the weakened soil shielding effect and increased rock weathering caused by the frequent freeze–thaw effect (Hartmann, et al., 2014), thereby leading to N and water possibly becoming more limiting factors.”\n[Comment 10] In your statement on“a higher N deposition rate can stimulate phosphatase activities and thus increase the inorganic P supply, which could indirectly promote the bioavailability of soil P….” take into account that this also increases the potential to be leached.\n[Response 10] Thank you for this suggestion. According to your comments, we\nhave added the leaching process in the revised version. Please see Supplementary Information page 27, line 685.\nThe added sentence in the revised MS is as follows:\n“However, the possibility of AP being leached has also increased.”\n[Comment 11] Fig R4: It should be noted that in forested areas, the available phosphorus seems to be very high, while in pasture areas it is very low. This should be stated and discussed, likely due to the more favorable climate and soil protection, as well as having a larger reservoir in the vegetation, which has resulted in less loss of phosphorus.\n[Response 11] Thanks for the valuable comment. Inspired by your comment, we used one-way analysis of variance to compare the AP densities among different ecosystems. As you guessed, the soil AP of the forest ecosystem is significantly higher than that of the grassland ecosystem in the 2020s. We have supplemented the discussion on this issue. Please see Supplementary Information page 27, line 674-677; Supplementary Fig. 25.\nThe added sentence in the revised MS is as follows:\n“The soil AP in forest was significantly higher than that in grassland in the 2020s, which may be related to the more favorable climate, higher vegetation coverage and litter yield, as well as the better soil protection effect in forest.”\n14\n[Comment 12] Regarding the response letter itself:\nThe response letter demonstrates a good effort to address the reviewers' comments.\nThe authors have made significant revisions to the manuscript, which is\ncommendable. However, some responses could be more detailed and specific. For\nexample:\n① Reviewer 1, Response 2: This response is a bit vague. What specific sentence did you rephrase, and what additional information did you add?\n[Response 12] Thank you very much for the careful review. Sorry for the confusion.\nAccording to your comments, here we have clarified the added and rephrased\nsentences in the Introduction. Please see Revised Manuscript page 2-3, line 40-42,\nline 44-46, line 48-52.\nThe added sentence in the revised MS is as follows:\n“The unidirectional flow of P from terrestrial ecosystems to freshwater bodies and,\neventually, the oceans is considered to be the greatest sustainability challenge in the P\ncycle (Alewell, et al., 2020).”\nThe rephrased sentences in the revised MS are as follows:\n“…In addition, the availability of P alters the trajectory of the permafrost carbon (C)\ncycle by influencing gross primary productivity and net ecosystem productivity (Zimov,\net al., 2006; Yang, et al., 2021)…”\n“Therefore, quantifying the temporal changes of the P stored in permafrost is crucial\nand urgently needed in guiding strategies to tackle the grand challenge of future P\nlimitation (Alewell, et al., 2020; Barbieri, et al., 2022), as well as the environmental\ncrisis caused by huge losses of P into water bodies that pushes the Earth's environment\nto its boundary, beyond which come abrupt environmental changes with catastrophic\nimpacts (Chen and Graedel, 2016; Brownlie, et al., 2022).”\n② Reviewer 1, Response 7: What specific studies did you compare your results to?\nWhat specific details did you compare (e.g., magnitude of soil P loss, spatial patterns)?\n[Response 12] Thank you very much for the careful review. We are very sorry for the\n15\nabsence of some details in the reply letter, but we have made modifications and\nsupplements in the revised version. We reorganized the statements and responded\nto this issue. In order to verify the reliability of the results, we compared the direction\nand magnitude of soil P changes with previous studies. Please see Revised\nManuscript page 4-5, line 99-101.\nThe added comparative analysis in the revised MS is as follows:\n“The rapid decrease in soil P (14.9 kg ha−1 yr−1) across the Tibetan Plateau was within\nthe range of the global soil P depletion rate (4–19 kg ha−1 yr−1)(Alewell, et al., 2020).”\nSince the largest part of soil P loss is caused by water erosion, we also compared the\ndata of P loss caused by water erosion on the Tibetan Plateau with other studies. Please\nsee Revised Manuscript page 7, line 162-165 and Revised Supplementary\nInformation Supplementary Table 6.\nThe outflow of P caused by water erosion between the results of this study and\nthose of previous studies in the revised MS is as follows:\n“Our area-related calculations resulted in an average outflow of P from the alpine\npermafrost region caused by water erosion of about 16.1 kg P ha−1 yr−1 (Supplementary\nFig. 10a), which is slightly higher than that of previous calculations worldwide\n(Supplementary Table 6).”\nThe table in the Revised Supplementary Information (Supplementary Table 6)\nis as follows:\nTable R1 P loss from soil systems to waters as discussed in recent literature. For\ncomparability, all values have been normalized to kg P ha−1 yr−1.\nSoil P outflow (kg P\nRegion Description Reference\nha−1 yr−1)\n(Quinton, et al.,\nGlobal 10.2–18.5 Agricultural land\n2010)\n(Chen and\nGlobal 8.2–12 Agricultural land\nGraedel, 2016)\nGlobal 10 Cropland (Smil, 2000)\n(Alewell, et al.,\nAsia 7.8 Cropland soils\n2020)\n(Alewell, et al.,\nChina 12.3 Cropland soils\n2020)\nGrassland, forest and\nChina 5.2–14.02 (Song, et al., 2024)\ncropland soil\nTibetan Grassland, forest and\n16.14 Present study\nPlateau cropland soil\n[Comment 13] Recommendation: Accept after considering the suggested revisions.\n[Response 13] Thanks for the positive comment.\nOverall, we appreciate the editor and reviewers’ insightful comments. These\ncomments enabled us to have a comprehensive overview of published literature in P\ndynamic in permafrost regions, and a deeper consideration of the P dynamics in\ndifferent ecosystems, field sample methods, driving factors and future directions, and\nguided us to thoroughly revise the manuscript. By addressing these comments, we feel\nthat the revised manuscript has been greatly improved. Thank you!\nintensification through rotations with grain legumes in Sub-Saharan Africa: A\nreview. Agr. Ecosyst. Environ. 261, 172-185 (2018).\nImportance of biological pumping. Annu. Rev. Earth Pl. Sc. 29, 135-163 (2001).\nRev. Earth Env. 3, 701-716 (2022).\nchemical weathering and associated P-release - The role of lithology,\ntemperature and soil properties. Chem. Geol. 363, 145-163 (2014).\nwarming. Nature. 540, 104-108 (2016).\nphosphorus imbalances across the world's croplands. P. Natl. Acad. Sci. USA.\n108, 3086-3091 (2011).\nNat. Commun. 11, 4546 (2020).\n4 3\nan alpine agroecosystem. Field Crop. Res. 201, 192-199 (2017).\nOECD. Search. 15, 50 (1984).\nQinghai-Tibetan Plateau. Environ. Pollut. 272, 116032 (2021).\nmicrobial communities in a glacial lake on the Tibetan Plateau. Water Res. 242,\n120213 (2023).\nbudget. Science. 312, 1612-1613 (2006).\ndynamics after permafrost thaw. Global Change Biol. 27, 5818-5830 (2021).\nsystem resilience to phosphorus shortages on a telecoupled planet. Nat. Sustain.\n5, 114-122 (2022).\nproduction, consumption, recycling, and environmental impacts. Global\nEnviron. Chang. 36, 139-152 (2016).\nPhosphorus Future: Towards global phosphorus sustainability. 371 (Centre for\nEcology & Hydrology, Edinburgh, UK, 2022).\nagricultural soil erosion on biogeochemical cycling. Nat. Geosci. 3, 311-314\n(2010).\nenvironmental impacts from a consumption perspective. J. Ind. Ecol. 12, 229-\n247 (2008).\nChina. Global Change Biol. 30, e17108 (2024).\n19"}]

Decision letter