## 1 Introduction  
Phosphorus (P) availability in soils is known to be controlled by its parent material, the stage of soil development, weathering, and erosion [1], as well as soil microbial activity [2]. Many forms of P exist in soils, but plants can acquire only simple inorganic orthophosphate (Pi) [3]. Since Pi is easily fixed by soil particles in soils with a range of chemical properties [4], its levels are usually limiting, which in turn limits the productivity of a variety of ecosystems [5]. It is widely accepted that nutrients limit terrestrial plant production and carbon (C) [6, 7, 8]. Nitrogen (N) has traditionally been considered the most important nutrient that limits plant production in terrestrial ecosystems [8, 9], whereas P has been considered less important, perhaps with exception of lowland tropical forest regions [6, 7, 10]. Thistraditional view, however, has been challenged by mounting evidence of significant P limitation in areas other than the lowland tropical regions, including tundra regions [11; 12] and temperate region areas with strongly weathered soils [13; 14]. Despite the important role of P in terrestrial plant production, we still lack a quantitative understanding the mechanisms underlying terrestrial P limitation on a global scale, its spatial heterogeneity, or interactions of P with C and N [15; 16]. For example, none of the models in the fifth phase of the Coupled Model Intercomparison Project consider P as a limiting nutrient, resulting in substantial uncertainty in estimates of the strength of the terrestrial carbon (C) sink through the 21 st century [15; 16; 17; 18; 19]. Therefore, identifying the P cycle specific processes and quantifying the distribution and spatial heterogeneity of P limitation, as well as its potential effects on ecosystem responses to future climate change, remain a high priority. Integrating the P cycle into terrestrial ecosystem models should be essential, given that laboratory and field experiments, especially empirical stoichiometry observations [20; 21; 22], have already provided theoretical and data foundations [2] that allowed scientists to develop terrestrial ecosystem models that take into account C, N and P processes, interactions and feedback in entire ecosystems. In fact, the early process-based models developed, focused on C cycles of terrestrial ecosystems under climate change and elevated levels of atmospheric CO\({}_{2}\)[15]. A growing number of models that couple C-N interactions have been applied to understanding the effects of global change on multiple scales. The models that most widely applied are the biogeochemical model CENTURY [23; 24; 25; 26], ecosystem process-based model, BGC (BIOMASS-BGC, FOREST-BGC) [27; 28; 29; 30; 31], DNDC [32; 33; 34; 35; 36; 37], and land biogeochemistry model, CLM [38; 39; 40]. These models have been applied to different ecosystems but have been limited primarily to forests, grasslands, and croplands. At global scale, the land carbon models, including LPJ [41], LPJ_GUESS [42], ORCHIDEE [43], CLM4 CN [44], TRIFFID [45], VEGAS [46], SDGVM [47], and Hyland [48], encapsulate natural and human induced processes through which land ecosystems absorb and emit CO\({}_{2}\), and carbon is stored by plants and soils, and have been widely used by the TRENDY modeling project for simulating the energy, water, nitrogen and carbon balances at daily, monthly and yearly time step [49]. Unfortunately, no model has included P cycle and limitation in their model simulations. Efforts have been made to integrate co-limiting N and P on simulations of the productivity of terrestrial ecosystems [6; 50; 51], but the most of second-generation models do not take P limitation into account [15]. In part since this omission creates uncertainty in estimates of vegetation productivity and changes in terrestrial carbon sinks, interest is growing in understanding the effect of P on terrestrial ecosystems, how P limitation affects C and N, and how P and N interact globally [52].