The distribution of soil phosphorus for global biogeochemical modeling

Phosphorus (P) is a major element required for biological activity in terrestrial ecosystems. Although the total P content in most soils can be large, only a small fraction is available or in an organic form for biological utilization because it is bound either in incompletely weathered mineral particles, adsorbed on mineral surfaces, or, over the time of soil formation, made unavailable by secondary mineral formation (occluded). In order to adequately represent phosphorus availability in global biogeochemistry–climate models, a representation of the amount and form of P in soils globally is required. We develop an approach that builds on existing knowledge of soil P processes and databases of parent material and soil P measurements to provide spatially explicit estimates of different forms of naturally occurring soil P on the global scale. We assembled data on the various forms of phosphorus in soils globally, chronosequence information, and several global spatial databases to develop a map of total soil P and the distribution among mineral bound, labile, organic, occluded, and secondary P forms in soils globally. The amount of P, to 50cm soil depth, in soil labile, organic, occluded, and secondary pools is 3.6 ± 3, 8.6± 6, 12.2± 8, and 3.2± 2 Pg P (Petagrams of P, 1 Pg = 1× 1015g) respectively. The amount in soil mineral particles to the same depth is estimated at 13.0 ± 8 Pg P for a global soil total of 40.6± 18 Pg P. The large uncertainty in our estimates reflects our limited understanding of the processes controlling soil P transformations during pedogenesis and a deficiency in the number of soil P measurements. In spite of the large uncertainty, the estimated global spatial variation and distribution of different soil P forms presented in this study will be useful for global biogeochemistry models that include P as a limiting element in biological production by providing initial estimates of the available soil P for plant uptake and microbial utilization.


Introduction
Phosphorus (P) is a primary macronutrient that often limits the productivity and growth of terrestrial ecosystems. Capturing the dynamics of P cycling in ecosystem models is important for understanding and predicting terrestrial carbon dynamics. Although the total P content of soils can be large, in most soils, only a small fraction is available for biological utilization because it is either bound in incompletely weathered mineral particles, adsorbed on mineral surfaces, or, over the time of soil formation, made unavailable by secondary mineral formation (occluded). The lack of information on the spatial distribution of different forms of P has hampered our efforts to incorporate P cycling into coupled biogeochemistry-climate models.
Since soil P transformations occur on geological timescales, an approach that appropriately estimates P status for model initialization is more efficient than modeling P processes at these timescales to arrive at present day conditions. Building on existing knowledge of soil P processes and soil P measurements, we develop a spatially explicit estimate of different forms of P on the global scale and provide a data informed method of initializing P pools for global models. Recently available measurements of different soil P forms and other data such as global lithology maps and rock P concentration have provided us the opportunity for an initial attempt to construct a global map of different P forms. Specifically, 2 maps of surface lithology or rock types are combined with a database of P measurements from surface rock materials to estimate P concentration of the parent material. USDA (United States Department of Agriculture) soil orders are grouped into three categories indicating increasing degrees of weathering intensity. We use 29 soil chronosequences and vertical pedon measurements to Published by Copernicus Publications on behalf of the European Geosciences Union. 2526 X. Yang et al.: The distribution of soil phosphorus estimate a phosphorus pedogenic depletion index (PPDI), an estimate of the accumulated loss of parent material P during pedogenesis. Finally, a recent synthesis of Hedly P fractions from over 170 soil profiles (Yang and Post, 2011) grouped by the USDA soil order allows us to estimate the magnitude of the 5 different forms of P for each combination of parent material and soil order globally. This approach allows a quantitative utilization of currently available data to produce an estimate of initial available P conditions for modeling purposes in a fashion that can readily be improved and uncertainty reduced as more data becomes available.

Background
Our method relies on combining several global datasets and our understanding of P transformations during pedogenesis.
Here we provide some background information about the datasets and the conceptual model we used.
In unmanaged ecosystems, soil P is initially supplied by the weathering of parent material. Parent material of soils includes not only primary bedrock, but also secondary material such as alluvial, aeolian, and colluvial deposits. P concentration of parent material varies considerably, from 140 ppm in carbonate rocks to more than 1000 ppm in volcanic materials (Gray and Murphy, 2002). It follows that the distribution of different parent material exerts a strong control over the soil P status of terrestrial ecosystems (Buol and Eswaran, 2000). Walker and Syers (1976) developed a conceptual model of soil P transformations during pedogenesis based on four choronosequence studies in New Zealand. The model postulates that, at the beginning of soil development, all soil P is in the primary mineral form, mainly as calcium phosphate minerals. With time, primary P minerals weather rapidly, giving rise to phosphorus contained in various other forms (organic P, secondary mineral P, occluded P) with an overall decline of total P due to leaching. The dissolved P from primary minerals can be taken up by biota, entering the organic P reservoir, or sorbed onto the surface of secondary minerals in soils to become "non-occluded P". The non-occluded P is slowly but continually converted to occluded P with time, as P is physically encapsulated or surrounded by secondary minerals (such as Fe and Al oxides). Therefore at the late stages of soil development, soil P is dominated by organic P and occluded P. The Walker and Syers conceptual model suggests that a relationship exists between P content and forms in soils and stage of soil development.
Although Walker and Syers' conceptual model has been well accepted, few studies have focused on the quantitative evaluation of total P depletion during pedogenesis using this conceptual model. Chronosequence studies show decreasing total P content during pedogenesis, while available P first increases from weathering and then decreases as a result of various soil processes (Crews et al., 1995;Chadwick et al., 1999;Vitousek, 2004;Selmants and Hart, 2010;Walker and Syers, 1976). Most chronosequence studies include soils ranging from slightly weathered Entisols to highly weathered Ultisols and Oxisols. We use measurements from a number of chronosequence studies to estimate the percentage loss of parent material P for soils in different weathering stages. Additionally, we use measurements of total soil P in different vertical soil horizons to estimate the percentage loss of parent material P during weathering on the premise that soils are less weathered with soil depth and that the profiles have been formed from the underlying parent material (St Arnaud et al., 1988).
Several previous studies have estimated soil P in terrestrial ecosystems. The most spatially complete global soil P database, the WISE-ISRIC (World Inventory of Soil Emission Potentials-International Soil Reference and Information Centre) pedon database (Batjes, 2010), contains plant available P data based on the Olsen method (Olsen, 1954) for 1044 soil profiles. However, available P is only a very small fraction of total P in the terrestrial ecosystems (less than 6 %) (Cross and Schlesinger, 1995). Cross and Schlesinger (1995) summarized P data based on the Hedley fractionation method from 88 studies that cover 10 major USDA soil orders and found that the distribution of different P forms in soils is related to the stage of soil development. The Hedley fractionation method fractionates soil P into various inorganic and organic pools by sequentially removing different forms of P with successively stronger extraction agents (Hedley and Stewart, 1982;Tiessen and Cole, 1984) and gives a comprehensive picture of different forms of P in soils. In our recent study, we expanded Cross and Schlesinger (1995) to create a larger Hedley P database, which will help us gain a better understanding of the relationships between the distributions of different P forms and the stages of soil development (Yang and Post, 2011) .
Previous studies have used soil orders in USDA's soil taxonomy as an index for soil development stages and have demonstrated its use in relating soil P and soil development stages (Cross and Schlesinger, 1995;Smeck, 1985;Johnson et al., 2003). Although soil genesis is an important component of the definition of most soil orders (Wilding et al., 1983), the relationship between several soil orders and weathering intensity can be ambiguous. Here we take a more general approach by considering three soil weathering categories, e.g. slightly, intermediately and highly weathered soils. By combining the conceptual model of P transformations with available datasets, we provide a first attempt to estimate soil P status on the global scale for use in initialization of global biogeochemistry-climate models.

Lithology data
The global distribution of the soil parent material is derived by combining two global surficial lithology maps. The fundamental one is the global surficial lithology map of Dürr et al. (2005), which includes 15 rock types plus water and ice. The original map in Dürr et al. (2005) is in vector format with about 8300 polygons (minimal size of polygon about 10 km). We transform it into raster format with the resolution of 0.5 • × 0.5 • . We follow the transformation method of Dürr et al. (2005) by first transforming the vector data into a raster map with a 1 km × 1 km resolution in order to retain as much detail as possible. We then summarize the 1km resolution data at 0.5 • resolution to derive the dominant rock type for each grid cell.
One shortcoming of the Dürr et al. (2005) data is that the sedimentary categories in this dataset do not differentiate between shales and pure sandstones, which may have very different P concentration (Gray and Murphy, 2002). In this study, each 0.5 • grid cell in the sedimentary rock category from Dürr et al. (2005) is further assigned as either shale or sandstone by overlaying the Amiotte Suchet et al. (2003) map, which consists of global distribution of 6 main rock types (shales and sandstones as separate categories) at 0.5 • resolution.

P concentration in parent material of soils
P concentration (in ppm) of each lithology unit is assigned values from Hartmann et al. (2012), which is based on a literature review of typical rock P concentration and composition of rock types per lithology classes (Table 1). P concentration of loess in Hartmann et al. (2012) was based on Chinese loess. We assign P concentration of glacial loess in North America and other regions based on measurements from previous studies (Runge et al., 1974).

Quantification of soil P transformations
The Walker and Syers conceptual model of P transformations during pedogenesis provides a useful tool to link soil development stages and soil P amount and forms (Johnson et al., 2003;Smeck, 1985). Quantification of this conceptual model needs two kinds of information. The first is percentage loss of parent material P for soils in different weathering stages, for which we rely on soil P data from chronosequence studies and soil profile measurements. Since soils are open systems that not only gain or lose mass but also volume as well, volumetric soil strain needs to be considered in order to calculate soil P content to certain depth (Brimhall and Dietrich, 1987;Brimhall et al., 1992;Chadwick et al., 1990). Volumetric strain is the change in volume of mineral mate-rial divided by the original volume from which it is derived. This is 0 if there is no change in volume, positive if there is an increase in volume (through reduction in bulk density for example), or negative if there is a decrease in volume (loss trough dissolution for example). We group 9 of 12 USDA soil order into three soil categories reflecting the relative weathering stages of soils. Entisol and Inceptisol orders are slightly weathered soils. Aridsol, Vertisol, Mollisol, and Alfisol comprise intermediately weathered soils. Spodosol, Ultisol, and Oxisol represent strongly weathered soils. Gelisol, Histosol, Andisol are not considered in this study since they are considered as special soil orders that have unique genetic properties unrelated to the stage of soil development. The second kind of information is the fraction of different forms of P for each USDA soil order, which are derived based on Hedley P data that we compiled from the literature (Yang and Post, 2011).

The depletion of total P during soil development
We use both measurements of total P from chronosequence studies and vertical sequence of soils to estimate the loss of total P (See supplementary, Table S1). The "pedogenic index", was originally used to quantitatively express the changes that occurred in soils during their development (Santos et al., 1986). Here we apply this idea by introducing the "phosphorus pedogenic depletion index" (PPDI) to quantitatively describe the cumulative total P loss during soil development with mass balance considered (Eq. 1).
TP is the current total soil P content (in kg P ha −1 or g P m −2 ) and TP 0 is the total P content (in kg P ha −1 or g P m −2 ) of the least weathered soil in the chronosequence or the total P in the unweathered horizon for the soil profile. ε is the volumetric soil strain (Brimhall and Dietrich, 1987;Brimhall et al., 1992) (see Sect. 3.3). We found eight chronosequence studies (Table S1) in the literature that provide measurement of total P and soil order at each site in the sequences, allowing quantification of the loss of parent material P for soils with different weathering intensity along the chronosequence. The humid tropical chronosequence in Hawaii, in Chadwick et al. (1999), encompasses soils ranging from slightly weathered soils (Entisols and Inceptsols) to highly weathered soils (Ultisols and Oxisols). The semiarid chronosequence in northern Arizona covers soil orders from Entisols to Mollisols to Alfisols (Selmants and Hart, 2010). Five chronosequences (Parfitt et al., 2005;Richardson et al., 2004;Lichter, 1998;Eger et al., 2011) provided the parent material P and total P for podzolised soils (Spodosols). One desert soil chronosequence in southern New Mexico provided the parent material P and total P for Aridsols (Lajtha and Schlesinger, 1988). In addition, we utilize the measurements of total P along the vertical sequence of soil horizons to determine PPDI for Mollisols and Alfisols on the assumption that the parent material from which the soils developed is the material found at the deepest part of the profile (St Arnaud et al., 1988). Most of the chronosequence and soil profile studies provided the estimate of PPDI. We calculate PPDI for the remaining chronosequence studies and soil profile measurements using Eq. 1 (P content was either provided in literature cited or calculated by multiplying P concentration (in ppm) with bulk density and the depth sampled, see Table  S1 for detail). Using these chronosequence and soil profile measurements in Table S1, we derive the average PPDI for slightly and intermediately weathered soil categories (Table   2). We calculate PPDI for highly weathered Spodosols, Ultisols, and Oxisols individually since the distinction of PPDI between Ultisols and Oxisols is crucial to capture the spatial variation of soil weathering stages and total P in Amazonia, where P limitation is pronounced.

Hedley fractionation data
The Hedley sequential fractionation method (Hedley and Stewart, 1982;Tiessen and Moir, 1993), which first removes labile inorganic and organic P and then the more stable inorganic and organic P using sequentially stronger extracting agents (first anion exchange resin, followed by 0.5 M NaHCO 3 , 0.1 M NaOH and 1 M HCl), has gained considerable attention as a useful tool to examine different forms of soil P and provides a comprehensive assessment of available P in soils (Johnson et al., 2003). We expand an earlier study (Cross and Schlesinger, 1995) that summarizes Hedley P data and created a larger Hedley P database that includes 178 published Hedley fractionation measurements (Yang and Post, 2011). We summarized the compiled literature data by soil order for fractions of total P in different P forms (see supplementary, Table S2).

Distributions of different P forms in soils
The spatial distribution of soil P forms is developed in three main steps (Fig. 1). First, we combine the parent material map with the rock P concentration database to generate a map of parent material P concentration (see Sect. 3.1). Second, we derive the map of total P content in top 50 cm soils by overlaying the map of soil order on the map of parent material P, and applying PPDI and soil strain (see below) using Eq. 2.
Where TP S (g P m −2 ) is total P in the top 50 cm soil, D is the soil depth (50 cm), ρ p (g cm −3 ) is the bulk density of parent material (see Table 1), C P (ppm) is parent material P concentration, and ε is the volumetric soil strain. Lastly, we derive the maps of different forms of P in soils by applying the relationship between soil order and the fractions of total P held in different P forms derived based on literature Hedley data (see Sect. 3.2.2). The global distribution of soil orders is obtained from the USDA website (http://soils.usda.gov/use/worldsoils/ mapindex/order.html). For Latin America, we use the soil order map based on the Soil and Terrain database for Latin America and the Caribbean (SOTERLAC, http://www.isric.org/. We convert the (Food and Agriculture Organization) FAO90 soil code used in the SOTERLAC database to USDA soil taxonomy following Quesada (2011). The SOTERLAC database was used to replace the USDA soil order map for Latin America because the USDA map is based on a reclassification of the FAO-UNESCO Soil Map of the World combined with a soil climate map, which has a considerable bias towards the dominance of Oxisols (Ferralsols in FAO map) in Amazonia (Richter and Babbar, 1991). We collected strain values of top 50 cm soils from the literature and group them in the three soil weathering categories (Table S3 in supplementary and Table 3). Soil strain ε is calculated based on the assumption that a relatively immobile element (usually Ti or Zr) is conserved during pedogenesis (Brimhall and Dietrich, 1987;Brimhall et al., 1992;Chadwick et al., 1990). Positive strain values represent soil dilation due to biological activity, organic matter accumulation and other physical processes. Negative strain values represent soil collapse through weathering and leaching. Generally soil is dilated (positive strain) at earlier stages of soil development and then collapsed (negative strain) for highly weathered soils. It is worth mentioning that soil strain values not only depend on the extent of soil weathering, but also on soil parent material (Schroeder and West, 2005;Merkli et al., 2009). In particular, calcareous minerals are very soluble and the strong leaching of carbonate can lead to strongly negative strains even in slightly weathered Entisols and Inceptisols (Merkli et al., 2009). Soil strain values are usually positive for Entisols and Inceptisols developed from silicate minerals (Table S3). For highly weathered Ultisols and Oxisols in tropical soils, surface soil strain ranges from -0.25 to 0 due to a progressive upward re-expansion, which may be caused by root regrowth and burrowing animals (Colin et al., 1992).

Uncertainty estimates
The uncertainty of our soil P estimates depends greatly on the uncertainty of PPDI, Hedley fractions of different P forms, and soil strain. The mean and standard deviation of PPDI (Table 1), Hedley fractions (Table S2) and soil strain (Table 3) were used here to quantify the uncertainties of estimated soil P. Here we use coefficient of variation (%, standard derivation divided by mean then multiplied by 100) to describe the uncertainty in our estimate. Uncertainties are propagated in quadrature (Eq. 3). The assumption here is that the variables  involved are uncorrelated and normally distributed.
where σ is standard deviation and x 1, x 2, x 3 are PPDI, Hedley faction and soil strain respectively.

The global distribution of P in parent material
The parent material of a soil determines the initial amount of P, which is released by weathering and then either retained in soils through a number of chemical and biological processes or lost through leaching and erosion. Globally parent material P ranges from less than 300 ppm to more than 1300 ppm (Fig. 2). As shown in Fig. 2, parent material P is highest in the west coastal regions of North and South America, as well as in some parts of south Asia and Africa, where volcanic rocks are prevalent. Parent material P is also high in central parts of North America, northern Europe and Asia, where shale rocks are common. In the Sahara region of Africa, southern Africa, central Asia, northwestern Russia, and central part of Australia, where sand or sandstones occur, parent material P is relatively low, around 350 ppm. We found no latitudinal pattern of parent material P distribution.

The global distribution of different P forms in soils
Using the mean PPDI and Hedley fraction values, we generated the spatial distribution of different forms of P and total P in soils (Fig. 3). Total P is lowest in tropical regions, where soils have gone through millions of years of soil development and have lost most of their original P through leaching or erosion. The warm and humid climate in tropical regions enhances the weathering of parent material and leaching of P from soils. Total soil P is high in the slightly weathered Inceptisols and Entisols along the east coast of North America and South America. Total P in the top half meter of soils ranges from less than 100 g P m −2 in highly weathered Oxisols to around 500 g P m −2 in Inceptisols and Entisols. Global total P to 50 cm soil depth on the global scale is 40.6 Pg P. Estimates of total soil P on the global scale in the literature vary greatly, from 30.6 Pg P (Wang et al., 2010) and 40 Pg P (Smil, 2000) to 200 Pg P (Jahnke, 1992). Our estimate of total soil P is near the lower end of previous estimates. Labile inorganic P represents the most readily available form of P for plants, but generally is a small proportion of total P. Labile inorganic P is generally regarded as biologically available. It is defined through the sequential extraction procedure in the Hedley method (extractions with exchange resin followed by 0.5 M NaHCO3, each for 16 h). Because of the length of time required for this extraction, from modeling perspective, labile inorganic P from Hedley fractionation method may or may not be treated as plant available P depending on the model structure and time step. It is reasonable to treat labile inorganic P as plant available P in models with daily time steps (Goll et al., 2012;Wang et al., 2010). However if the model has a much shorter time step (hourly or less), it is not appropriate to treat labile inorganic P as plant available P. As shown in Fig. 3, of all forms of P, labile inorganic P has the lowest content in soils. Labile inorganic P is higher in slightly weathered Entisols and Incepitsols, with the average content of around 52.0 g P m −2 , mainly due to the rapid weathering of parent material in these soils. However there is not much variation in labile inorganic P among soils with different weathering stages. Labile inorganic P is not only controlled by parent material and weathering stage, but also biological processes such as plant uptake, microbial uptake, immobilization and mineralization. It has been suggested that in tropical regions, mineralization of organic P is the major source for labile inorganic P (Vitousek, 1984). In addition, fluctuating redox conditions in highly weathered wet tropical soils can lead to the release of labile inorganic P through the reduction of Fe(III) minerals (Chacon et al., 2006;Liptzin and Silver, 2009). Therefore there is no clear relationship between labile inorganic P and soil order. This has been also suggested in previous studies (Cross and Schlesinger, 1995;Zhang et al., 2005). Total labile inorganic P in the top half meter of soil on the global scale is estimated as 3.6 Pg P.
Organic P (Po) shown here contains both easily mineralized Po, which can contribute to plant available P in the short term, and more stable forms of Po, which is usually involved in the long-term transformations of P in soils. Overall, organic P is lower in tropical regions and higher in temperate regions. The low organic P in tropical regions is the result of several processes: loss of P through leaching, adsorption of P on secondary aluminum and iron oxide minerals, and increasing occlusion of P by Al and Fe minerals under low soil pH conditions, all of which increase with soil development. In addition, the warm and humid climate in tropical regions enhances soil organic matter decomposition and mineralization of organic P. While in temperate regions, where most of soils are of intermediate stages, a continuing supply of inorganic P from weathering of parent material and moderate loss of P due to leaching and erosion leads to the ample supply of P for plants and microbes and highest Po in soils. The relatively slow organic matter decomposition in temperate regions also contributes to the accumulation of Po. Globally, Po in the top half meter soil is estimated about 8.6 Pg P. Smil (2000) estimated that Po in soils ranged from 5 to 10 Pg P; so our estimate is near the higher end of their estimates. Our estimate is higher than the estimates of Mackenzie et al. (2002) (5 Pg P), Goll et al. (2012) (5.7 Pg P) and Wang et al. (2010) (4.4 Pg P).
Secondary P represents inorganic P that is held on the surface of soil minerals by various reactions such as anion adsorption. Secondary P can dissolve in soil solution and become available for plants, but at a much slower rate compared to labile inorganic P. Secondary P is higher in slightly weathered soils and lower in highly weathered soils. The low concentration of secondary P in Oxisols, Ultisols, and Spodosols results from leaching losses or occlusion within secondary minerals. The forms of secondary P are highly dependent on soil pH. In acid soils of tropical regions, secondary P is bound by iron and aluminum oxide minerals, while in alkaline soils P is bound with calcium minerals. The total global secondary P in the top half meter soil is about 3.2 Pg P. Our estimate is higher than the estimates of Goll et al. (2012) (1.3 Pg P) and Wang et al. (2010) (1.7 Pg P).
Occluded P refers to P physically encapsulated inside Alor Fe-or Ca-oxide minerals, which makes P unavailable for plant use. Although the proportion of occluded P in total P generally increases with soil development, occluded P is higher along the west coast of North America and South America, mainly because of the higher P content of parent material in these regions and minimal loss of total P. Globally, total occluded P in the top half meter soil is about 12.2 Pg P.
The distribution of apatite P clearly shows the weathering stage of the soils. There is almost no apatite P in strongly weathered Oxisols and Utisols, while apatite P is higher in slightly weathered Entisols and Inceptisols. The high apatite P in Aridsols despite the moderate weathering stage results from dry climate condition and the generation of secondary calcium phosphate in Aridsols, where the high concentrations of calcium carbonate leads to the chemical reaction between released P from primary apatite P and calcium minerals under neutral or alkaline conditions (Cross and Schlesinger, 1995;Lajtha and Schlesinger, 1988). Globally there is about 13.0 Pg P in the form of apatite in the top half meter soils. Zhang et al. (2005) investigated the distributions of soil P in the top 50cm of soil in China based on total P measurements of more than 2400 soil profiles. In particular, they provided the changes of total P for three different weathering stages: slight (Andisols, Entisols, Gelisols, Histosols, Inceptisols), intermediate (Aridsols, Vertisols, Alfisols, Mollisols), and strong (Spodosols, Ultisols, Oxisols). Table 4 compares our mean estimates of average total P for the three different weathering stages in China with the Zhang et al. (2005) results. Our estimate of total P for slightly and highly weathered soils in China is comparable to the Zhang et al. (2005) estimate. However our estimates of total soil P for intermediately weathered soils in China are lower than Zhang et al. (2005). This may be due to the fact that soil profiles used in Zhang et al. (2005) not only cover natural ecosystems but also agriculture land where P fertilizer has been applied. Our estimate of the total national soil P in the surface half meter is 2.8 Pg P, close to the 3.5 Gt P reported in Zhang et al. (2005).

Comparison with the absolute value of Hedley data from the global collection of pedon measurements
We converted our estimated total P content (g P m −2 ) to total P concentration (ppm) using a global map of soil bulk density (Global Soil Data Task Group 2000) and 50 cm soil depth and then averaged them for USDA soil orders. We compared our estimated mean total P concentration with total P concentration based on Hedley P measurements from pedon data assembled from the literature and averaged for each USDA soil order (Yang and Post, 2011). As shown in Fig. 4, overall our estimate of total P is higher than that based on the Hedley database, especially for highly weathered Spodosols and Ultisols. Our estimated total soil P for intermediately weathered soils is very close to that from data averaged by soil orders. However, there are large discrepancies between our estimate of total soil P and that from the literature for slightly weathered Inceptisols and highly weathered Spodosols and Ultisols. The discrepancy could be due to the fact that our estimate accounts for varying P content in the distribution of parent material, while total P from literature data are based on limited number of observations, which do not account for the spatial heterogeneity of total P on the global scale. The discrepancy could also be the result of our underestimate of P loss from leaching during development of these soils. Unfortunately, our estimate of PPDI for slightly and highly weathered soils is poorly constrained across lithology types and climate, relying on averaging PPDI from a limited number of chronosequence studies and soil profile measurements. In spite of the discrepancy for slightly weathered soils, the general agreement of total P for different soil orders between our estimate and the measurements based on the Hedley P database indicates that the approach we take here is reasonable.   Fig. 4. The comparison of estimated total P with field measurement of total P concentration (mean with one standard error, based on Hedley fractionation data from the literature (Yang and Post, 2011)) for soil orders grouped into three weathering categories.

Uncertainties and limitations
Using the method described in Sect. 3.4, we are able to quantify the uncertainties of our estimated soil P due to the uncertainties from PPDI, Hedley fractions and soil strain. The uncertainty of total soil P is about 17 % for slightly weathered soils, 65 % for intermediately weathered soils and 68 % for highly weathered soils. The uncertainties for different forms of soil P are unavoidably large (most of them around or larger than 50 %, Fig. 5) because of the large variations in PPDI (Table 2), Hedley fractions (Table S2) and soil strain ( Table  3). As further data provide better estimates of PPDI, Hedley fractions and soil strain, these uncertainties will be reduced and we can increase the confidence in estimates of soil P. In addition to the uncertainty of PPDI and Hedley fractions, there are other sources of uncertainty that contribute to the uncertainty of our estimate of soil P. For example, there is uncertainty in the lithology map we used in this study, as discussed in detail by Dürr et al. (2005). In brief, the map area for ice, major water bodies and basic volcanic rocks in most regions has an uncertainty of less than 10 %; the mapping for most very common (>15 % global occurrence) and common rock types (5-15 % global occurrence) has an uncertainty of between 10 and 20 %, except for the differentiation of intermediate rock classes such as complex lithology, mixed sedimentary, Precambrian basement, which has an uncertainty of 20 to 50 %. The mapping for evaporites has a large uncertainty -greater than 50 %. While the lithology map used here is suitable for global-scale studies, uncertainty increases greatly when used at regional and local scales (Dürr et al., 2005). For example, the surficial lithology map we used here failed to capture the loess cover on nearly the entire Columbia Plateau. If a much more detailed world surficial lithology map with higher resolution becomes available, this source of uncertainty can be reduced.
The other source of uncertainty comes from the estimate of P content for each lithology category. Since each lithology category contains various types of rocks that have different P concentrations, the P content of each lithology cat-egory depends not only on P concentrations of rock types, but also on the rock class composition in each category. In this study, P concentration of each lithology category is from Hartmann et al. (2012), which is based on the geochemical composition of typical rocks and the assumptions about the relative proportions of each rock class in each lithology category. The P concentration of typical rock types used in Hartmann et al. (2012) is consistent with previous studies (Anderson, 1988) and earlier rock P database (Gray and Murphy, 2002). However, the assumptions regarding rock-class composition per lithological classes are based on reviews of limited sources and knowledge from regional literature and can have a large uncertainty. For example, it was assumed in Hartmann et al. (2012) that mixed sedimentary consolidated rocks (SM in Table 1) are comprised of 15 % carbonates, 60 % shales, and 25 % sandstone. Dürr et al. (2005), however, suggested that this lithology class has higher carbonate contents, ranging from 30 to 70 %. Clearly significant improvements can be achieved if local or regional geochemical compositions per lithology category are combined with regional lithology maps of improved resolution.
Here we take a very general approach by grouping soil orders into three weathering categories due to the lack of measurement data. It will be ideal if we can conduct this study at a more detailed level, for example the great group level or less preferably at the suborder level. However, P data is presently sparse at either suborder or great group levels, precluding us from pursuing this study at these more detailed levels. Nonetheless, on the global scale this general approach can provide an overall view of weathering status of soils and by using best available data, our study will be useful for estimating the global distribution of available soil P needed for initialization of coupled climate-biogeochemical models.
The sampling depth in our PPDI data (Table S1) ranges from 10 to 50 cm. We applied the averaged PPDI (Table  2) to a 50 cm depth assuming 50 cm depth is appropriate to represent the active mineral soil layer that have undergone approximately the same extent of weathering. Although this is a reasonable place to start, in reality, variations in pH,   5. Estimated soil P content (g P m −2 ) to 50 cm depth of soil for soil orders in three soil weathering categories. The uncertainty of these estimates is indicated by the standard deviations, which is calculated based on the uncertainty of PPDI (Table 2), the uncertainty of Hedley fractions (Table S2), and the uncertainty of soil strain (Table 3). PPDI (phosphorus pedogenic depletion index) is an estimate of the accumulated loss of parent material P during pedogenesis.
redox, soil mineralogy and plant-soil interactions can cause PPDI to change with depth in soil profiles. For example, most Spodosols exist in acidic environments where weathering happens rapidly and P could be leached to deeper soil layers along with Fe and Al through the action of dissolved organic carbon (DOC). Thus we could overestimate PPDI for Spodosols by applying PPDI from top 10-15 cm mineral soil to 50cm depth. Another example is the redistribution of soil P from deep horizon to surface through biological activity, whereas P taken up by plant roots accumulates in surface horizons through litter decomposition and mineralization. This process may lead to the underestimation of PPDI when we apply PPDI from top 10-15 cm mineral soil to 50 cm depth.
In this study we are not considering the effects of dust transfer on P status of soils, which could be fairly important in some parts of the world. For example, in Hawaii (Chadwick et al., 1999) and the Amazon region (Swap et al., 1992), transported dust is reported to be an important source of P in these highly weathered soils. Dust from arid central Asia is the dominant source of P input (relative to P release from parent material) on the highly weathered site from the Hawaii chronosequence. We are probably underestimating PPDI of highly weathered soils by using data from Hawaii chronosequence. However we do not anticipate dust will significantly impact the overall pattern of total soil P for most soils because the importance of dust input varies depending on the dust deposition rate and the amount of total soil P. Porder and Hilley (2010) showed that addition of dust did not substantially alter the trends in their estimate of fraction of parent material P remaining in soil. Furthermore little is known about the large variation in dust fluxes over the timescale of soil development, making it difficult to incorporate a definitive dust effect here. For any global models that use our estimate for initialization, it will be important to consider dust input during model spinup.
We have not considered the impacts of human activity on soil P in terrestrial ecosystems. Studies have suggested that fertilizer application accounts for 42 % of total P inputs to agricultural ecosystems, or 25 % of total plant P uptake in fertilized soils (Wang et al., 2010). Therefore, our analysis underestimates soil P content in agricultural regions. The estimates for agricultural soil P will be improved if the maps derived here were used in conjunction with process-based models and global fertilizer data (Bouwman et al., 2009;MacDonald et al., 2011). We are also not considering the effects of tectonic uplift and erosion, specifically, since we assume that on geological timescales they are balanced with each other. This assumed geological equilibrium may be inaccurate for some soils (Porder et al., 2007).
Here we are not explicitly modeling various P processes in terrestrial ecosystems over geological timescales to provide initial conditions for global models, as it requires unrealistic computational time. We acknowledge that although our data-based approach can provide a reasonable initial condition for global biogeochemical models, they may not be accurate at local or regional scales due to lack of consideration of specific processes at those scales. For example, biological processes including plant uptake, immobilization, biological and biochemical mineralization, biologically enhanced weathering, etc., can influence soil P distribution. These processes are not explicitly accounted for in this study but are considered in most biogeochemical models that include the P cycle (Wang et al., 2010;Goll et al., 2012;Zhang et al., 2011). We recommend the maps provided here to be used in conjunction with process-based models as a reasonable starting point for describing soil P distribution and availability at regional and global scales and to improve our understanding of P dynamics in terrestrial ecosystems.

Conclusions and future research
We derived the global distribution of different forms of P based on surface lithology maps of the Earth, distribution of soil development stages, fraction of parent material P remaining for different soil weathering stages using chronosequence studies, and the distribution of P in different forms for each soil order based on the analysis of Hedley P data from a global literature collection of profile data. The general agreement between our estimated total P and measured total P indicates that initial P concentration in parent material and geochemical processes such as weathering of parent material and the loss of P via leaching are dominant in regulating soil phosphorus cycles in the long term. Our estimates of different forms of P are derived from estimated total P and the fraction of total P for each P form based on the Hedley P dataset. The amounts of P, to 50 cm soil depth, in soil labile, organic, occluded, and secondary pools are 3.6 ± 3, 8.6 ± 6, 12.2 ± 8, and 3.2 ± 2 Pg P respectively. The amount of P in soil mineral particles to the same soil depth is estimated at 13.0 ± 8 Pg P and global soil total P is 40.6 ± 18 Pg P.
Our study attempts to estimate the global distribution of different forms of P in soils for the initialization of globalscale biogeochemical models that incorporate P cycling. We find that sufficient spatial data and soil process measurements are available to develop a spatially explicit map of the magnitude and distribution of different P forms in soils. These estimates are consistent with other regional and global estimates of soil P amounts and distribution. This spatial and process specific detail is important for global terrestrial biogeochemical model initialization. The analysis reveals significant uncertainties. The uncertainty of total soil P is about 17 % for slightly weathered soils, 65 % for intermediately weathered soils and 68 % for highly weathered soils. The uncertainties for different forms of soil P are unavoidably large (most of them around or larger than 50 %) because of the large variations in PPDI, Hedley fractions and soil strain. The reduction of uncertainties will require improved spatial data resources, additional sampling in under-represented soil types, and further improvements in process level understanding. In particular our analysis used a global lithology map that aggregated parent material into 15 types. Actual global lithology distribution is much more heterogeneous (Dürr et al., 2005) and so analysis at finer spatial resolution will require a more detailed lithology description. There is uncertainty in the estimate of P content of each lithology type used in this study, mainly due to our limited understanding of the rock composition of each lithology type. The improvement of this requires a better understanding of local or regional geochemical compositions per lithology category and regional lithology maps with improved resolution. Our estimates of total P loss during soil development are based on limited number of chronosequence studies and soil profile measurements. Improvement of this estimate requires more data that can provide a better understanding of the relationships between soil forming factors and pedogenic P depletion. Although our general approach of using three soil weathering categories is useful on the global scale, it may be too coarse at local or regional scales. Additional data at suborder or great group level is needed to refine our estimates.
In spite of the large uncertainty of our estimates, we believe a global soil P assessment for global biogeochemical models is needed, using the best available data and approaches. This is an initial attempt to tackle soil P estimates on the global scale by combining the most complete and comprehensive data on lithology, chronosequence studies, and soil P measurement data with our understanding of P transformations during pedogenesis. The estimates in this study can be refined and improved when more data become available. To improve estimates, research on understanding the mechanisms and processes controlling P weathering, transformation and loss need to be expanded with emphasis on understanding the relative contributions of geochemical and biological processes leading to different forms of P in soils. This would also improve our capability of modeling P dynamics in soils and predicting the role of P in terrestrial productivity. In addition, field measurements of different forms of P (particularly using the Hedley P fractionation method) should be expanded to additional climate regimes, parent materials and soil development stages in order to provide critical information on the contemporary P stocks in soils and determine the climatic, parent material, and soil age controls over the distributions of different forms of P in soils.