scientific **data**

OPEN

DATA DESCRIPTOR

<figure><img src="image_1.png" /></figure>

Check for updates

# A global dataset on phosphorus in agricultural soils

Bruno Ringeval ID<fnref n="1" />✉, Josephine Demay<fnref n="1" />, Daniel S. Goll², Xianjin He ID², Ying-Ping Wang ID³, Enqing Hou ID⁴, Sarah Matej⁵, Karl-Heinz Erb ID⁵, Rong Wang ID⁶, Laurent Augusto ID⁷, Fei Lun⁷, Thomas Nesme ID<fnref n="1" />, Pasquale Borrelli⁸,⁹, Julian Helfenstein¹⁰, Richard W. McDowell ID¹¹,¹², Peter Pletnyakov ID¹² & Sylvain Pellerin<fnref n="1" />

Numerous drivers such as farming practices, erosion, land-use change, and soil biogeochemical background, determine the global spatial distribution of phosphorus (P) in agricultural soils. Here, we revised an approach published earlier (called here GPASOIL-v0), in which several global datasets describing these drivers were combined with a process model for soil P dynamics to reconstruct the past and current distribution of P in cropland and grassland soils. The objective of the present update, called GPASOIL-v1, is to incorporate recent advances in process understanding about soil inorganic P dynamics, in datasets to describe the different drivers, and in regional soil P measurements for benchmarking. We trace the impact of the update on the reconstructed soil P. After the update we estimate a global averaged inorganic labile P of 187 kgP ha⁻<fnref n="1" /> for cropland and 91 kgP ha⁻<fnref n="1" /> for grassland in 2018 for the top 0–0.3 m soil layer, but these values are sensitive to the mineralization rates chosen for the organic P pools. Uncertainty in the driver estimates lead to coefficients of variation of 0.22 and 0.54 for cropland and grassland, respectively. This work makes the methods for simulating the agricultural soil P maps more transparent and reproducible than previous estimates, and increases the confidence in the new estimates, while the evaluation against regional dataset still suggests rooms for further improvement.

## Background & Summary

Soil phosphorus (P) availability limits crop yields in many regions of the World<fnref n="1" />,², while anthropogenic application of P in excess triggers aquatic eutrophication in other regions³. While P in unmanaged soils can be predicted to some degree from local soil properties and climate⁴⁻⁶, this approach cannot be used for agricultural systems because humans have altered the P cycle in these ecosystems. E.g. the change in soil P input/output related to farming practices has modified the spatial distribution of P in soils⁷⁻⁹ with long-lasting legacy effects of past management¹⁰. Nevertheless, an analysis of the relationships between crop productivity and soil P availability or between P losses and soil P on regional to global scale requires a spatially explicit agricultural soil P dataset. At the global scale, such a dataset has only recently been developed by Ringeval *et al.*¹¹ and Zhang *et al.*¹².

Both Ringeval *et al.*¹¹ and Zhang *et al.*¹² combined datasets describing the different drivers of the global spatial distribution of soil P with a process-based model of soil P dynamics to reconstruct the past and present distributions of P in agricultural soils. Hereafter, the term “drivers” encompasses: variables related to soil P input and output (farming practices, deposition, erosion, sludge), land-use change, variables that have an effect on soil P dynamics (soil properties, climate), and biogeochemical background (i.e. the initial soil P content at the time

<footnote n="1">¹ISPA, Bordeaux Sciences Agro, INRAE, 33140, Villenave d'Ornon, France. ²Université Paris Saclay, CEA-CNRS-UVSQ, LSCE/IPSL, Gif-sur-Yvette, France. ³CSIRO Oceans and Atmosphere, Aspendale, Vic., Australia. ⁴Key Laboratory of Vegetation Restoration and Management of Degraded Ecosystems, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, China. ⁵Institute of Social Ecology, University of Natural Resources and Life Sciences, Vienna, Austria. ⁶Department of Environmental Science and Engineering, Fudan University, Shanghai, 200438, China. ⁷College of Land Science and Technology, China Agricultural University, Beijing, 100193, China. ⁸Department of Science, Roma Tre University, 00146, Rome, Italy. ⁹Department of Biological Environment, Kangwon National University, Chuncheon, 24341, Republic of Korea. ¹⁰Soil Geography and Landscape Group, University of Wageningen, Wageningen, 6700AA, The Netherlands. ¹¹AgResearch, Lincoln Science Centre, Private Bag 4749, Christchurch, 8140, New Zealand. ¹²Faculty of Agriculture and Life Sciences, Lincoln University, Lincoln, PO Box 84, 7647, Christchurch, New Zealand. ✉e-mail: bruno.ringeval@inrae.fr</footnote>
of conversion to cropland). The soil P dynamics model was used to reconstruct how the different soil P pools characterized by different bio-availabilities evolved in time from the biogeochemical background following i) soil P input/output directly (e.g. farming practices) or indirectly (land-use change) derived from some drivers and ii) as function of some drivers that modulate the fluxes between soil P pools (e.g. climate).

Refs. ¹¹,¹² used similar datasets to describe the main drivers. In particular, both used the estimates by Bouwman *et al.*¹³ for soil P input/output corresponding to farming practices, and unmanaged soil P data by Yang *et al.*⁴ to approach the biogeochemical background. These dataset have some limitations: e.g. the soil P input/output in grassland provided by Bouwman *et al.*¹³ was poorly described with the plant uptake forced to be equal to a constant fraction of soil P input. In the work of Yang *et al.*⁴ how the total soil P is distributed among different soil pools was only a function of soil order while it is now known that the soil order is a poor predictor¹⁴. The representation of soil P pools and fluxes between these pools differed between Ringeval *et al.*¹¹ and Zhang *et al.*¹², reflecting the gaps in our understanding of the soil P dynamics and data to constrain process-based models. Zhang *et al.*¹² simulated only two soil pools (one labile and one stable pool) with poor representation of soil P dynamics, whereas Ringeval *et al.*¹¹ used a 6-pool P model based on Hedley fractionation method. Despite a more apparent mechanistic representation of soil P transformation, the model used by Ringeval *et al.*¹¹ was based on a very limited knowledge about soil P dynamics: e.g. the soil P buffering capacity was represented as function of soil order only and was poorly constrained by measurements¹⁵.

Because of the above-mentioned limitations, soil P maps simulated by both Ringeval *et al.*¹¹ and Zhang *et al.*¹² had low confidence and none of these studies really investigated the realism of the spatial distribution of the soil P simulated: Zhang *et al.*¹² focused on the P plant uptake while Ringeval *et al.*¹¹ assessed the contribution of each driver to the spatial variation of soil P. The simulated soil P maps were not evaluated and this was in partly explained by missing large-scale datasets of P measurements at that time. Since the initial release in 2017, the knowledge about soil P dynamics has improved, datasets representing some drivers of soil P have been updated, and data on measured soil labile P in agricultural soils started to be available, offering the opportunities to generate more reliable global soil P maps. In particular, the comparison between different measurements (isotopic exchange kinetics and Hedley fractionation) improves our understanding of the forms extracted¹⁶,¹⁷ and our ability to represent the fluxes between these pools¹⁸. Also, parameterizations of exchanges between inorganic P pools have been improved by considering the effects of soil properties on these fluxes¹⁸. Dataset used to describe some drivers has been also improved. In particular, both total soil P and its fractions within different pools have been improved by including more predictors and more sites to train the machine learning algorithms in these studies, improving our understanding of variables driving soil P distributions in unmanaged soils⁶. Finally, regional datasets on measured soil available P based on soil monitoring networks are now available (e.g. LUCAS database¹⁹).

Here we updated the approach described in Ringeval *et al.*¹¹ with the aim to improve the global soil P dataset in five ways: i) use new theoretical understanding of inorganic soil P dynamics, ii) use new and better constrained datasets to describe the different drivers, iii) benchmark the new soil P maps with available global and regional estimates, iv) extend the period of simulations over time (simulations ends in 2005 in Ringeval *et al.*¹¹), v) clarify the different assumptions used throughout all steps of our approach and v) make the whole approach (i.e. scripts to download the datasets, scripts to generate the input from these datasets, soil P dynamic model) available to the community. Points iv) and v) would make further updates easier.

The update of Ringeval *et al.*¹¹ described in this article provides the simulated size of different soil P pools for the period 1900–2018, with a spatial resolution of 0.5° for both croplands and grasslands. We also provided estimates of the uncertainty related to the description of the different drivers. Spatially gridded estimates of soil P pools will be useful to quantify the current limitation of agricultural production by P and to inform assessments of strategies aiming at reduced fertilizer input by making better use of P present in soils.

## Method

**General approach.** Similar to Ringeval *et al.*¹¹, the approach used here modelled the global distribution of soil P in croplands and grasslands, with a 0.5° × 0.5° spatial resolution. No vertical discretization was considered as the approach considers only one soil layer, which is the top 0–0.3 m soil layer. For each grid-cell, soil P was distinguished into four land-cover type fractions: (i) cropland (crop), (ii) grassland (grass) defined as managed pasture + rangeland in Chini *et al.*²⁰, (iii) non-agricultural vegetation (nonagri), and (iv) urban land (urban). While the soil P of the last two land-cover fractions were not explicitly simulated (see below), their consideration was necessary to account for the effect of the land-use and cover change on the cropland/grassland soil P. Land-use and cover change and soil P input/output were prescribed at yearly time-step while soil P dynamics work at daily time-step. Simulations cover the 1900–2018 period.

Our modelling approach combined several global datasets describing the drivers of agricultural soil P with a soil P dynamics model to simulate the temporal evolution of the agricultural soil P from the beginning of the 20th century to the present time period. Drivers are named with four capital letters hereafter. Starting from the soil biogeochemical background used as initial conditions for the year 1900 (BIOG), our modelling approach compute the temporal evolution of soil P pools for cropland and grassland by, i) computing annually the effect of land-use and land-cover change (LUCC), ii) considering the annual soil P input/output corresponding to different drivers (chemical fertilizer, manure, plant uptake and plant residues, encompassed under the term “farming practices”, FARM; input from atmospheric deposition, DEPO; P losses through water erosion, LOSS, P input from sludges, SLUD), iii) simulating at daily time-step the soil P dynamics (i.e. how the different soil P pools exchanges together). Fluxes between soil P forms depends on two drivers: near-surface air temperature and soil temperature/water content (CLIM), and soil properties (SPRO) In total, eight drivers were considered (Tables 1–8). The datasets used to represent the different drivers resulted from the combination of different measurements (satellites, on site measurements, etc.) and/or model simulations. All data sets were regridded to a half-degree resolution.
<table><caption>Table 1. FARM driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>FARM (farming practices)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>Soil P input/output related to farming practices</td></tr><tr><td>Corresponding variables in our approach</td><td>Soil P input/output corresponding to: chemical fertilizer (F<sup>out→i</sup><sub>chem</sub>), manure reaching cropland/grassland soils (F<sup>out→m</sup><sub>manu</sub> with m in {i-lab, o-lab, o-sta}), crop/grass residues (F<sup>out→m</sup><sub>res</sub> with m in {i-lab, o-lab, o-sta}), and plant uptake (F<sup>i→lab</sup><sub>upta</sub>), Here, 'residues' refer to plant biomass that remains on/within the soil after the harvest (for cropland) or harvest/grazing (for grassland). It includes root biomass.</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>Chemical fertilizer and manure were taken from Bouwman et al.<sup>13</sup>. Residues and plant uptake were derived from P harvest/ withdrawal provided by Bouwman et al.<sup>13</sup> following basic assumptions.</td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td>The approach used in Bouwman et al.<sup>13</sup> (based on IMAGE<sup>67</sup>) excluded ~20% of cropland area (due to some boolean treatment about land-use). This lead to some inconsistency with LUCC in GPASOIL-v0. In Bouwman et al.<sup>13</sup>, the treatment in grassland was simple: P withdrawal was estimated as ~90% of the total applied fertilizer. In GPASOIL-v0, we used few global numbers, without crop distinction, to derive uptake and residues from withdrawal/harvest (see Supp.Inf of Ringeval et al.<sup>11</sup>).</td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>For both cropland and grassland, P in manure reaching the soil was derived from a combination between N in manure reaching grassland soil given by Xu et al.<sup>51</sup>, a global P:N ratio (0.20), and P manure produced at country-scale based on livestock category population and manure production rate<sup>9</sup>. Grassland: P in chemical fertilizer was derived from N in chemical fertilizer applied to grassland given by Xu et al.<sup>31</sup> and a global P:N ratio (0.22). P plant uptake was derived from the actual NPP (including above and belowground NPP) in grassland and some spatially constant parameters (P concentration, etc.). P in residues was derived from P plant uptake and (1-GI) with GI: grazing intensity given by Kastner et al.<sup>41</sup> and corresponding to the ratio of the Human appropriation of net primary production through harvest and the actual NPP. The spatial distribution of NPP provided by Kastner et al.<sup>41</sup> was used but global average NPP was made varying between the values of Kastner et al.<sup>41</sup> and Wang et al.<sup>42</sup>. Cropland: P in chemical fertilizer was provided by Lu and Tian<sup>34</sup>. The latter ref.<sup>34</sup> assumed that all chemical fertilizer is applied to cropland and here, we corrected the chemical fertilizer applied to cropland by subtracting the chemical fertilizer applied on grassland described above. P uptake and P in residues were computed as functions of yield and crop specific parameters. Crop-dependent, spatially explicit and temporal varying yields are get by combining spatially-explicit yield per crop for year 2000 by Monfreda et al.<sup>39</sup> and country-averaged yield (but temporal-varying) given by FAOSTAT. The yield of FAOSTAT at country-scale is extended for years before 1961 by scaling the 1961 yield with the country-sale human population provided by HYDE 3.2<sup>23</sup>. Yield is then combined with different crop dependent parameters (concentration of P in different organs - root, harvest, aboveground biomass excluding the harvest, root:shoot ratio and harvest index) and fraction of total residue remaining on/within soil (equal to 1/2 or 0, following Smil<sup>44</sup>).</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>Xu et al.<sup>31</sup> provided for grassland, chemical N fertilizer for years between 1961 and 2016, and N in manure for years between 1860 and 2016, both in kgN.km<sup>-2</sup>.yr<sup>-1</sup>, and at half-degree spatial resolution<sup>68</sup> (https://doi.org/10.1594/PANGAEA.892940). Demay et al.<sup>9</sup> provided P manure produced at country-scale per livestock category for years between 1950 and 2017 (<sup>69</sup>https://doi.org/10.57745/ LEPJCS). Kastner et al.<sup>41</sup> provided NPP (in gC m<sup>-2</sup> yr<sup>-1</sup>) and GI (no unit) for grassland for some years (1910, 1930, 1950, 1960, 1970, 1980, 1990, 2000, 2010) at 0.083° spatial resolution(<sup>70</sup>https://doi.org/10.5281/zenodo.7313791). Lu and Tian<sup>34</sup> provided chemical P fertilizer applied on cropland for any years between 1900 and 2013, in gP m<sup>-2</sup> yr<sup>-1</sup>, at half-degree spatial resolution, (<sup>71</sup>https://doi.org/10.1594/PANGAEA.863323). Zhang et al.<sup>30</sup> provided N in manure applied on cropland for any years between 1860 and 2014, in kgN.km<sup>-2</sup>.yr<sup>-1</sup> at half-degree spatial resolution (<sup>72</sup>https://doi.org/10.1594/PANGAEA.871980). Monfreda et al.<sup>39</sup> provided crop-dependent and spatially explicit (half-degree resolution) yield for years ~2000 (Earthstat dataset, http://www.earthstat.org/ harvested-area-yield-175-crops/). FAOSTAT provided country-averaged and crop-dependent yield for any years over 1961-2019 (FAOSTAT. License: CC BY-NC-SA 3.0 IGO. Extracted from: https://www.fao.org/faostat/en/#data/QCL and https://www.fao.org/faostat/en/#definitions Data of Access: 20-01-2023.). Klein Goldewijk<sup>23</sup> provided the country-sale human population (HYDE 3.2<sup>73</sup>, https://doi.org/10.17026/dans-25g-gez3).</td></tr></tbody></table>

<table><caption>Table 2. DEPO driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>DEPO (P in atmospheric deposition)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>Soil P input resulting from deposition of atmospheric P</td></tr><tr><td>Corresponding variables in our approach</td><td>P deposition falling within labile (F<sup>out→i</sup><sub>depo</sub>) and apatite (F<sup>out→i</sup><sub>depo</sub>) soil P pools</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>P deposition falling within labile and apatite soil P pools were derived from the Wang et al.<sup>49</sup> dataset that provided the atmospheric P deposition resulting from mineral dust, primary biogenic aerosol particles, sea salt, natural combustion and anthropogenic combustion, averaged over different time-periods. Year-to-year variability of P deposition from anthropogenic combustion was computed for information about emissions. Following Mahowald et al.<sup>50</sup>, we considered that 10% of P deposition from dust and 50% of deposition from other sources falls within P<sub>i-lab</sub> while rest falls within P<sub>i-prim</sub></td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td>None</td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>P deposition falling within labile and apatite soil P pools were derived from the combination of the dataset of Wang et al.<sup>49</sup> and the one provided by Wang et al.<sup>51</sup>. For year before 2007, the same approach as in GPASOIL-v0 was used while the dataset of Wang et al.<sup>51</sup> was used to compute the year-to-year variability of P deposition from anthropogenic combustion after the year 2007. The estimates of Wang et al.<sup>51</sup> were corrected to ensure equality with Wang et al.<sup>49</sup> for the years in common (1997-2007). See Eqs. 54, 55 and Text S4.</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>Estimates corresponding to Wang et al.<sup>49,51</sup> were provided by Rong Wang, pers.comm (2020). The monthly estimates of Wang et al.<sup>49</sup> were representative to different time-periods between 1960 and 2011, depending on the variables considered (dust deposition, etc.). Total deposition at monthly time-scale was provided for any years between 1997 and 2013 in Wang et al.<sup>51</sup>. The different estimates of deposition are spatially explicit without land-use (cropland, grassland, other) distinction. Emissions related to natural and anthropogenic combustion were provided at annual time-scale for any years between 1960 and 2007 for big World regions in Wang et al.<sup>49</sup>.</td></tr></tbody></table>
<table><caption>Table 3. SLUD driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>SLUD (P in sludges)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>P in sludges from sewage treatment that are spread on cropland soils</td></tr><tr><td>Corresponding variables in our approach</td><td>fP<sup>out-m</sup><sub>slud</sub> with m in {i-lab, o-lab, o-sta}</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>Not considered in Ringeval et al.<sup>11</sup></td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td></td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>fP<sup>out-x-tot</sup> was estimated by combining human P excretions with the fraction of sewage sludge that is treated and the P removal efficiency of treatments plants following van Puijnenbroek et al.<sup>52</sup> and Demay et al.<sup>9</sup>. Human excretions were computed by using excretion rates and human population. Computations were performed at country-scale then all crop within the same country receive the same P (in kgP ha<sup>-1</sup>). For the composition of P in sludges, we used the same labile vs stable contribution as for manure. See Eq. 56 and Text S5.</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>Human P excretion rates for big World regions, as well as fraction of sewage sludge that is treated, both at country scale, were provided by van Puijnenbroek et al.<sup>52</sup>. Human P excretion rates were available for the years 1970 and 2010 while the fraction of sewage sludge that is treated was available for years 1990, 2000, 2010. Human population was provided by HYDE 3.2<sup>25</sup> and Area at country scale was computed from the dataset used in LUCC (Table 5). The corresponding between big World regions and iso country code was provided by the IMAGE framework region classification (https://models.pbl.nl/image/index.php/Region_classification_map).</td></tr></tbody></table>

<table><caption>Table 4. LOSS driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>LOSS (P losses through soil erosion)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>Losses of P due to water erosion and runoff processes</td></tr><tr><td>Corresponding variables in our approach</td><td>Losses of P for each soil P pool (fP<sup>out</sup><sub>loss</sub> with m in {i-lab, o-lab, o-sta, i-prim, i-sec, x-occ}).</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>Losses were computed using the soil P content simulated in GPASOIL-v0, the sediment gross fluxes resulting from erosion provided by van Oost et al.<sup>74</sup> and the weight of top 0-0.3 m soil computed from bulk density of the same horizon provided by Soilgrids database (ISRIC - World Soil Information, 2016)</td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td>We did not consider temporal change in the sediment gross fluxes resulting from erosion (but we considered historical changes in agricultural land area, see LUCC)). Also, gross (and not net) erosion rates were considered. Finally, we considered that the erosion rates concerns exclusively the top soil layer (top 0-0.3 m).</td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>The same approach as used in (data = GPASOIL-v0) was used. Sediment gross fluxes resulting from erosion with a cropland/grassland distinction were computed following the approach described in Borrelli et al.<sup>54</sup> Losses through erosion of P in soil solution (P<sub>i-stol</sub>) was set to 0. We made the same assumptions as done in GPASOIL-v0: no temporal change in the gross rate of erosion, focus on gross (instead of net) erosion rate, effect of erosion exclusively on the top layer considered in our study. See Eqs. 57, 58.</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>Dataset were provided by Pasquale Borrelli, pers.comm, (2020) and corresponded to the use of the approach described in Borrelli et al.<sup>54</sup> here applied to the land-use dataset described in LUCC for year 2000 to simulate the gross soil losses (in (kg of soil) yr<sup>-1</sup> ha<sup>-1</sup>) through erosion with a cropland/grassland distinction at 250mx250m resolution. Bulk density and volumetric fraction of coarse fragment were provided by Soilgrids 2.0<sup>21</sup> at half-degree resolution.</td></tr></tbody></table>

<table><caption>Table 5. LUCC driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>LUCC (land use and land cover change)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>Decrease/increase of agricultural soil P resulting from land conversion within the same grid-cell (see Eq. 18).</td></tr><tr><td>Corresponding variables in our approach</td><td>Fractions of cropland (frac(y, crop, g)) and of grassland (frac(y, grass, g)) for any year y and grid-cell g as well as transitions between the 4 land-use categories considered in our study (e.g. Δ<sub>i</sub>(y, j, g) is the conversion from j to i with both i and j in {crop, grass, nonagri, urban}).</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>Fractions and transitions were provided by the Land Use Harmonization dataset<sup>75</sup></td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td></td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>Fractions and transitions were provided by Chini et al.<sup>20</sup>, the update of the Land Use Harmonization dataset (the so-called LUH2-GCB2019 dataset). We made a corresponding between our land-use categories and Chini et al.<sup>20</sup> categories. See Eqs. 59–62</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>Chini et al.<sup>20</sup> provided LUH2-GCB2019 (https://doi.org/10.3334/ORNLDAAC/185176), an update of the Land Use Harmonization dataset used in GPASOIL-v0. The LUH2-GCB2019 dataset provided fractions and transitions between crops (C3/C4, annual/perennial, C3-nitrogen fixing), grassland (managed pasture, rangeland), forest and urban land, at 0.25°x0.25° spatial resolution for any years between 850 and 2018.</td></tr></tbody></table>

Independent of any driver, three other soil properties were used in our approach to perform change in unit, namely the bulk density of the fine earth fraction, the volumetric fraction of coarse fragments (>2 mm) and the soil depth. These three variables were either used to convert concentration (i.e. per kg of soil) into quantity
<table><caption>**Table 6.** CLIM driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>CLIM (near-surface air temperature, soil temperature and soil water content)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>Few effects can be distinguished:<br/>- Effect of soil temperature and relative soil water content on P weathering and P mineralization in both (model = GPASOIL-v0) and (model = GPASOIL-v1) (See Eqs. 9–11)<br/>- Near-surface atmospheric temperature and soil water content involved in parameterizations of (model = GPASOIL-v1) (see Table 10)<br/>- Soil water content used to translate P<sub>t-soil</sub> into soil solution P concentration in (model = GPASOIL-v1) (see e.g. Equations 3 and 6)</td></tr><tr><td>Corresponding variables in our approach</td><td>Relative soil liquid water content (W<sub>rel</sub>, in fraction of field capacity) or soil liquid water content (W<sub>abs</sub>, in L(kg of soil)<sup>-1</sup>) soil temperature (T<sub>soil</sub>, in °C), T<sub>a</sub>, near-surface atmospheric temperature (in °C). W<sub>rel</sub>, W<sub>abs</sub>, T<sub>soil</sub> are representative to the 0–0.3 m soil layer.</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>The variables were computed by averaging corresponding variables simulated by two Dynamic Global Vegetation Models (ISBA<sup>27</sup> and ORCHIDEE<sup>24</sup>). The annual average of the climatology computed for the 1979–2010 period (i.e. no year-to-year variability) was used.</td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td>Only two Dynamic Global Vegetation Models were considered. No temporal change was considered (the 1970–2010 climatology was used instead).</td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>We used the average among 9 simulations (combination between 4 land-surface models and 3 climate data used as input of the land-surface models for the historical period) performed for the CMIP-6 exercise. These simulations have been considered as they provide the variables at the basis of our computation of W<sub>rel</sub>, W<sub>abs</sub>, and T<sub>soil</sub>. For each variable, we computed the annual averages over the 1850–2012 period. Variables used are representative to the top 0–0.3 m soil layer. Only liquid content was considered to compute W<sub>rel</sub> and W<sub>abs</sub>. T<sub>a</sub> is given by forcing files used as input of the land-surface models.</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>Original resolution and time-period considered varied among the 9 simulations considered. The link (with doi) of each simulation is provided below:<br/>CNRM-ESM2-1 x land-hist: https://doi.org/10.22033/ESGF/CMIP6.9599<br/>CNRM-ESM2-1 x land-hist-cruNcep: https://doi.org/10.22033/ESGF/CMIP6.9600<br/>CNRM-ESM2-1 x land-hist-princeton: https://doi.org/10.22033/ESGF/CMIP6.9601<br/>CNRM-CM6-1 x land-hist: https://doi.org/10.22033/ESGF/CMIP6.4095<br/>CNRM-CM6-1 x land-hist-cruNcep: https://doi.org/10.22033/ESGF/CMIP6.4100<br/>CNRM-CM6-1 x land-hist-princeton: https://doi.org/10.22033/ESGF/CMIP6.4101<br/>IPSL-CM6A-LR x land-hist: https://doi.org/10.22033/ESGF/CMIP6.5205<br/>MIROC6 x land-hist-cruNcep: https://doi.org/10.22033/ESGF/CMIP6.5627<br/>MIROC6 x land-hist-princeton: https://doi.org/10.22033/ESGF/CMIP6.5628</td></tr></tbody></table>

<table><caption>**Table 7.** SPRO driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>SPRO (soil properties)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>Soil properties involved in the computation of the soil P dynamics</td></tr><tr><td>Corresponding variables in our approach</td><td>Variables considered within the SPRO drivers vary between GPASOIL-v0 and GPASOIL-v1:<br/>- in (model = GPASOIL-v0), SPRO encompasses parameters involved in the Langmuir equation used to describe the equilibrium between P<sub>t-lab</sub> and P<sub>t-sec</sub> (so-called K<sub>s</sub> and S<sub>max</sub> parameters)<br/>- in (model = GPASOIL-v1), SPRO encompasses soil properties involved in parameterizations for soil P dynamics: soil texture (sand, clay, silt percentages: s<sub>a</sub>, s<sub>c</sub>, s<sub>i</sub> respectively, in %), soil water pH (pH, no unit), C is soil carbon concentration (in gC.(kg of soil)<sup>-1</sup>)</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>A coupled of values for (K<sub>s</sub>, S<sub>max</sub>) was used for each soil order. These values were provided by Wang et al.<sup>15</sup>. The global distribution of soil orders was similar to the one used in Yang et al.<sup>1</sup>.</td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td>Soil orders are likely a poor predictor of soil P buffering capacity in (model = GPASOIL-v0). Other properties (texture, pH, carbon) were not used in (data = GPASOIL-v0).</td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>K<sub>s</sub>, S<sub>max</sub> are not used any more in (model = GPASOIL-v1).<br/>Soil texture, soil water pH, and soil carbon concentration for top 0.3 m were get from Soilgrids 2.0<sup>21</sup>.<br/>We assumed that that soil properties at half-degree resolution could be applied to the cropland/grassland (soils) fraction. Soil properties were involved in equations described in Table 10.</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>Soil texture, soil water pH, and soil carbon concentration for top 0.3 m were computed by averaging values for 0–0.05, 0.05–0.15 and 0.15–0.30 m soil layers provided by Soilgrids 2.0<sup>21</sup>. The Soilgrids procedure allows to download the data at the resolution needed here (half degree spatial resolution).</td></tr></tbody></table>

per soil surface for both P or water, or to compute soil eroded in fraction of soil mass lost per year. In particular, soil P input/output were mainly given in the different dataset in kgP ha<sup>-1</sup> yr<sup>-1</sup> and soil P dynamics worked with P concentration in mgP (kg soil)<sup>-1</sup>. Given the scarcity of dataset about the thickness of the plough layer at the global scale, we considered a globally uniform thickness of 0.3 m as in Ringeval et al.<sup>11</sup>. This soil layer was also considered as encompassing a major proportion of crop roots. The soil bulk density at 0.5° latitude and longitude resolution from Soilgrids 2.0<sup>21</sup> was used to approach the density of cropland/grassland of each grid-cell given the unavailability of global datasets focusing on agricultural soils, even though it is known that soil treatment has an effect on soil physical properties<sup>22</sup>. Same reasoning applied to the fraction of coarse fragments, also provided by
<table><caption>**Table 8.** BIOG driver and its representation in (data = GPASOIL-v0) and (data = GPASOIL-v1).</caption><thead><tr><th>Driver name</th><th>BIOG (natural soil biogeochemical background)</th></tr></thead><tbody><tr><td>Description of its effect on the soil P</td><td>P inherited from natural soils at the time of conversion to agriculture (Eq. 20). P in natural soils was also used to approach P in agricultural soils at the beginning of the simulation (initial conditions, Eqs. 22, 23) and soil P pools at steady-state used to compute the parameter describing the exchanges between pools (Eq. 17).</td></tr><tr><td>Corresponding variables in our approach</td><td>P content of natural soils for any soil P pools considered, i.e. P<sup>NA</sup><sub>i-sol</sub>, P<sup>NA</sup><sub>i-lab</sub>, P<sup>NA</sup><sub>o-lab</sub>, P<sup>NA</sup><sub>o-sta</sub>, P<sup>NA</sup><sub>i-prim</sub>, P<sup>NA</sup><sub>i-sec</sub>, P<sup>NA</sup><sub>i-occ</sub>, in kgP ha<sup>-1</sup> for top 0-0.3 m.</td></tr><tr><td>Dataset ref used in (data = GPASOIL-v0)</td><td>Natural soils P pools for top 0-0.3 m were approached by estimates of current P in unmanaged soils for top 0-0.5 m provided by Yang et al.<sup>4</sup>. P<sub>i-sol</sub> was not considered in GPASOIL-v0.</td></tr><tr><td>Issues related to the use of this dataset in GPASOIL-v0</td><td>Soil orders used to compute how total P is held in different fractions in Yang et al.<sup>4</sup> is likely a poor predictor of soil P pools. In our approach, we assumed that P concentration provided for top 0-0.5 m was representative to the considered top 0-0.3 m soil layer.</td></tr><tr><td>Representation in (data = GPASOIL-v1)</td><td>We used the dataset of He et al.<sup>14</sup> that provides the current soil P distribution for P<sup>NA</sup><sub>i-sol</sub>, P<sup>NA</sup><sub>i-sec</sub>, P<sup>NA</sup><sub>i-prim</sub>, P<sup>NA</sup><sub>i-occ</sub>, P<sup>NA</sup><sub>i-occ</sub>, P<sup>NA</sup><sub>o-sta</sub>, P<sup>NA</sup><sub>o-sta</sub>. The values we used are representative to the top 0 - 0.3 m soil layer. P<sub>i-sol</sub> is derived from P<sub>i-∞</sub> that was prescribed to 0.1 mgP L<sup>-1</sup> (but sensitivity to this value was tested).</td></tr><tr><td>Characteristics of the spatially explicit dataset used in (data = GPASOIL-v1)</td><td>He et al.<sup>14</sup> provided the current soil P pools distribution (in mgP (kg of soil)<sup>-1</sup>) for different soil horizons at half-degree spatial resolution for natural ecosystems (78, https://doi.org/10.6084/m9.fishshare.16988029.v2).</td></tr></tbody></table>

<table><caption>**Table 9.** Corresponding between pools name and Hedley fractions in (model = GPASOIL-v0) and (model = GPASOIL-v1). Squared brackets in columns 3 and 6 are used for Hedley fractions which are not systematically quantified. Columns 2 and 3 were adapted from the Table 2 of Yang and Post<sup>79</sup>. Please, note that all occluded P was considered as inorganic P in Ringeval *et al.*<sup>11</sup> while this was not stated in Yang and Post<sup>79</sup>.</caption><thead><tr><th>Pools name used in (model = GPASOIL-v0)</th><th>Pools definition used in (model = GPASOIL-v0)</th><th>Hedley fractions considered on sites used in Ringeval et al.<sup>11</sup></th><th>Pools name used in (model = GPASOIL-v1)</th><th>Pools definition used in (model = GPASOIL-v1) (following He et al.<sup>14</sup>)</th><th>Corresponding with Hedley fraction used in He et al.<sup>14</sup></th></tr></thead><tbody><tr><td rowspan="2">P<sub>ILAB</sub></td><td rowspan="2">Labile inorganic P</td><td rowspan="2">[H<sub>2</sub>O Pi/Resin Pi+] Bicarbonate Pi</td><td>P<sub>i-sol</sub></td><td>Inorganic P in soil solution</td><td rowspan="2">H<sub>2</sub>O Pi/Resin Pi + Bicarbonate Pi</td></tr><tr><td>P<sub>i-lab</sub></td><td>Labile inorganic P</td></tr><tr><td>P<sub>SEC</sub></td><td>Inorganic P bound on secondary minerals</td><td>Hydroxide Pi [+Sonic Pi]</td><td>P<sub>i-sec</sub></td><td>Moderately inorganic P</td><td>Hydroxide Pi</td></tr><tr><td>P<sub>OLAB</sub></td><td>Labile organic P</td><td>[H<sub>2</sub>O Po + Resin Po+] Bicarbonate Po</td><td>P<sub>o-lab</sub></td><td>Labile organic P</td><td>Bicarbonate Po</td></tr><tr><td>P<sub>OSTA</sub></td><td>Stable organic P</td><td>Hydroxide Po [+Sonic Po + HCl Po]</td><td>P<sub>o-sta</sub></td><td>Moderately organic P</td><td>Hydroxide Po</td></tr><tr><td>P<sub>APA</sub></td><td>Apatite</td><td>HCl Pi (HCl diluted or not)</td><td>P<sub>i-prim</sub></td><td>Primary inorganic P</td><td>HCl Pi (HCl diluted only)</td></tr><tr><td>P<sub>Occ</sub></td><td>Occluded inorganic P</td><td>Residual P [+Hot HCl]</td><td>P<sub>x-occ</sub></td><td>Occluded P (inorganic + organic)</td><td>Residual P [+Hot HCl P + Sonic P]</td></tr></tbody></table>

Soilgrids 2.0. The uncertainty associated to the bulk density and fraction of coarse fragments was not considered in our study as we restricted the uncertainty analysis to the one of the drivers. In the present study, all P pools were finally expressed in kgP ha<sup>-1</sup> for top 0-0.3 m soil layer. Large-scale soil P pools are computed by using cropland and grassland fraction (see LUCC driver below) and grid-cell land area (computed from regriding the land fraction provided by HYDE 3.2<sup>23</sup>, https://doi.org/10.17026/dans-25g-gez3).

The new estimates of the soil P distribution in cropland and grasslands are called GPASOIL-v1 and we used the name GPASOIL-v0 to describe the estimates calculated in the previous study of Ringeval *et al.*<sup>11</sup>. Estimates result from a coupling between dataset describing drivers (called “data” in the following) and a soil P dynamics model (called “model” hereafter), and v0 or v1 can be attributed to each component (data or model) as each one has been updated from the study of Ringeval *et al.*<sup>11</sup> to the current one.

**Soil P dynamics model.** *Pools and fluxes design.* The soil P pools considered in (model = GPASOIL-v1) were named with the following nomenclature: P<sub>a-b</sub> with a in {i,o,x} (i for inorganic, o for organic and x for inorganic + organic (i + o)); and b referring to different types of soil P pools. In total, 7 pools were considered following the merging of Hedley fractions (Table 9, Fig. 1b): P<sub>i-sol</sub> (inorganic P in soil solution), P<sub>i-lab</sub> (labile inorganic P), P<sub>i-sec</sub> (moderately labile inorganic P), P<sub>i-prim</sub> (primary inorganic P), P<sub>o-lab</sub> (labile organic P), P<sub>o-sta</sub> (stable organic P), P<sub>x-occ</sub> (occluded P). P<sub>x-tot</sub> is the sum of all pools considered.

Hereafter, the following nomenclature was used to name the flux mediated by the process “proc” from pool m1 to m2: $fP_{proc}^{m1→m2}$. The superscripts m1 and m2 correspond to soil P pools in case of soil P dynamics or could be used to describe soil P input/output (« out » for outside is used in that case). The name “proc” was used to describe the process that lead to the flux considered and either referred to a soil process (“occl” for occlusion, “deoccl” for deocclusion, “sorp” for sorption, “desorp” for desorption, “weat” for weathering, “mine” for mineralization), or to a soil P input/output (“chem” for chemical fertilizer, “manu” for manure, “resi” for plant residues, “upta” for plant uptake, “loss” for losses through water erosion, “depo” for atmospheric deposition, “slud” for sludges from waste treatment). As examples, $fP_{mine}^{o-sta→i-lab}$ corresponded to the mineralization of P<sub>o-sta</sub> to P<sub>i-lab</sub> while $fP_{chem}^{out→i-lab}$ corresponded to soil P input of chemical fertilizer to P<sub>i-lab</sub>, $fP_{lucc}$ was used to name the effect of land-use change on soil P (section “Effect of land-use and land-cover change (LUCC) on cropland/grassland soil P”).
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Difference in design between soil P pool dynamics model used in this study (model = GPASOIL-v1, panel (b)) and used in Ringeval *et al.*¹¹ (model = GPASOIL-v0, panel (a)). Inorganic pools are in blue, organic ones are in orange and grey pools correspond to pools which encompass both inorganic and organic P forms. Double arrows means that an equilibrium is considered. Nomenclature used to name the pools changed between (model = GPASOIL-v0) and (model = GPASOIL-v1).</figcaption></figure>

<table><caption>**Table 10.** Original parametrizations provided by Wang *et al.*¹⁸ and modified parametrizations used in this study after i) exclusion of oxalate concentrations from explanatory variables and ii) setting *k* parameters as constant in time. All parameterizations are built from a stepwise multivariate linear regression against the database used in Wang *et al.*¹⁸. r² provided in the last column characterize either the original parametrizations (and thus they are equal to values given in Table 2 of Wang *et al.*¹⁸) or parametrizations after excluding oxalate from explanatory variables (this study). The parameters *k* were kept constant in time through i) the use of a constant-in-time value for *f*<sub>m</sub> for any soil pool *m* and, ii) the use of *P*<sub>*o*-tot,<sub>∞</sub> and *P*<sub>*i*-tot\prim,<sub>∞</sub> in the parametrizations instead of *P*<sub>*o*-tot and *P*<sub>*i*-tot\prim</sub>, respectively. Note: C is soil carbon concentration in gC (kg of soil)<sup>-1</sup>, O<sub>x</sub> is oxalate metal oxide (Al and Fe) concentration in mmol(kg of soil)<sup>-1</sup>, C<sub>x</sub> is the ratio of total soil C and oxalate metal oxide concentrations in g/mmol (C<sub>x</sub> = 0.001 C/O<sub>x</sub>), s<sub>s</sub>, s<sub>c</sub> and s<sub>i</sub> are sand, clay and silt percentages, respectively, T<sub>a</sub> is the mean annual near-surface air temperature in °C, pH is soil pH measured in water (no unit), P<sub>*o*-tot is the total soil organic P in mgP (kg of soil)<sup>-1</sup> and *P*<sub>*i*-tot\prim is the total inorganic P minus primary P in mgP (kg of soil)<sup>-1</sup>, *f*<sub>m</sub> denotes the ratio (no unit) at steady-state of pool *m* and *P*<sub>*i*-tot\prim,<sub>∞</sub> Subscript ∞ denotes pools at the steady-state.</caption><thead><tr><th>Parameters</th><th>Name in Wang et al.<sup>18</sup></th><th>unit</th><th>Equation</th><th>r<sup>2</sup></th></tr></thead><tbody><tr><td rowspan="2">P<sub>c,∞</sub></td><td rowspan="2">C<sub>∞</sub></td><td rowspan="2">mgP L<sup>-1</sup></td><td><sup>18</sup> P<sub>c,∞</sub> = 10<sup>[1.7425+1.7782(1-exp(-1.6653C<sub>x</sub>))]</sup></td><td rowspan="2">0.54</td></tr><tr><td>This study P<sub>c,∞</sub> = 0.1</td></tr><tr><td rowspan="2">b</td><td rowspan="2">b</td><td rowspan="2">dimensionless</td><td><sup>18</sup> b = exp(0.532 - 0.13O<sub>x</sub> - 37.49f<sub>i-sol</sub> + 0.81f<sub>i-lab</sub> + 0.003s<sub>s</sub>)</td><td>0.43</td></tr><tr><td>This study b = exp(-0.628 - 36.702f<sub>i-sol</sub> + 1.102f<sub>i-lab</sub> + 0.0024s<sub>s</sub>)</td><td>0.23</td></tr><tr><td rowspan="2">k<sup>i-lab→i-sol</sup></td><td rowspan="2">k<sub>WL</sub></td><td rowspan="2">day<sup>-1</sup></td><td><sup>18</sup> k<sup>i-lab→i-sol</sup> = -4.82 + 209f<sub>i-sol</sub> + 14.64f<sub>x-occ</sub> + 9.26f<sub>i-sec</sub> - 0.008C - 0.0003P<sub>i-tot</sub> + 0.018s<sub>i</sub></td><td>0.28</td></tr><tr><td>This study Same as<sup>18</sup> but use of P<sub>i-tot</sub> in instead of P<sub>i-tot</sub></td><td></td></tr><tr><td rowspan="2">k<sup>i-sol→i-sec</sup></td><td rowspan="2">k<sub>SW</sub></td><td rowspan="2">mgP (kg soil)<sup>-1</sup> day<sup>-1</sup> (mg P/L)<sup>-b</sup></td><td><sup>18</sup> k<sup>i-sol→i-sec</sup> = exp(0.002 + 4.0f<sub>i-sec</sub> + 0.0008P<sub>o-tot</sub> + 0.012C + 0.108T<sub>a</sub> - 0.0002P<sub>i-tot</sub>)</td><td>0.64</td></tr><tr><td>This study Same as<sup>18</sup> but use of P<sub>o-tot</sub> and P<sub>i-tot</sub> instead of P<sub>o-tot</sub> and P<sub>i-tot</sub>, respectively</td><td></td></tr><tr><td rowspan="2">k<sup>i-sec→x-occ</sup></td><td rowspan="2">k<sub>OS</sub></td><td rowspan="2">day<sup>-1</sup></td><td><sup>18</sup> k<sup>i-sec→x-occ</sup> = 1.6e<sup>-5</sup> + 0.0001f<sub>x-occ</sub> - 0.0001f<sub>i-sec</sub> + 1.3e<sup>-5</sup>pH - 9.0e<sup>-6</sup>Ox - 4.48e<sup>-7</sup>s<sub>c</sub> + 1.1e<sup>-7</sup>C</td><td>0.90</td></tr><tr><td>This study k<sup>i-sec→x-occ</sup> = 3.68e<sup>-5</sup> + 9.60e<sup>-5</sup>f<sub>x-occ</sub> - 1.47e<sup>-4</sup>f<sub>i-sec</sub> + 1.22e<sup>-5</sup>pH - 4.07e<sup>-7</sup>s<sub>c</sub> + 5.16e<sup>-8</sup>C</td><td>0.90</td></tr></tbody></table>

In (model = GPASOIL-v1), the equations describing soil inorganic P dynamics between P<sub>i-sec</sub>, P<sub>i-lab</sub>, P<sub>x-occ</sub>, P<sub>i-sol</sub> were given by Wang *et al.*¹⁸. Weathering and mineralization were simulated as in (model = GPASOIL-v0) but other mineralization rates values were here tested (see section “Parameter estimates”). In the present section, fP was expressed in mgP (kg soil)<sup>-1</sup> day<sup>-1</sup>, soil P pools (P<sub>i-sol</sub>, P<sub>i-sec</sub>, etc.) in mgP (kg soil)<sup>-1</sup> and soil water content (W<sub>abs</sub>) in L (kg of soil)<sup>-1</sup> to make parameters unit consistent to the units used in Wang *et al.*¹⁸ (Table 10). The different unit conversions were done by using the bulk density (in kg m<sup>-3</sup>), volumetric fraction of coarse fragments (no unit) and spatially constant soil depth (0.3 m). Fluxes and pools varied as functions of year (y), land-cover (lu, either equal to crop or grass) and space (g, for grid-cell) (i.e. fP = fP(y,lu,g)). Parameters involved in Eqs. 1–8 (e.g. k<sup>i-sec→x-occ</sup> and b) varied as a function of grid-cell only and W<sub>abs</sub> as a function of year and grid-cell; but for the sake of readability, we omit these in the below equations.
As in Wang et al.¹⁸, occlusion/deocclusion are given by Eqs. 1, 2 and sorption/desorption between Pᵢ-sol and Pᵢ-sol are computed with Eqs. 3, 4:

$$fP_{occl}^{i-sec \to x-occ} = k^{i-sec \to x-occ} \cdot P_{i-sec} \quad (1)$$

$$fP_{deoccl}^{x-occ \to i-sec} = k^{x-occ \to i-sec} \cdot P_{x-occ} \quad (2)$$

$$fP_{sorp}^{i-sol \to i-sec} = k^{i-sol \to i-sec} \cdot (P_{i-sol} / W_{abs})^b \quad (3)$$

$$fP_{desorp}^{i-sec \to i-sol} = k^{i-sec \to i-sol} \cdot P_{i-sec} \quad (4)$$

where Pᵢ-sol/Wₐbs corresponds to the P concentration in the soil solution (in mgP L⁻¹), k (in day⁻¹) and b (unitless) are parameters. Sorption follows a Freundlich equation.

In Wang et al.¹⁸, sorption/desorption fluxes between Pᵢ-sol and Pᵢ-lab are computed as follows:

$$fP_{desorp}^{i-lab \to i-sol} = k^{i-lab \to i-sol} \cdot P_{i-lab} \quad (5)$$

$$fP_{sorp}^{i-sol \to i-lab} = k^{i-sol \to i-lab} (P_{i-sol} / W)^b \quad (6)$$

To simplify the application of the¹⁸ model at the global scale in the current study, an equilibrium between Pᵢ-sol and Pᵢ-lab is assumed at daily time-step, following Wang et al.¹⁵. Thus, kᵢ-lab→i-solPᵢ-lab = kᵢ-sol→i-lab(Pᵢ-sol/Wₐbs)ᵇ i.e.:

$$P_{i-lab} = \frac{k^{i-sol \to i-lab}}{k^{i-lab \to i-sol}} (P_{i-sol} / W_{abs})^b \quad (7)$$

Numerical resolution of Eq. 7 was done by substituting Pᵢ-sol,lab = Pᵢ-sol + Pᵢ-lab in Eq. 7 and rearranging to get:

$$k^{i-sol \to i-lab} (P_{i-sol})^b + k^{i-lab \to i-sol} (W_{abs})^b P_{i-sol} - k^{i-lab \to i-sol} (W_{abs})^b P_{i-sol,lab} = 0 \quad (8)$$

where the unique unknown is Pᵢ-sol. Equation 8 was used to redistribute Pᵢ-sol,lab between Pᵢ-sol and Pᵢ-lab by assuming a steady-state between Pᵢ-sol and Pᵢ-lab. Equation 8 was solved using scipy.optimization root solver in python3.6.

The Wang et al.¹⁸ model excludes inorganic primary P and organic pools, and their exchanges with other inorganic pools. To model weathering and mineralization in (model = GPASOIL-v1), the same equations as in (model = GPASOIL-v0) were used. The weathering is computed as follows:

$$fP_{weat}^{i-prim \to i-lab} = k^{i-prim \to i-lab} * g_1(T_{soil}) * g_2(W_{rel}) * P_{i-prim} \quad (9)$$

where g₁ and g₂ (unitless) described the sensitivity to soil temperature (Tₐoil, in °C) and relative soil water content (Wₐrel, unitless), respectively. As detailed in Ringeval et al.¹¹, g₁(T) is described with a Q₁₀ of 2.4 and a reference temperature of 15 °C and g₂(Wₐrel) = Wₐrel. The relative soil water content is defined by:

Wₐrel = Wₐbs/Wₐfc with Wₐfc representing the field water capacity (in L (kg of soil)⁻¹). Dataset used to get values of Wₐbs, Wₐfc, Tₐoil are described in the CLIM driver.

Fluxes from organic to inorganic pools (mainly assimilated to mineralization) are given by:

$$fP_{mine}^{o-sta \to i-lab} = k^{o-sta \to i-lab} * h_1(T_{soil}) * h_2(W_{rel}) * P_{o-sta} \quad (10)$$

$$fP_{mine}^{o-lab \to i-lab} = k^{o-lab \to i-lab} * h_1(T_{soil}) * h_2(W_{rel}) * P_{o-lab} \quad (11)$$

where h₁ and h₂ represent the sensitivity to soil temperature and soil water content and kᵢ-sta→i-lab and kᵢ-sta→i-lab are turnover rates for stable and labile pools (in day⁻¹), respectively. The same functions as the ones used in the ORCHIDEE²⁴ dynamic global vegetation model to characterize the carbon mineralization were prescribed to h₁(Tₐoil) and h₂(Wₐrel) (i.e. Q₁₀ = 2 and reference temperature = 30 °C for h₁; h₂(Wₐrel) = -1.1(Wₐrel)² + 2.4Wₐrel - 0.29).

At the 1st day of each year, soil P input/output (chemical fertilizer, manure, plant residues, plant uptake, atmospheric deposition, sludges, losses by water erosion) were added/removed to/from the different soil P pools. Each soil input/output had a given composition that defines which pools were enriched or depleted, as provided later. The only exception to this rule was Pᵢ-sol which was neither enriched nor depleted by soil P input/output in the 1st day. Any P fluxes that should have reached or leave Pᵢ-sol was set to occur with the Pᵢ-lab pool instead. For instance, while it is well known that plant only takes up P in soil solution, our modeling approach did not allow to simulate daily plant growth and related P uptake from P in solution. Instead, our methodological choice was to remove annual plant uptake from Pᵢ-lab at the 1st day of the year, then to make Pᵢ-lab and other soil pools interact for each following day of the rest of the year. Within each time-step (day), the following order of processes was considered: exchanges between pools/forms were first computed, then pools were updated according to these changes and finally, the equilibrium between Pᵢ-sol and Pᵢ-lab was computed.
To avoid over depletion of a given soil P pool, i.e. the net output flux at the model time step exceeding the pool size, we restricted this flux to the pool size. The effect of this restriction will be evaluated in the section « Technical validation ». After removing $fP_{upta}^{i-lab→out}$ from P$_{i-lab}$ and before computing the equilibrium between P$_{i-sol}$ and P$_{i-lab}$, if (P$_{i-sol}$ + P$_{i-lab}$) was negative, (P$_{i-sol}$ + P$_{i-lab}$) was set to 0 and the uptake not satisfied by the soil P pools simulated was saved to a variable called $fP_{upns}$. This variable corresponded to the part of the uptake which was prescribed to the soil P dynamics model but that the model was not able to sustain. In that case, $fP_{upta}^{i-lab→out}$ was corrected to the maximum uptake that was determined by the size of (P$_{i-sol}$ + P$_{i-lab}$).

*Parameter estimates.* Equations 1–8 are based on Wang *et al.*¹⁸. To estimate the flux parameters b, $k^{i-sol→i-sol}$, $k^{i-sol→i-sec}$, and $k^{i-sec→x-occ}$, Wang *et al.*¹⁸ proceed in two steps by first calibrating these parameters against isotopic exchange kinetics experiments from 147 sites. Then, they regressed such optimized parameters against pedo-climatic variables. They found that b, $k^{i-lab→i-sol}$, $k^{i-sol→i-sec}$, $k^{i-sec→x-occ}$ correlate with soil P fractions of different inorganic pools (the term “fraction” is used here to define the ratio between a given pool and the sum of all pools considered in Wang *et al.*¹⁸, see the variables called f in Eqs. 12–14 and definition in Eq. 15), soil organic carbon and oxalate extractable metal oxide concentrations. Parametrizations found in Wang *et al.*¹⁸ are given in Table 10. Here, in order to estimate these parameters at the global scale at which oxalate extractable metal oxide concentrations is not available, we re-compute the regressions after excluding oxide concentrations from explanatory variables (Table 10).

Wang *et al.*¹⁸ derived other parameters ($k^{i-sol→i-lab}$, $k^{i-sec→i-sol}$, $k^{x-occ→i-sec}$) from relationships with b, $k^{i-lab→i-sol}$, $k^{i-sol→i-sec}$, and $k^{i-sec→x-occ}$. These relationships were also used in our study:

$$ \frac{k^{i-sol→i-lab}}{k^{i-lab→i-sol}} = \frac{f_{i-lab} \overline{W}_{abs}}{b(P_{c,\infty})^{b-1} f_{i-sol}} \qquad (12) $$

$$ \frac{k^{i-sol→i-sec}}{k^{i-sec→i-sol}} = \frac{f_{i-sec} \overline{W}_{abs}}{b(P_{c,\infty})^{b-1} f_{i-sol}} \qquad (13) $$

$$ \frac{k^{x-occ→i-sec}}{k^{i-sec→x-occ}} = \frac{f_{i-sec}}{f_{x-occ}} \qquad (14) $$

where $P_{c,\infty}$ is the concentration of P in soil solution at steady-state, overlined $W_{abs}$ means temporal average of $W_{abs}$, $f_m$ denotes the ratio (no unit) at steady-state of pool m and the sum of inorganic pools considered in Wang *et al.*¹⁸ called here $P_{i-tot\text{prim}}$, e.g. for m = i-sec,

$$ f_{i-sec} = \frac{P_{i-sec;\infty}}{P_{i-tot\text{prim};\infty}} = \frac{P_{i-sec;\infty}}{P_{i-sec;\infty} + P_{i-sol;\infty} + P_{x-occ;\infty} + P_{i-lab;\infty}} \qquad (15) $$

with subscript ∞ denotes pools at the steady-state. Note that Wang *et al.*¹⁸ made distinctions between measured (f*) and modeled (f) fractions but this distinction does not make sense in this study where fractions are only simulated. Wang *et al.*¹⁸ also considered $W_{abs}$ equal to 10 L (kg of soil)⁻¹ representative to water content of isotopic dilution while we used here the averaged over the whole simulation (1900–2018) of $W_{abs}$. The parameter $P_{c,\infty}$ that describes the concentration of P in soil solution at steady-state depends only on oxalate in Wang *et al.*¹⁸ and we prescribed here an arbitrary constant value:

$$ P_{c,\infty} = 0.1 \text{ mgPL}^{-1} \qquad (16) $$

which is equal to the median value of Helfenstein *et al.*²⁵. The sensitivity to the value chosen for $P_{c,\infty}$ was assessed in the “Technical validation” section. Dataset for calculating soil parameters using the relationships in Table 10 are given in the SPRO driver description (Table 7).

The rate parameters (different k) and parameter b are only grid-cell dependent and do not differ between cropland and grassland. Following parameterizations found in Wang *et al.*¹⁸, these parameters could theoretically vary in time with temporal variations of the different soil P fractions at equilibrium (f), organic P ($P_{o-tot}$) and the sum of total inorganic pools considered in Wang *et al.*¹⁸ (i.e. $P_{i-tot\text{prim}}$) (see Table 10). We choose here to keep them constant through i) the use of time-invariant value for $P_{m,\infty}$ (and thus for $f_m$) for any soil pool m and, ii) the use of $P_{o-tot,\infty}$ and $P_{i-tot\text{prim},\infty}$ in Table 10 instead of $P_{o-tot}$ and $P_{i-tot\text{prim}}$, respectively. $P_{m,\infty}$ is assumed to have the current value of unmanaged soil P pools $P_m^{NA}$ given by He *et al.*¹⁴ for the same grid-cell g, as follows:

$$ P_{m,\infty}(g) = P_m^{NA}(g) \qquad (17) $$

in a similar way to the computation of initial conditions (see the section about LUCC).

To summarize, the main differences between the original Wang *et al.*¹⁸ model and the model used here for inorganic soil P dynamics are the following:
<table><thead><tr><td>Name of the soil P estimates</td><td>data</td><td>model</td></tr></thead><tbody><tr><td>GPASOIL-v0<sup>11</sup></td><td>GPASOIL-v0</td><td>GPASOIL-v0</td></tr><tr><td>Not named in the manuscript</td><td>GPASOIL-v1</td><td>GPASOIL-v0</td></tr><tr><td>GPASOIL-v1.0</td><td>GPASOIL-v1</td><td>GPASOIL-v1</td></tr><tr><td>GPASOIL-v1.1</td><td>GPASOIL-v1</td><td>GPASOIL-v1.1 corresponds to GPASOIL-v1 with increased mineralization rates</td></tr></tbody></table>

**Table 11.** Name of the soil P estimates used in this study. The soil P estimates results from a coupling between dataset describing drivers (“data”) and a soil P dynamics model (“model”), and different combinations (data x model) were tested in the manuscript. “v0” refers to Ringeval *et al.*¹¹ while “v1” refers to the update presented in the current draft.

- we assumed an equilibrium between $P_{i-sol}$ and $P_{i-lab}$ (Eq. 8) instead of (Eqs. 5-6)

- our parametrization excluded oxalate as explanatory variables (thus simplified equations for $b, k^{i-lab→i-sol}, k^{i-sol→i-sec}$, and $k^{i-sec→x-occ}$, as well as a constant in space $P_{c,∞}$; Table 10) and we considered that all k parameters as constant in time.

The weathering rate ($k^{i-prim→i-lab}$) was set to 2.7e-7 day⁻¹ following Buendia *et al.*²⁶ as in (model = GPASOIL-v0). Regarding mineralization rates, $k^{o-sta→i-lab} = 2.7e^{-5} day^{-1}$ and $k^{o-lab→i-lab} = 2.7e^{-4} day^{-1}$ were prescribed for organic stable and organic labile pools, i.e. a residence time of 100 and 10 years, respectively, as in (model = GPASOIL-v0). Given the uncertainty related to such rates²⁷, another values were tested: $k^{o-sta→i-lab} = 1.8e^{-4} day^{-1}$ and $k^{o-lab→i-lab} = 1.4e^{-3} day^{-1}$; i.e. a residence time of respectively 15 and 2 yr. Such smaller residence time of those two parameters allow the simulated soil P pool for grassland to become stable during the 1st half of the 20th century (see the “Technical validation” section). (model = GPASOIL-v1.0) was used to name the model with $k^{o-sta→i-lab} = 2.7e^{-5} day^{-1}$ and $k^{o-lab→i-lab} = 2.7e^{-4} day^{-1}$ while (model = GPASOIL-v1.1) was used for the model with $k^{o-sta→i-lab} = 1.8e^{-4} day^{-1}$ and $k^{o-lab→i-lab} = 1.4e^{-3} day^{-1}$ (Table 11).

**Summary of the main differences with (model = GPASOIL-v0).** Soil P pools considered in (model = GPASOIL-v0) were listed in Table 9 and both pools and fluxes considered were plotted in Fig. 1a. Main differences with (model = GPASOIL-v1) were represented in Fig. 1 and listed hereafter:

(i). (model = GPASOIL-v1) considered an additional inorganic P in soil solution ($P_{i-sol}$) that exchanges with both $P_{i-lab}$ and $P_{i-sec}$. $P_{i-sol}$ is an intermediary pool between $P_{i-sec}$ and $P_{i-lab}$ in Wang *et al.*¹⁸, and thus is key to represent processes of sorption/desorption. $P_{i-sol}$ would also allow a better coupling with P uptake model in further studies.

(ii). a Freundlich equation was used to describe the fluxes leaving $P_{i-sol}$ to either $P_{i-lab}$ or $P_{i-sec}$ (vs. a Langmuir equilibrium between $P_{i-sec}$ and $P_{i-lab}$ in (model = GPASOIL-v0)). Freundlich equations were assumed to be better (see the introduction of Wang *et al.*¹⁸).

(iii). in (model = GPASOIL-v0), all parameters describing soil P fluxes – except those involve in the Langmuir equilibrium, which vary as function of soil orders- were constant in space (Text S1). In (model = GPASOIL-v1), parameters involved in the soil inorganic P dynamics varied in space following Wang *et al.*¹⁸.

(iv). the flux from $P_{x-occ}$ to $P_{i-sec}$ was omitted in (model = GPASOIL-v0) but is resolved in (model = GPASOIL-v1)

(v). (model = GPASOIL-v0) was based on yearly time-step, while (model = GPASOIL-v1) was based on daily one, and

(vi). while soil organic P dynamics was similar between the two model versions, one alternative couple of residence time were tested in (model = GPASOIL-v1) (same residence time as in model = GPASOIL-v0 for GPASOIL-v1.0 while another one for GPASOIL-v1.1, Table 11).

**Effect of land-use and land-cover change (LUCC) on cropland/grassland soil P.** Within a given grid-cell, the change in the different land-cover fractions from one year to the other modified the soil P of cropland (or grassland), if there were some gross conversions from other land-cover to cropland (grassland respectively). Thus, the effect of LUCC on cropland soil P pools for year y and grid-cell g (called $P_m(y, crop, g)$ with subscript m referring to a given soil P pool) was simulated using Eq. 18 as follows:

$$ \frac{\text{frac}(y, \text{crop}, g) * P_m(y, \text{crop}, g)}{\text{frac}(y-1, \text{crop}, g) - \sum_{i \in lu\{\text{crop}\}} \Delta_i(y-1, \text{crop}, g)} * P_m(y-1, \text{crop}, g) + \sum_{j \in lu\{\text{crop}\}} (\Delta_{\text{crop}}(y-1, j, g) * P_m(y-1, j, g)) \quad (18) $$

with $frac(y, crop,g)$ was the fraction of the grid-cell g covered by cropland, $\Delta_i(y,j,g)$ was the change in cover area from land-use j to land-use i, $lu$ is the list of land-use considered (i.e. $lu = \{\text{crop, grass, nonagri, urban}\}$), $lu\{n\}$
means the list of land-use *lu* after removing {n}. A similar equation where replacing *crop* by *grass* described the effect of LUCC on grassland soil P pools. The increase/decrease in a given soil pool Pm of a given land-cover (cropland or grassland) through LUCC was called fP<sub>lucc</sub> and was positive when Pm increases from one year to the other. The definition of frac was the following:

$$frac(y, lu, g) = Area(y, lu, g)/land(g) \qquad (19)$$

with Area for the area of the grid-cell covered by *lu* and land the total land area in the grid-cell *g*. The soil P content of nonagri land-cover fraction was not explicitly simulated and was defined as follows:

$$P_m(y, nonagri, g) = P_m^{NA}(g) \qquad (20)$$

where P<sub>m</sub><sup>NA</sup> (g) is the current soil P pool *m* of natural (unmanaged) soils and grid-cell *g*. P<sub>m</sub><sup>NA</sup> (g) was provided by He *et al.*<sup>14</sup> for *m* within {P<sub>i-lab</sub>, P<sub>i-sec</sub>, P<sub>i-prim</sub>, P<sub>o-lab</sub>, P<sub>o-sta</sub>, P<sub>x-occ</sub>} . P<sub>i-sol</sub><sup>NA</sup> (g) was not provided by He *et al.*<sup>14</sup> and for this pool, we used the following equation:

$$P_{i-sol}^{NA}(g) = P_{c,\infty} \cdot \overline{W}_{abs}(g) \qquad (21)$$

with P<sub>c,∞</sub> defined in (Eq. 16) and W<sub>abs</sub> the soil water content, as defined earlier.

Through Eq. 20, we assumed that the P inherited from natural soils at the conversion to agriculture (the so-called biogeochemical background, BIOG) could be represented by prescribing the current P in unmanaged soils to all soils converted to agricultural soils over the last ~120 years. As explained in Ringeval *et al.*<sup>11</sup>, this could be limiting, in particular in regions where shifting cultivation occurred, leading to modification of the P content of soils covered by natural vegetation. We also neglected soil P input corresponding to forest biomass left on soil at the time of conversion<sup>28</sup>. All these simplifications had to be taken unless we explicitly simulated the change of soil P in non-agricultural areas, which would require the use of global vegetation models that rely on their own assumptions (e.g. Sun *et al.*<sup>29</sup>).

In addition, P in natural soils was also used as initial conditions, i.e.:

$$P_m(y_0, \text{crop}, g) = P_m^{NA}(g) \qquad (22)$$

$$P_m(y_0, \text{grass}, g) = P_m^{NA}(g) \qquad (23)$$

with y<sub>0</sub> corresponds to the 1<sup>st</sup> year of the simulation, chosen here equal to 1900.

The soil P content of urban land-cover fraction was not explicitly simulated and was defined by:

$$P_m(y, \text{urban}, g) = 0 \qquad (24)$$

Assumption of neglected soil P pools for urban land-cover was likely wrong especially for urban green areas but this assumption had no effect on cropland/grassland soil P as conversion from urban areas to agricultural areas was extremely rare.

**Drivers of agricultural soil P distribution: datasets and computation.** As explained earlier, the term “drivers” encompasses variables related to soil P input and output (farming practices - FARM, atmospheric deposition - DEPO, sludges - SLUD, losses through water erosion - LOSS), land-use change (LUCC), variables that have an effect on soil P dynamics (soil properties - SPRO, climate - CLIM), and biogeochemical background (BIOG). The different drivers and their representations were summarized in Tables 1–8. The datasets used in this study are called (data = GPASOIL-v1) while (data = GPASOIL-v0) is used to name the datasets used in Ringeval *et al.*<sup>11</sup>. Tables 1–8 provide the difference between (data = GPASOIL-v0) and (data = GPASOIL-v1) for each driver. BIOG and SPRO were constant in time while other drivers varied in time. Drivers whose the computation required more information than the one given in Tables 1–8 were described below.

**FARM.** P in manure (for both cropland and grassland) Two estimates of P in manure applied to cropland soil and two estimates for grassland soil were used to assess the uncertainty associated with this soil input (Section “Error estimate…” and Supplementary Table S1). For each land-use category (cropland or grassland), the first estimate relied on a half-degree resolution dataset for N in manure reaching the soil (Zhang *et al.*<sup>30</sup> for cropland, fN<sub>manu,Zhang</sub><sup>out→x-tot</sup>, Eq. 25, and Xu *et al.*<sup>31</sup> for grassland, fN<sub>manu,Xu</sub><sup>out→x-tot</sup>, Eq. 26). A spatially constant P:N ratio of 0.20 (kgP kgN<sup>-1</sup>) was used to convert N fluxes to P fluxes of manure.

$$fP_{manu,1}^{out→x-tot}(y, \text{crop}, g) = 0.20 * fN_{manu,Zhang}(y, \text{crop}, g) \qquad (25)$$

$$fP_{manu,1}^{out→x-tot}(y, \text{grass}, g) = 0.20 * fN_{manu,Xu}(y, \text{grass}, g) \qquad (26)$$

fN<sub>manu,Xu</sub><sup>out→x-tot</sup> correspond to the sum of N in what was called deposition by grazing animals and application of manure by Xu *et al.*<sup>31</sup>.

The second estimate was computed thanks to Eqs. 27, 28 and relied on the P in manure produced by the livestock at the country-scale given by Demay *et al.*<sup>9</sup>, which combined data on of livestock population and P excretion rates per livestock category, without any distinction between cropland and grassland soils:
$$fP_{manu,2}^{out \to x-tot}(y, \text{crop}, g) = fP_{manu,1}^{out \to x-tot}(y, \text{crop}, g) * \frac{fP_{manu,Demay}(\text{crop} + \text{grass}, \text{country}(g))}{fP_{manu,1}^{out \to x-tot}(\text{crop} + \text{grass}, \text{country}(g))} \quad (27)$$

$$fP_{manu,2}^{out \to x-tot}(y, \text{grass}, g) = fP_{manu,1}^{out \to x-tot}(y, \text{grass}, g) * \frac{fP_{manu,Demay}(\text{crop} + \text{grass}, \text{country}(g))}{fP_{manu,1}^{out \to x-tot}(\text{crop} + \text{grass}, \text{country}(g))} \quad (28)$$

with the overbar there corresponds to a temporal average over the 1950–2017 period (i.e. the years in common for Demay et al.⁹ and $fP_{manu,1}^{out \to x-tot}$), country(g) is the country where the grid-cell g can be found, $fP_{manu,Demay}^{out \to x-tot}$ is the P manure produced at country scale provided by Demay et al.⁹ and $fP_{manu,1}^{out \to x-tot}$ ($\text{crop} + \text{grass}, \text{country}(g)$) is the country average of variables computed thanks to Eqs. 25, 26 without any distinction crop vs grass. For a given country, $fP_{manu,2}^{out \to x-tot}$ has the same relative spatial distribution among grid-cells as $fP_{manu,1}^{out \to x-tot}$ but its spatial variabilities between countries is different. The ratio of 0.20 used in Eqs. 25, 26 was given by Sun et al.²⁹ and corresponds to a P:N ratio averaged among ruminant manure (mean P:N value of 0.165 between value of 0.15 for cattle and 0.18 for goat and sheep³²) and monogastric manure (mean P:N value of 0.26 between value of 0.24 for chicken and 0.28 for swine³²), weighted by the global amount of manure for each livestock super-category (ruminants: 14.4 TgN yr⁻¹ and monogastrics: 10.1 TgN yr⁻¹) given by FAOSTAT for the year 2000. The use of a spatially and temporal invariant P:N ratio for manure was needed given the lack of gridded data about P manure applied on cropland and grassland soil. It was however a simplistic assumption as N and P were partly decoupled due to N losses in manure through NH₃ volatilization. While applied to the same data of manure application in Sun et al.²⁹ as used here, the ratio of Lun et al.³² characterized manure excretion (i.e. before volatilization). Also, the P:N values used to compute the value of 0.20 was representative to a given country (USA)³² while such ratio could change as function of the country and farming practices, following its relationship with animal nutrition³³. Finally, the use of the same ratio for both cropland and grassland was also an approximation as the proportion of monogastric vs ruminant likely differs between cropland and grassland.

P in chemical fertilizer (for both cropland and grassland) Given the lack of global dataset about P fertilizer applied to grassland, a spatially constant P:N ratio of 0.22 for chemical fertilizer was used to derive P in chemical fertilizer for grassland ($fP_{chem}^{out \to i-lab}(y, \text{grass}, g)$) from N in chemical fertilizer applied to grassland ($fN_{chem,Xu}(y, \text{grass}, g)$) given by Xu et al.³¹:

$$fP_{chem}^{out \to i-lab}(y, \text{grass}, g) = 0.22 * fN_{chem,Xu}(y, \text{grass}, g) \quad (29)$$

The ratio of 0.22 corresponds to the ratio of global P chemical fertilizer for managed grassland in 2002–2010 (0.4 TgP yr⁻¹ in Lun et al.³²) and the global N chemical fertilizer for grassland for same years (1.8 TgN yr⁻¹ in Xu et al.³¹). This ratio was computed and used in Sun et al.²⁹.

P in chemical fertilizer for cropland ($fP_{chem}^{out \to i-lab}(y, \text{crop}, g)$) was derived from Lu and Tian³⁴, which assumed that all chemical fertilizer was applied to cropland (i.e. grassland did not receive any chemical fertilizer). Here, we corrected the chemical fertilizer applied to cropland by subtracting the chemical fertilizer applied on grassland given by Eq. 29 from Lu and Tian³⁴ estimates, i.e.:

for any grid-cell g and year y:

$$fP_{chem}^{out \to i-lab}(y, \text{crop}, g) = 1/\text{frac}(y, \text{crop}, g) * (fP_{chem,LuTian}(y, \text{crop}, g) * \text{frac}(y, \text{crop}, g) - fP_{chem}^{out \to i-lab}(y, \text{grass}, g) * \text{frac}(y, \text{grass}, g)) \quad (30)$$

with frac is the grid-cell fraction occupied by cropland or grassland, and $fP_{chem,LuTian}$ was taken from Lu and Tian³⁴ (in kgP ha⁻¹ yr⁻¹). This correction lead to very small change (at the global scale, $fP_{chem}^{out \to i-lab}$ over 2008–2018 is 18.1 TgP yr⁻¹ for cropland vs 0.5 TgP yr⁻¹ for grassland).

P uptake and P residues for grasslandP uptake ($fP_{upta}^{i-lab \to out}$) and P in plant residues ($fP_{resi}^{out \to i-lab}$) in grassland were derived from the actual NPP, the part of NPP which is human-appropriated through grazing/mowing (the so-called grazing intensity, GI), as well as from some spatially constant parameters as follows:

$$fP_{upta}^{i-lab \to out}(y, \text{grass}, g) = \frac{1}{\gamma} * \frac{1}{\text{dry}} * \left( \frac{P_{\%\text{above}}}{100} * \frac{\text{lifespan}_{\text{above}}}{L} * \text{NPP}_{\text{above}}(y, \text{grass}, g) + \frac{P_{\%\text{below}}}{100} * \frac{\text{lifespan}_{\text{below}}}{L} * \text{NPP}_{\text{below}}(y, \text{grass}, g) \right) \quad (31)$$

$$fP_{resi}^{out \to i-lab}(y, \text{grass}, g) = (1 - \text{GI}(y, \text{grass}, g)) * fP_{upta}^{i-lab \to out}(y, \text{grass}, g) \quad (32)$$

where γ is the carbon content of dry matter (in kgC (kgDM)⁻¹), dry is the ratio between dry-matter biomass and fresh-matter biomass (kgDM (kgFM)⁻¹), P% is the P concentration of grass biomass (in gP (100gFM)⁻¹), lifespan is the life duration of plant organs (in year), L is the length of the growing season (in years), NPP is the annual NPP (expressed here in kgC ha⁻¹ yr⁻¹), GI is the so-called grazing intensity, corresponding to the part of NPP
which is grazed/mowed (no unit). A distinction aboveground vs belowground was done for P%, lifespan, and NPP. Contrary to estimates based on carbon stock (such as used for cropland below), NPP-based estimates should consider a recycling through biomass turnover³⁵, that is why we rely here on leaf lifespan (lifespan above, in yr) and root lifespan (lifespan root, in yr) derived from observations divided by the length of the growing season (L).

Lifespan for root show large variations among climatic zones and root diameter (fine vs coarse)³⁶. Here, we used a lifespan below of 0.7 yr, close to the value provided by Zhang and Wang³⁷ for graminoids. The parameter lifespan above was set to 0.8 yr slightly below the leaf longevity used in absence of stress in LPJ-GUESS for C3-C4 grasses (1 yr in Smith et al.³⁸). We set the length of the growing season (L) equal to 1 yr. γ is set to 0.45 kgC (kgDM)⁻¹, as commonly used (see e.g Monfreda et al.³⁹). The parameter dry is set to 0.20 kgDM (kgFM)⁻¹ which is the value of dry fraction of economic yield given by Monfreda et al.³⁹ for the crop category « alfalfa ».

Terms involved in Eq. 31 are prone to uncertainty, especially NPP⁴⁰ and P concentrations. Here we used the spatial distribution of the actual NPPtot (defined as NPPtot = NPPabove + NPPbelow) given by Kastner et al.⁴¹ as they also provided the spatial distribution of GI used in the current study to estimate the P in residues. Kastner et al.⁴¹ provided different estimates and we used here one of them (the so-defined run number 1_3_2_1_1_1 in the Readme file of Kastner et al.⁴¹). In the chosen run, LPJ-GUESS without nitrogen limitation is used to approach NPPtot and its global average reaches ~460 G C m⁻² yr⁻¹, well above the value found in Wang et al.⁴² for temperate/C3 grassland (273 G C m⁻² yr⁻¹) but well below global value derived from fields (979 ± 78 G C m⁻² yr⁻¹) found in Sun et al.⁴⁰. Thus, we used for each grid-cell g, two estimates for NPPtot:

$$NPP_{tot,1}(y, \text{grass}, g) = NPP_{Kastner}(y, g) \qquad (33)$$

$$NPP_{tot,2}(y, \text{grass}, g) = NPP_{Kastner}(y, g) * \frac{NPP_{\overline{Sun}}}{NPP_{Kastner}} \qquad (34)$$

with NPP Sun and NPP Kastner equal to 979 and 460 G C m⁻² yr⁻¹, respectively. Globally averaged GI is 0.12. NPP below and NPP above are computed as 0.54*NPP tot and 0.46*NPP tot respectively, as found in Sun et al.⁴⁰.

Lun et al.⁴³ provided a value of 1.5e⁻¹ gP (100gFM)⁻¹ for the P concentration of the economic yield for the crop category « other crops » encompassing « forage and silage (maize, grasses, alfalfa, clover, ...) », which is close to the value provided by Smil⁴⁴ for the harvest of the crop category « Forages » (2.0e⁻¹ gP (100gFM)⁻¹) but much higher than the P concentration found in GOLUM⁴² (and taken from Zechmeister-Boltenstern et al.⁴⁵) or found in Spohn et al.⁴⁶. Wang et al.⁴² provided a molar P:C ratio for foliage of 1/753 for temperate C3 grass and 1/1728 for tropical C4 grass (Table S1 of Wang et al.⁴²), i.e. a C3-C4 averaged molar P:C ratio of 1/1241, very close to the molar P:C ratio of 1/1260 given by Spohn et al.⁴⁶ for senescent leaves of grassland. These values correspond to a P concentration of 2.5e⁻² gP (100gFM)⁻¹, following the below equation:

$$P_{\% above} = 100 * \gamma * \text{dry} * \frac{M(P)}{M(C)} * P_{\text{molar, above}} \qquad (35)$$

with P molar above the P:C molar ratio (no unit), P% the P concentration in gP (100gFM)⁻¹, M(P) and M(C) the atomic mass of P and C. In our study, we used this value (2.5e⁻² gP (100gFM)⁻¹) as lower boundary and the Lun et al.⁴³ value (1.5e⁻¹ gP (100gFM)⁻¹) as upper boundary, thus a mean value of 8.8e⁻² gP (100gFM)⁻¹.

Finally, we assumed the following relationship between P% above and P% below based on the order of magnitude of leaf:root ratio of molar P:C ratio used in⁴⁵ for C3 or C4 grassland:

$$P_{\% below} = \frac{1}{2} P_{\% above} \qquad (36)$$

P uptake and P residues for croplandP uptake at grid-cell scale for the cropland fraction (fP i-lab-out uptake (y, crop, g)) was computed as the sum of the uptake by different crops weighted by the harvested area of each crop:

$$fP_{upta}^{i-lab-out}(y, \text{crop}, g) = \frac{1}{\sum_{c \in \text{list}_{\text{crop}}} \text{Area}(y, c, g)} * \sum_{c \in \text{list}_{\text{crop}}} (\text{Area}(y, c, g) * fP_{upta}^{i-lab-out}(y, c, g)) \qquad (37)$$

P in plant residues (fP out→x-tot resid was computed as:

$$fP_{resi}^{out→x-tot}(y, \text{crop}, g) = \frac{1}{\sum_{c \in \text{list}_{\text{crop}}} \text{Area}(y, c, g)} * \sum_{c \in \text{list}_{\text{crop}}} (\text{Area}(y, c, g) * fP_{resi}^{out→x-tot}(y, c, g)) \qquad (38)$$

with c a crop within the list of crop considered in Monfreda et al.³⁹ (called here list crop) and Area(y,c,g) the area covered by the crop (in ha). Among 172 crops listed in list crop, 14 were not considered as we did not find any corresponding crop in FAOSTAT. It is particularly the case for some crops of the group « forage » in Monfreda et al.³⁹ (alfalfa, beets for fodder, clover, rye grass for forage and silage, etc.) but not for all, as we made some corresponding between non-forage and forage crop in FAOSTAT. E.g. we used the “maize” crop in FAOSTAT to get some information needed to be combined with the crop “maize for forage” of Monfreda et al.³⁹. The corresponding between non-forage and forage crop was disputable but it concerned here only the relative change in yield and area as compared to the year 2000 (Eqs. 43, 46). We also excluded carob as carob area in Spain shows very weird pattern in FAOSTAT (from ~130 000 ha on average over 1961–1989 to 450 ha on average over 1990–2017).
For each crop, P uptake and P residues are computed as follows:

$$fP_{uptake}^{i \to lab \to out}(y, c, g) = \frac{\text{Yield}(y, c, g)}{100} \cdot \left[ P_{\%_{\text{root}}}(c) \cdot \frac{\text{dry}_{\text{harvest}}(c)}{\text{dry}_{\text{root}}(c)} \cdot \frac{\text{RSR}(c)}{\text{HI}(c)} + P_{\%_{\text{abov}\text{harvest}}}(c) \cdot \frac{\text{dry}_{\text{harvest}}(c)}{\text{dry}_{\text{abov}\text{harvest}}(c)} \cdot \left( \frac{1}{\text{HI}(c)} - 1 \right) + P_{\%_{\text{harvest}}}(c) \right] \quad (39)$$

$$fP_{\text{resid}}^{\text{out} \to x-\text{tot}}(y, c, g) = \frac{\text{Yield}(y, c, g)}{100} \cdot \left[ P_{\%_{\text{root}}}(c) \cdot \frac{\text{dry}_{\text{harvest}}(c)}{\text{dry}_{\text{root}}(c)} \cdot \frac{\text{RSR}(c)}{\text{HI}(c)} + \text{frac}_{\text{resid}}(c) \cdot P_{\%_{\text{abov}\text{\textbackslash} \text{harvest}}}(c) \cdot \frac{\text{dry}_{\text{harvest}}(c)}{\text{dry}_{\text{abov}\text{\textbackslash} \text{harvest}}(c)} \cdot \left( \frac{1}{\text{HI}(c)} - 1 \right) \right] \quad (40)$$

with P%_m the P concentration (in gP (100gFM)⁻¹) for the organ m with m being {root, abov\harvest, harvest} where « abov\harvest » defines the aboveground biomass excluding the harvest (also called total residues in Smil⁴⁴), dry_m the dry proportion of organ m expressed in fresh matter (in gDM gFM⁻¹), frac_resid the fraction of the aboveground biomass excluding the harvest that remains on the soil (the rest being exported from the field), Yield the economic yield (expressed in kgFM ha⁻¹), RSR the root:shoot ratio (no unit), HI the harvest index (no unit). Text S2 shows how we get these equations as well as how we get corresponding equations for crops whose the harvested part being root. Following Bentsen et al.⁴⁷, we assumed that at harvest, the different organs have the same dry proportion (i.e. dry_grain = dry_abov\harvest = dry_root). For all crops but root crops, we considered P%_abov\harvest equal to the P concentration of aboveground residues given by Lun et al.⁴³ (and completed for few crop categories by Smil⁴⁴). P%_harvest is provided by Lun et al.⁴³. RSR is derived from the aboveground fraction (frac_above) given in Monfreda et al.³⁹ (RSR = 1/frac_above - 1) and HI is directly given by Monfreda et al.³⁹. frac_resid is equal to 0 for forage crop (group « forage » in Monfreda et al.³⁹) and equal to 1/2 for other crops, following Smil⁴⁴. We did not find any database available providing the P concentration of roots for different crop categories. Instead we assumed, following Ye et al.⁴⁸ (focusing on rice), that for any crop c:

$$P_{\%_{\text{root}}}(c) = 0.75 \cdot P_{\%_{\text{abov}\text{\textbackslash} \text{harvest}}}(c) \quad (41)$$

The specificities of root crop computation are given in Text S2.

The spatially explicit yield (Yield (y, c, g) in kgFM ha⁻¹) for the crop category c is computed as follows:

$$\text{Yield}(y, c, g) = \text{Yield}_{\text{Monfreda}}(2000, c, g) \cdot \frac{\text{Yield}_{\text{FAO,bis}}(y, c, \text{country}(g))}{\text{Yield}_{\text{FAO}}(2000, c, \text{country}(g))} \quad (42)$$

with Yield_Monfreda the spatially explicit (half-degree resolution) distribution of yield given by Monfreda et al.³⁹ for years 2000 (and expressed as weight of economic yield, i.e. with standardized water content), Yield_FAO the country-scale FAO yield and Yield_FAO,bis defined as follows:

$$\text{Yield}_{\text{FAO,bis}}(y, c, \text{country}(g)) = \text{Yield}_{\text{FAO}}(y, c, \text{country}(g)) \quad \text{if } y \ge 1961 \quad (43)$$

$$\text{Yield}_{\text{FAO,bis}}(y, c, \text{country}(g)) = \text{Yield}_{\text{FAO}}(1961, \text{country}(g), c) \cdot \frac{\text{pop}(y, \text{country}(g))}{\text{pop}(1961, \text{country}(g))} \quad \text{if } y < 1961 \quad (44)$$

with pop the country-scale population defined by HYDE 3.2²³. With Eq. 42, we combined the spatially explicit yield distribution of Monfreda et al.³⁹ (but only one value for year 2000) with the temporal varying FAOSTAT yield (but at country-scale). For each crop, the relative intra-country distribution of yield given by Monfreda et al.³⁹ for year 2000 was kept constant in time. With Eq. 44, we scaled the country-scale yield to the country-scaled human population before 1961, as it was assumed in Bouwman et al.¹³. As country spatial boundaries can change in time, we made some computations to allow consistency between each grid-cell and its belonging to a country whose the boundaries can change in time (Text S3).

The harvested area (Area (y, c, g) in ha) for the crop category c for any grid-cell g and any year y is computed as follows:

$$\text{Area}(y, c, g) = \text{area}_{\text{Monfreda}}(2000, c, g) \cdot \frac{\text{Area}_{\text{FAO,bis}}(y, c, \text{country}(g))}{\text{Area}_{\text{FAO}}(2000, c, \text{country}(g))} \quad (45)$$

with Area_Monfreda the spatially explicit distribution of harvested area given by Monfreda et al.³⁹ for years 2000, Area_FAO the country-scale FAO harvested area and Area_FAO,bis defined as follows:
$$Area_{FAO,bis}(y, c, country(g)) = Area_{FAO}(y, c, country(g)) \text{ if } y > 1961 \quad (46)$$

and

$$Area_{FAO,bis}(y, c, country(g)) = Area_{FAO}(1961, c, country(g)) \text{ if } y < 1961 \quad (47)$$

Composition of manure and residues for cropland and grassland The composition (inorganic labile, organic labile, organic stable) of manure and residues were considered constant in space and among cropland/grassland:

$$fP_{manu}^{out \to i-lab}(y, lu, g) = frac_{i-lab,manu} * fP_{manu}^{out \to x-tot}(y, lu, g) \quad (48)$$

$$fP_{manu}^{out \to o-lab}(y, lu, g) = frac_{o-lab,manu} * fP_{manu}^{out \to x-tot}(y, lu, g) \quad (49)$$

$$fP_{manu}^{out \to o-sta}(y, lu, g) = (1 - frac_{i-lab,manu} - frac_{o-lab,manu}) * fP_{manu}^{out \to x-tot}(y, lu, g) \quad (50)$$

$$fP_{resi}^{out \to i-lab}(y, lu, g) = frac_{i-lab, resi} * fP_{resi}^{out \to x-tot}(y, lu, g) \quad (51)$$

$$fP_{resi}^{out \to o-lab}(y, lu, g) = frac_{o-lab, resi} * fP_{resi}^{out \to x-tot}(y, lu, g) \quad (52)$$

$$fP_{resi}^{out \to o-sta}(y, lu, g) = (1 - frac_{i-lab, resi} - frac_{o-lab, resi}) * fP_{resi}^{out \to x-tot}(y, lu, g) \quad (53)$$

with $lu$ in {crop,grass}, $y$ the year and $g$ the grid-cell considered, $frac_{m,l}$ corresponding to the fraction of $fP_{manu}^{out \to x-tot}$ reaching the pool $m$. Following Ringeval et al.¹¹, we prescribed: $frac_{i-lab,manu} = 0.8$, $frac_{o-lab,manu} = 0.1$, $frac_{i-lab, resi} = 0.4$ and $frac_{o-lab, resi} = 0.4$. We considered an uncertainty associated to the composition of residues. To do so, we use a random value between $frac_{i-lab, resi}$ - 50% and $frac_{i-lab, resi}$ + 50% as well as as between $frac_{o-lab, resi}$ - 50% and $frac_{o-lab, resi}$ + 50% (see the section about error estimates).

DEPO. Different classes of atmospheric deposition of P are often considered based on the sources of P in the atmosphere: mineral dust, sea salt, primary biogenic aerosol particles (PBAP), natural combustion and anthropogenic combustion⁴⁹. We considered that P deposition input enriches both the inorganic labile pool ($fP_{depo}^{out \to i-lab}$) and the primary inorganic pool ($fP_{depo}^{out \to i-prim}(y, lu, g)$), as follows:

$$fP_{depo}^{out \to i-lab}(y, lu, g) = frac_{i-lab,dust} * D_{dust}(y, g) + frac_{i-lab,other} * (D_{seasalt}(y, g) + D_{PBAP}(y, g) + D_{natcomb}(y, g) + D_{anthcomb}(y, g)) \quad (54)$$

$$fP_{depo}^{out \to i-prim}(y, lu, g) = (1 - frac_{i-lab,dust}) * D_{dust}(y, g) + (1 - frac_{i-lab,other}) * (D_{seasalt}(y, g) + D_{PBAP}(y, g) + D_{natcomb}(y, g) + D_{anthcomb}(y, g)) \quad (55)$$

with $D_{dust}$, $D_{seasalt}$, $D_{PBAP}$, $D_{natcomb}$ and $D_{anthcomb}$ are atmospheric P deposition fluxes (in kgP ha⁻¹ yr⁻¹) of mineral dust, sea salt, PBAP, natural combustion and anthropogenic combustion, respectively. $frac_{i-lab,dust}$ and $frac_{i-lab,other}$ are labile fractions of atmospheric deposition of P (no unit) for respectively dust and all other sources (sea salt, PBAP, natural and anthropogenic combustion). $frac_{i-lab,dust}$ was chosen to be equal to 0.1 while $frac_{i-lab,other}$ was chosen equal to 0.5⁵⁰.

For any year $y$ of the simulation and for any grid-cell $g$, we computed $D_{dust}(y, g)$, $D_{seasalt}(y, g)$, $D_{PBAP}(y, g)$, $D_{natcomb}(y, g)$ and $D_{anthcomb}(y, g)$ based on variables provided by Wang et al.⁴⁹,⁵¹. Wang et al.⁴⁹ provided, for any grid-cell $g$, $D_{dust}^A$ (2000 - 2011, $g$), $D_{seasalt}^A$ (2000 - 2011, $g$), $D_{PBAP}^A$ (2000, $g$), $D_{totcomb}^A$ (1960 - 2007, $g$) and for any World regions reg, $E_{anthcomb}^A$ (1960 ≤ $y$ ≤ 2007, reg) and $E_{natcomb}^A$ (1960 ≤ $y$ ≤ 2007, reg) while Wang et al.⁵¹ provided, for any grid-cell $g$, $D_{all}^B$ (1997 ≤ $y$ ≤ 2013, $g$), with $D$: atmospheric deposition of P, $E$: P emissions to the atmosphere, the upper letter corresponding to the version of the dataset of deposition (A: Wang et al.⁴⁹ and B: Wang et al.⁵¹), the line over years means that only the data for an averaged time-period is available while $a \le y \le b$ means that data are available for each year between $a$ and $b$, $totcomb$ means ($natcomb + anthcomb$) while all means: ($dust + (seasalt + PBAP + natcomb + (anthcomb)$.

To summarize our computation, the strategy was to consider $D_{dust}$, $D_{seasalt}$, $D_{PBAP}$, $D_{natcomb}$ as static in time, and temporal average given by Wang et al.⁴⁹ were used for these variables. $D_{natcomb}$ (respectively $D_{anthcomb}$) were derived from $D_{totcomb}$ and the ratio between natcomb and totcomb (resp. anthcomb and totcomb) in emissions. $D_{anthcomb}(y, g)$ was considered varying in time and its interannual variability was estimated from the interannual variability of emissions from anthropogenic combustion. For years after 2007, Wang et al.⁵¹ were used and corrected to ensure equality with Wang et al.⁴⁹ for the years in common (1997–2007). Equations are given in Text S4.
SLUD. P in sludges from sewage treatment that reaches cropland soils P ($fP_{slud}^{out→x-tot}$ in kgP ha⁻¹ year⁻¹) is estimated by combining human P excretions with the fraction of sewage sludge that is treated and the P removal efficiency of treatments plants following van Puijenbroek et al.⁵² and Demay et al.⁹:

$$fP_{slud}^{out→x-tot}(y, \text{crop}, g) = \sum_{\text{type}=1,2,3} (\text{frac}_{\text{treat}}(y, \text{country}(g), \text{type}) * \text{NR}(\text{type})) * \text{human}_{\text{excr}}(y, \text{reg}(\text{country}(g))) * \text{pop}(y, \text{country}(g)) / \text{Area}(y, \text{crop}, \text{country}(g)) \quad (56)$$

with human$_{excr}$ the total P in human excretion (in kgP capita⁻¹ year⁻¹) for big World regions, frac$_{treat}$ is country-scale fraction of total population (no unit) with access to different types of wastewater treatment installations (types 1, 2 or 3), each one being characterized by a different removal fraction of P (NR, no unit with NR = 0.10, 0.45, 0.90 for respectively wastewater treatment installations type 1, 2, 3), pop the country-scale human population (in inhabitants), Area the country-scale cropland area (in ha).

NR, human$_{excr}$, frac$_{treat}$ are provided by van Puijenbroek et al.⁵², pop is provided by HYDE 3.2²³ and Area at country scale is computed from the dataset used in LUCC. The variable human$_{excr}$ is available for the years 1970 and 2010. The variable frac$_{treat}$ is available for years 1990, 2000, 2010 in van Puijenbroek et al.⁵² and we computed it for the year 1970 (see Text S5 for further technical details). Then, linear regressions are used to derive human$_{excr}$ and frac$_{treat}$ for any year between 1970 and 2010. Values for year 1970 (resp. 2010) were used for any years before (resp. after) 1970 (resp. 2010). The corresponding between big World regions and countries was provided by the IMAGE framework region classification (https://models.pbl.nl/image/index.php/Region_classification_map).

The Eq. 56 relies on the following assumptions:

- P in sludges is solely sourced from human food consumption (thereby excluding P release from detergents)

- all sludges produced are used in agriculture. This is a strong assumption as a few countries have specific rules banning the use of sludge in agriculture (e.g. Switzerland incinerated sludges⁵³), but this assumption was required as we did not find any database that compiled this information per country at the global scale

- all sludges used in agriculture are spread on croplands (not on grasslands).

For the composition of P in sludges, we used the same labile vs. stable contribution as for manure without considering any uncertainty related to the composition (Supplementary Table S1).

LOSS. The annual lateral flux of P lost from soil erosion by water (in kgP ha⁻¹ yr⁻¹) was computed by combining the fraction of soil mass of the top 0–0.3 m layer lost each year through erosion (frac$_{loss}$, in yr⁻¹) with each soil P pools (kgP ha⁻¹):

$$fP_{loss}^{m→out}(y, lu, g) = P_m(y, lu, g) * frac_{loss}(y, lu, g) \quad (57)$$

with m in {i-lab, o-lab, o-sta, i-prim, i-sec, i-occ}, lu in {crop, grass}, and P$_m$ the soil P content of pool m simulated by our approach.

frac$_{loss}$ was computed by combining the gross soil losses by water erosion with a cropland/grassland distinction (loss, in (kg of soil) ha⁻¹ yr⁻¹), the soil bulk density of the fine earth fraction (bdod) and the volumetric fraction of coarse fragment (cfvo) to estimate the soil mass of the top 0–0.3 m layer:

$$frac_{loss}(y, lu, g) = \frac{loss(2000, lu, g)}{bdod(g) * (1 - cfvo(g)) * \Delta z * 10^{+4}} \quad (58)$$

with bdod the bulk density for the top 0–0.3 m soil layer (in (kg of fine soil) m⁻³), cfvo is the volumetric fraction of coarse fragment (without unit), Δz = 0.3 m and 10⁺⁴ a converting factor in m² ha⁻¹. The variable loss was computed following the approach described in Borrelli et al.⁵⁴ but here applied to the land-use dataset described for year 2000 in LUCC. The geo-statistical approach proposed by Borrelli et al.⁵⁴ allowed for an accurate spatial definition of the land uses (native resolution at 250 m cell size at the equator aggregated to 0.5° latitude × 0.5° longitude) and the consideration of the effects of the different regional cropping systems.

As in Ringeval et al.¹¹, Eq. 58 is based on the following assumptions. First, we assumed that the soil losses of Borrelli et al.⁵⁴ corresponded to a loss of topsoil (i.e. 0–0.3 m soil layer). Second, the estimates of losses for year 2000 (in (kg of soil) ha⁻¹ yr⁻¹) were considered representative of erosion rates for the whole century. Last, the estimates from Borrelli et al.⁵⁴ correspond to a gross erosion rate and not to the net flux, which results from both mobilization and processes, such as deposition, storage and burial. However, we considered that these later processes were not relevant to the system represented here, either because they largely occur in non-agricultural fractions of the grid-cell (e.g. deposition on river banks) or because they concerned soil horizons below 0.3 m (burial).

LUCC. The following land-use categories are considered in our approach: cropland (crop), grassland (grass), non-agricultural vegetation (nonagri), and urban (urban). For each year y and grid-cell g, we needed the fractions of cropland (frac(y, crop, g)) and of grassland (frac(y, grass, g)) as well as the transitions between the above mentioned 4 land-use categories (i.e. Δᵢ(y, j, g)) for the conversion from j to i with both i and j in {crop, grass, nonagri, urban}). These variables were estimated based on Chini et al.²⁰ after we made the corresponding between our land-use categories (left hand side) and the ones considered in Chini et al.²⁰ (right hand side):
$$ cropland = c3ann + c3nfx + c3per + c4ann + c4per \quad (59) $$

$$ grassland = pastr + range \quad (60) $$

$$ nonagri = primf + primn + secdf + secdn \quad (61) $$

$$ urban = urban \quad (62) $$

with the following meaning for categories in Chini et al.²⁰: c3ann: C3 annual crops, c3nfx: C3 nitrogen-fixing crops, c3per: C3 perennial crops, c4ann: C4 annual crops, c4per: C4 perennial crops, past: managed pasture, range: rangeland, primf: forested primary land, primn: non-forested primary land, secdf: potentially forested secondary land, secdn: potentially non-forested secondary land, urban: urban land.

**CLIM.** We used the simulations performed by land-surface models for the CMIP6 exercise⁵⁵ to get information about the soil liquid water content ($W_{rel}$, in fraction of the field capacity and $W_{abs}$, in L(kg of soil)$^{-1}$) and the soil temperature ($T_{soil}$ in °C) for the top soil layer considered in our approach (0–0.3 m) over the historical period. All CMIP-6 simulations providing the following variables at the basis of our computation of $W_{reb}$, $W_{abs}$ and $T_{soil}$ have been used: *mrsll* (liquid water content of soil layers, in kg.m$^{-2}$), *mrsofc* or *fldcapacity* (field capacity, in kg m$^{-2}$), *tsl* (soil temperature, in °C). The choice of *mrsofc* or *fldcapacity* to approach the field capacity of a given land-surface model depends if the field capacity varies with depth in this land-surface model. Finally, the average among 9 simulations (combination between 4 land-surface models: CNRM-ESM2-1, CNRM-CM6-1, IPSL-CM6A-LR, MIROC6; and 3 climate data used as input of the land-surface models for the historical period: land-hist, land-hist-cruNcep, land-hist-princeton; see Table 6 for the reference of each simulation) has been computed and used. The period 1850–2012 was common to the 9 simulations. For years after 2012, to prevent any bias due to the lack of some simulations, we used the 1993–2012 climatology each year.

The annual average of each variable ($W_{abs}$, $W_{rel}$ and $T_{soil}$) was used in our approach. The influences of $T_{soil}$ and $W_{rel}$ on weathering and mineralization were modelled using different functions ($h_1$ and $g_1$ for sensitivity to $T_{soil}$, $h_2$ and $g_2$ for sensitivity to $W_{rel}$) in both (model = GPASOIL-v0) and (model = GPASOIL-v1). $W_{abs}$ was used in (model = GPASOIL-v1) in some parametrizations (Eqs. 12, 13) and to translate $P_{i-soil}$ into soil solution P concentration (Eqs. 3, 7). In addition, the near-surface atmospheric temperature needed in one parametrization of (model = GPASOIL-v1) ($T_a$ in Table 10) was taken from climate data sets used as input of the land-surface models quoted above.

**Error estimate from uncertainty associated with the datasets to describe the different drivers.** While we recognized the existence of different sources of uncertainty (dataset, model structure, parametrization⁵⁶), we mainly focus here on how the uncertainty related to the representation of the different drivers propagated to the simulation of current cropland/grassland soil P pools. Nevertheless, the uncertainty associated to SPRO has an effect on the model parameterizations. The uncertainty associated to mineralization rates (model parameter) was investigated separately.

To investigate the propagation of the uncertainty associated with the driver datasets, we performed 100 simulations by using a random value for each variable related to a given driver (see the list of variables in the 3rd line of Tables 1–8) with all grid-cells considered as independent. For each driver, how the random value was defined was described in Supplementary Table S1 and summarized here. The uncertainty associated with the driver LUCC was not considered given the difficulty to properly quantify it by keeping consistency between the different variables involved in LUCC (fraction and transitions between the different land-use). The coefficient of variation computed over the 100 simulations was used hereafter as uncertainty estimate.

For a given driver, when two estimates are considered (called hereafter E1 and E2), we choose to use a normal distribution with $4\sigma$ between E1 and E2, which means that ~96% of the random values chosen will fall between the two estimates. This means that we considered E1 and E2 as lower/upper boundaries and we considered the range they defines as encompassing most values possible. Instead of a normal distribution, a uniform distribution was considered for variable that we considered very uncertain (P in manure for cropland and grassland, residues composition) as we considered each estimate (E1 or E2) as reasonable as the other estimate. For a few variables related to some drivers for which only one estimate was available (called hereafter E) and no uncertainty was provided (i.e. P in chemical fertilizer for grassland, P uptake for cropland, residues composition, P in sludge), we derived E1 and E2 as $(1 - x)*E$ and $(1 + x)*E$ with x either equal to 30% or 50% depending on how we consider the driver uncertain, i.e. a standard-deviation equal to 15% or 25% of the mean, respectively.

Random values can lead to unmeaning values for few grid-cells. To prevent this, we set the minimum value allowed for soil input/output to zero. For variables related to other drivers (BIOG, CLIM, SPRO), we prescribed the minimum/maximum values by using the global spatial minimum/maximum of the mean value (Supplementary Table S1). For temporal varying drivers, we consider the uncertainty related to the different time-steps as independent.

In the following, we focused on the simulation performed with the mean value for each driver instead of the mean among the 100 simulations performed to assess the uncertainty. Both are close except for grid-cells with very small soil P pools. In these grid-cells, if the combination of random soil input/output lead to negative values of $P_{i-lab}$, $P_{i-lab}$ is set to 0 and thus, the mean of the 100 simulations was artificially increased, leading to values
larger than the simulation performed with mean drivers. As consequence, the coefficient of variation (instead of the standard deviation) among the 100 simulations was used to assess the uncertainty in soil P pools.

One additional set of 100 simulations was performed to estimate the role played by the uncertainty in BIOG on the uncertainty of P<sub>i-lab</sub>. In this set, we used the mean value for all variables related to BIOG and considered uncertainty for all other drivers. Difference between the CV computed with this set of 100 simulations and the one computed in the set where all drivers are uncertain allows to distinguish the role played by BIOG.

## Data Records

One tgz file called GPASOIL-v1_output.tgz was made available on recherche.data.gouv.fr<sup>57</sup> and contains the soil P pools simulated with (data = v1) and (model = v1.1), as described earlier, in a netcdf format. These soil P pools correspond to GPASOIL-v1.1 (cf. Table 11). To assess the uncertainty related to mineralization rates, we also provided the GPASOIL-v1.0 soil P pools, i.e. the pools simulated with (data = v1) and (model = v1.0). Concerning the coordinate systems, all variables are referenced to the World Geodetic System (WGS84). The cell size is 0.5 decimal degrees. Fluxes of P corresponding to soil P input/output and to soil P dynamics are also given in the same file. Output are provided for the period 1900-2018 at annual time-step for both cropland and grassland. Fluxes were expressed in kgP ha<sup>-1</sup> yr<sup>-1</sup> and pools were expressed in kgP ha<sup>-1</sup>, and are both representative to the top soil layer 0-0.3 m. Two files are provided: one corresponding to the simulation with mean estimate for each driver and another one corresponding to the coefficient of variation to assess the effect of the uncertainty in driver estimate on the simulated soil P pools. README_output_upload.txt describes the netcdf files. An additional netcdf file with the land area was also provided to allow users to compute global averages of soil P pools.

## Technical Validation

Our technical validation evaluated first each component of our approach, namely i) the driver estimates based on published datasets and ii) the model of soil P dynamic. Then it focused on how the update of each component (dataset, model) modified the soil P maps simulated from GPASOIL-v0 to GPASOIL-v1. This is part of the technical validation as it allows us to trace the impact of the update of each component on the reconstructed soil P. In this section, different combinations (data x model) were used (Table 11). Then, we discussed about how our best estimates (GPASOIL-v1.1) compared to previous published modelling products Zhang *et al.*<sup>12</sup> (McDowell *et al.*<sup>58</sup> and GPASOIL-v0). Finally, we compared the spatial variability simulated by GPASOIL-v1.1 to the one of published regional datasets on measured soil P.

**Evaluation of the computation of some selected drivers.** The estimate of the drivers considered in our approach (BIOG, FARM, etc. see Tables 1–8) was based on published datasets. While the estimate was straightforward from published datasets for some drivers (e.g. BIOG), some additional computations were required for others, in particular for P uptake (with strong differences in the computation between cropland and grassland), P in manure and P deposition. In the current section, we checked that our own computation performed well, i.e. that we succeed in applying the equations given in the section named "Drivers of agricultural soil P distribution: datasets and computation". If possible, we compared some intermediary variables (e.g. yield per crop category, which is an intermediary variable of the computation of P uptake for cropland) or final P-related variables (e.g. P in manure) to other independent datasets.

In Fig. 2, we compared our estimates of P in biomass grazed/mowed from grassland to the estimate computed by Demay *et al.*<sup>9</sup>. P in biomass grazed/mowed from grassland was defined as the difference $f_{upta}^{i-lab→out}(grass) - f_{resi}^{out→x-tot}(grass)$, or following Eq. 32 as $GI * f_{upta}^{i-lab→out}(grass)$. The two estimates used to assess the uncertainty ($GI * f_{upta,1}^{i-lab→out}(grass)$, $GI * f_{upta,2}^{i-lab→out}(grass)$) and their mean ($GI * f_{upta,mean}^{i-lab→out}(grass)$) were plotted in Fig. 2. P in biomass grazed/mowed from grassland in Demay *et al.*<sup>9</sup> (called P harvest from grassland in this latter reference) was computed at country scale based on forage P demand of livestock (Eq. S4 of Demay *et al.*<sup>9</sup>). We found that the independent estimate of Demay *et al.*<sup>9</sup> was very close to the mean of our estimates (black curve in Fig. 2).

P uptake for cropland was based on the spatially explicit temporal yield per cropland that we estimated using Eq. 42. Figure S1 shows how the Eq. 42 provided temporal estimates of yield for maize at grid-cell scale (Yield in kgFM ha<sup>-1</sup>, line 2 of Fig. S1) by combining Yield<sub>Monfreda</sub>, the spatially explicit yield given by Monfreda *et al.*<sup>39</sup> for circa the year 2000 (line 4) with Yield<sub>FAO</sub>, the country-scale annual FAO yield (line 3). Within each country, the relative spatial variability was kept constant for any year and equal to one of Yield<sub>Monfreda</sub>. Figure S1 focuses on the temporal averages over three time-periods (1982–1992, 1993–2003, 2004–2014) as it allows to compare our estimates with estimates of Iizumi and Sakai<sup>59</sup> (line 1) available over the 1982–2014 time-period. Overall, a relatively good agreement is found despite some discrepancies between our estimate and<sup>59</sup>, explained by some mismatch in country-scale variability (e.g. USA) or global distribution of the crop area (e.g. India) between Monfreda *et al.*<sup>39</sup> and Iizumi and Sakai<sup>59</sup>. Note that Iizumi and Sakai<sup>59</sup> assumed no temporal change in harvested area and focused only on few crops, that is why we only used this dataset for comparison in the current study.

Figure S2 focuses on the temporal variation at country scale of Yield and Area, the spatially explicit and crop-specific yield and harvested area, respectively. Yield and Area are computed thanks to Eqs. 42, 45. As Fig. S1, Figure S2 focuses on maize (c = maize in Eqs. 42, 45). We compared Yield to the dataset of Monfreda *et al.*<sup>59</sup> and Area to the dataset used to characterize LUCC<sup>20</sup> for the broad crop category encompassing maize in LUCC (namely “c4ann” in Chini *et al.*<sup>20</sup>). LUCC is used in our study but not within the computation of the cropland P uptake and thus is independent to its computation. Figure S2 shows that our procedure allows both Area and Yield to match Monfreda *et al.*<sup>39</sup> at country-scale for the year 2000 and to mimic the year-to-year variability of FAOSTAT (left columns of Fig. S2). When we put together all crops belonging to the category « c4ann » of
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Comparison at country-scale between P in biomass grazed/mowed from grassland estimated in this study and estimates provided by Demay *et al.*⁹. P in biomass grazed/mowed from grassland is defined as GI*fPuptake* in our study. Two estimates used in our study are plotted (subscript 1 in green, subscript 2 in blue) as well as the mean of these two estimates (subscript 'mean' in black). P in biomass grazed/mowed from grassland in Demay *et al.*⁹ (called P harvest from grassland in this reference) is plotted in orange and is an independent estimate. All fluxes are expressed in kgP (ha of grassland)⁻¹. Estimate 1 (respectively estimate 2) corresponds to the min (resp the max) value for all grid-cells within the country considered. “World” (panel (g)) corresponds to the sum of all countries considered in both our estimates and Demay *et al.*⁹ and thus excludes few countries not available in Demay *et al.*⁹. The number in the left-top corner of each panel corresponds to the root mean square error (kgP (ha of grassland)⁻¹) between orange and black curves for years in common.</figcaption></figure>

LUCC, we found that the simulated year-to-year variability of area differ to the one given by LUCC (except for Brazil) but that the order of magnitude of area per country relatively matches the one of LUCC. Crop composition of total P uptake is given in Fig. S3.

In our approach, P in manure applied to cropland and grassland soils was estimated following two different methods. While the 1st estimate was based on half-degree resolution dataset for N in manure and N:P ratios, the 2nd estimate was based on country-scale estimates of Demay *et al.*⁹ with the same inner-country spatial distribution and temporal variability as the 1st estimate. Figure 3 shows the estimate 1 and estimate 2 for few countries and at the global scale. Demay *et al.*⁹ did not distinguish cropland from grassland and we plotted in Fig. 3 the sum of P in manure for the two land-covers. By construction, the country-averages of estimate 2 and Demay *et al.*⁹ over 1950–2017 are equal, and the temporal variability of estimate 2 is similar to the one of estimate 1.
Figure 3 shows that estimates 1 and 2 greatly differed, with, e.g. 1.4 kgP/ha for estimate 1 vs 3.9 for estimate 2, over 1961–2017 at global scale (Fig. 3g). On Fig. 3, we also compared our estimates to the following independent country-scale computation ($fP_{manu,indirect}$, in kgP (ha of cropland + grassland)$^{-1}$) based on a livestock budget:

$$fP_{manu,indirect}(y, \text{crop} + \text{grass}, \text{country}) = \frac{1}{\text{Area}(y, \text{crop} + \text{grass}, \text{country})} \cdot \left[ \begin{aligned} & \text{GI}(y, \text{grass}, \text{country}) \\ & \quad \cdot fP_{upta,mean}^{i\text{-lab}\to\text{out}}(y, \text{grass}, \text{country}) \\ & \quad \cdot \text{Area}(y, \text{grass}, \text{country}) \\ & \quad + \text{grain}_{feed}(y, \text{country}) \\ & \quad + \text{miner}_{feed}(y, \text{country}) \\ & \quad - \text{product}_{animal}(y, \text{country}) \end{aligned} \right] \quad (63)$$

with $GI(y, grass, country) \cdot fP_{upta,mean}^{i\text{-lab}\to\text{out}}(y, grass, country) \cdot \text{Area}(y, grass, country)$ the P in biomass withdrawn from grassland (either mowed or grazed) to feed the livestock (in kgP), $\text{grain}_{feed}$ and $\text{miner}_{feed}$ the grain and mineral used to feed livestock, respectively, and product$_{animal}$ the P in animal products (P in milk, egg and slaughtered animals). The variables $\text{grain}_{feed}$, $\text{miner}_{feed}$ and product$_{animal}$ are expressed in kgP. The country-scale data required to compute $\text{grain}_{feed}$, $\text{miner}_{feed}$ and product$_{animal}$ were provided by FAOSTAT and Demay et al.⁹. Parameters needed to compute product$_{animal}$ (average carcass weight, carcass yield over live weight, P content of animal live weight, P concentration of egg and milk) were found in Senthilkumar et al.⁶⁰. We found that $fP_{manu,indirect}$ was close to our mean estimate in both temporal average value (2.6 kgP ha$^{-1}$ for our mean estimate over 1961–2018 at global scale vs 2.2 for $fP_{manu,indirect}$) and temporal correlation (e.g. pearson correlation of 1.0, p < 0.05). Interestingly, the year-to-year variability of our estimates and the one of Demay et al.⁹ show some consistent patterns for China, India and Russia but less for other countries.

Figure S4a shows that our procedure to compute P deposition allowed to get a consistent global P deposition from anthropogenic combustion over the whole time-period.

**Application of the Wang et al. soil P inorganic dynamics parametrizations at the global scale.** (model = GPASOIL-v1) was based on Wang et al.¹⁸ parametrization for soil inorganic P dynamics. Parameters used in our approach follow the relationships with pedo-climatic variables found in Wang et al.¹⁸ with some modifications (all k parameters constant in time, oxalate excluded from parametrization leading to both simplified equation for k parameters and $P_{c,\infty}$ constant in space). Soil P fractions used in the relationships to compute the different parameters were given by the initial conditions chosen for the soil P pools (i.e. $P_{c,\infty} = 0.1$ mgP L$^{-1}$ for $P_{i-sol}$ and P of unmanaged soil P given by He et al.¹⁴ for other pools). Figure 4 shows the upper value behind which the net flux is larger than the pool it leaves at daily time-step. The issue was particularly true for $k^{i\text{-sol}\to i\text{-sec}}$ whose the global median in (model = GPASOIL-v1) (12.6 mgP (kg soil)$^{-1}$ day$^{-1}$ (mg P/L)$^{-b}$) is well above the value consistent with a daily time-step ($W_{abs} (P_{c,\infty})^{1-b}$ has a median value of 0.05 mgP (kg soil)$^{-1}$ day$^{-1}$ (mg P/L)$^{-b}$) (Fig. 4). As an equilibrium was assumed between $P_{i-sol}$ and $P_{i-sec}$, this issue does not exist for neither $k^{i\text{-sol}\to i\text{-lab}}$ nor $k^{i\text{-lab}\to i\text{-sol}}$. To solve this issue, a shorter time-step should be chosen, but the time-step needed (e.g. ~5 minutes based on median values given above) is not consistent with global application. Our strategy was to not put any constraints on the parameters themselves nor on each gross flux but we prevent a net flux to be larger than the pool from which this net flux leaves.

Figure S5 shows the sensitivity to the time-steps chosen (1-day or 6-hours time-steps) on a given grid-cell. In such sensitivity tests, the flux parameters were modified by multiplying them by the ratios of the different time-steps. The effect is particularly visible on time-steps with positive soil P budget (Figs. S5, S6) for which the inconsistency between a long time-step (1 day) and $k^{i\text{-sol}\to i\text{-sec}}$ reduces the net flux from $P_{i-sol}$ to $P_{i-sec}$, letting the soil P input within $P_{i-lab}$. The effect of the time-step chosen could be significant at grid-cell scale (e.g. + 17% larger $P_{i-lab}$ for 1-day time-step than for 6-hours time-step in 2010, panels (b) and (c) of Fig. S5) or regionally (−15% lower $P_{i-lab}$ in China for 1-day time-step than for 2-days time-step, Fig. S6a) but remains small at the global scale (few percent of difference in the global average of both $P_{i-sec}$ and $P_{i-lab}$ between 1-day time-step vs 2-days time-step, Fig. S6).

**Effect of the update of each component (data and model) from GPASOIL-v0 on $P_{i-lab}$ simulated in 2018.** Effect of the update of the dataset to describe the drivers (from (data = GPASOIL-v0; model = GPASOIL-v0) to (data = GPASOIL-v1; model = GPASOIL-v0)). Here we focus on the effect of using new datasets to describe the different drivers with using the same model as the one of Ringeval et al.¹¹ (i.e. model = GPASOIL-v0). Figure S7 shows first that (data = GPASOIL-v1) modified the simulations of $P_{i-lab}$ and $P_{x-tot}$ for both cropland and grassland (comparison between 1st and last line of Fig. S7). In particular, $P_{i-lab}$ of grassland is strongly reduced with (data = GPASOIL-v1) (20 kgP ha$^{-1}$ at the global scale in 2005) as compared to (data = GPASOIL-v0) (208 kgP ha$^{-1}$). Sensitivity tests (Fig. S7) show that BIOG explains the most the difference between (data = GPASOIL-v1) and (data = GPASOIL-v0) for $P_{x-tot}$. Conversely, FARM is the main driver of the difference between (data = GPASOIL-v1) and (data = GPASOIL-v0) for $P_{i-lab}$ of grassland, and a combination of
<figure><img src="image_4.png" /><figcaption>Fig. 3 Comparison at country-scale between the two estimates of P in manure applied on cropland + grassland used in our approach ($fP_{manu,1}^{out→x-tot}$ in and $fP_{manu,2}^{out→x-tot}$). Temporal country-scale average of Demay et al.⁹ (red curve) is used to scale the country-scale average of $fP_{manu,2}^{out→x-tot}$ while the year-to-year variability of $fP_{manu,2}^{out→x-tot}$ follows the one of $fP_{manu,1}^{out→x-tot}$. The mean of our two estimates (subscript ‘mean’) is plotted in black. Independent estimate based on a livestock budget ($fP_{manu,indirect}$) is also plotted in orange. All fluxes are expressed in kgP (ha of cropland + grassland)⁻¹. Estimate 1 (respectively estimate 2) corresponds to the min (respectively the max) value for all grid-cells within the country considered. “World” in this figure corresponds to the sum of all countries considered in both our estimates and Demay et al.⁹ and thus excludes few countries not available in Demay et al.⁹. The number in the left-top corner of each panel corresponds to the root mean square error (in kgP (ha of cropland + grassland))⁻¹) between orange and black curves for years in common.</figcaption></figure>

BIOG and FARM explains the difference for $P_{i-lab}$ of cropland. The key roles played by BIOG and FARM in the spatial distribution of simulated soil P pools were already demonstrated in Ringeval et al.¹¹.

The biggest difference in FARM between (data = GPASOIL-v0) and (data = GPASOIL-v1) for grassland concerns the P uptake which is much larger in (data = GPASOIL-v1) than in (data = GPASOIL-v0) (Fig. S8). In (data = GPASOIL-v0), the plant uptake of grassland was constrained to be equal to ~90% of the total P input¹³.

The spatial extent of simulated $P_{i-lab}$ is larger in (data = GPASOIL-v1) (~36000 grid-cells and 1481 millions of ha for cropland; ~41000 grid-cells and 3237 millions of ha for grassland) than in (data = GPASOIL-v0) (~15000 grid-cells and 1048 millions of ha for cropland; ~20000 grid-cells and 2215 millions of ha for grassland). This is explained by a Boolean treatment of cropland in Bouwman et al.¹³ (with grid-cells with small cropland fraction not considered) and resulting to exclusion of around 1/3 of the real cropland area in GPASOIL-v0.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

**Fig. 4** Half-degree grid-cell distribution of each flux parameter involved in (model = GPASOIL-v1). Each panel corresponds to one parameter. In each panel, the red bar gives the min - max range of values for this parameter provided by Wang *et al.*¹⁸ (mainly from Fig. S2 of Wang *et al.*¹⁸), the magenta bar corresponds to the values required to make the output flux from a pool smaller than the size of this pool, the blue bar corresponds to the grid-cell distribution of this parameter in the current study and green bars correspond to the distribution of variables used to compute this parameter in the current study. The values required to make the output flux from a pool smaller than the size of this pool are given in magenta. For instance, the gross flux from $P_{i-sec}$ to $P_{i-sol}$ (called $fP_{desorp}^{i-sec \to i-sol}$) is computed with Eq. 4: $fP_{desorp}^{i-sec \to i-sol} = k^{i-sec \to i-sol} \cdot P_{i-sec}$. The inequality $fP_{desorp}^{i-sec \to i-sol} < P_{i-sec}$ at daily time-step is equivalent to $k^{i-sec \to i-sol} < 1$ (with $k^{i-sec \to i-sol}$ in day⁻¹) and thus we compared the distribution of $k^{i-sec \to i-sol}$ to 1 in the panel (d). Similarly, the gross flux from $P_{i-sol}$ to $P_{i-sec}$ (called $fP_{sorp}^{i-sol \to i-sec}$) is computed as follows: $fP_{sorp}^{i-sol \to i-sec} = k^{i-sol \to i-sec} \cdot (P_{i-sol} / W_{abs})^b$. Thus, the inequality $fP_{sorp}^{i-sol \to i-sec} < P_{i-sol}$ at daily time-step is equivalent to $k^{i-sol \to i-sec} < \frac{W_{abs}(P_{C,\infty})^{1-b}}{W_{abs}(P_{C,\infty})}$ with $k^{i-sol \to i-sec}$ expressed in mgP (kg soil)⁻¹ day⁻¹ (mg P/L)⁻¹. Thus, we compared the distribution of $k^{i-sol \to i-sec}$ to $\frac{W_{abs}(P_{C,\infty})^{1-b}}{W_{abs}(P_{C,\infty})}$ in the panel (c) Green bars correspond to the distribution of variables used to compute each parameter. For instance, panel (a) focuses on $k^{i-lab \to i-sol}$. The equation to compute this parameter is: $k^{i-lab \to i-sol} = -4.82 + 209f_{i-sol} + 14.64f_{x-occ} + 9.26f_{i-sec} - 0.008C - 0.0003P_{i-tot\prim,\infty} - 0.018s_p$ (cf. Table 10), thus the blue bar corresponds to the distribution of $k^{i-lab \to i-sol}$ and green bars correspond to the distribution of $-4.82$ (1st green bar), $+209f_{i-sol}$ (2nd green bar), $+14.64f_{x-occ}$ (3rd green bar), etc. This shows the contribution of each variable to the value of the parameter. The parameter b is without unit, $k^{i-sec \to x-occ}$, $k^{i-occ \to i-sec}$, $k^{i-lab \to i-sol}$ and $k^{i-sec \to i-sol}$ are in day⁻¹, $k^{i-sol \to i-sec}$ and $k^{i-sol \to i-lab}$ are in mgP (kg soil)⁻¹ day⁻¹ (mg P/L)⁻¹. $k^{i-occ \to i-sec}$, $k^{i-sol \to i-sec}$ and $k^{i-sec \to i-sol}$ are log transform to express them as a sum of another variables.
<figure><img src="image_12.png" /></figure>

**Fig. 5** Effect of using (model = GPASOIL-v1) instead of (model = GPASOIL-v0) on simulated soil P pools for cropland. First line shows the simulation with (model = GPASOIL-v0) while the second line shows (model = GPASOIL-v1). All simulations have been performed with (data = GPASOIL-v1). The last line shows the difference (model = GPASOIL-v1) - (model = GPASOIL-v0). All plots are in kgP ha⁻¹. The effect of changing the soil P dynamic model is provided for different soil P pools (Pᵢ-lab, Pᵢ-sec, Pₓ-oc) and on the variable fᵢ-upns, which corresponds to the P uptake prescribed by the data but that the soil P pools simulated by the model are not able to satisfy. Soil P pools plotted correspond to the year 2018 while fᵢ-upns corresponds to the annual average over 1900–2018.

<figure><img src="image_13.png" /></figure>

**Fig. 6** Same as Fig. 5 but applied to grassland.

Effect of the update of soil P dynamic model (from (data = GPASOIL-v1; model = GPASOIL-v0) to (data = GPASOIL-v1; model = GPASOIL-v1)). Here, we prescribed (data = GPASOIL-v0) for the whole section and assessed the effect of using (model = GPASOIL-v1) instead of (model = GPASOIL-v0) (Figs. 5, 6 for cropland and grassland, respectively). Overall, for cropland, compared to (model = GPASOIL-v0), (model = GPASOIL-v1) tends to increase vs. decrease Pᵢ-lab in grid-cells where values were very low vs. very high, respectively, making the spatial soil P distribution more homogeneous (Fig. S9). For grassland, (model = GPASOIL-v1) tends to increase Pᵢ-lab almost everywhere as compared to (model = GPASOIL-v0) (38 vs 23 kgP ha⁻¹, respectively) (Fig. 6). Note that, while Pᵢ-lab is not defined equally in (model = GPASOIL-v0) and (model = GPASOIL-v1) (as Pᵢ-sol is only distinguished from Pᵢ-lab in (model = GPASOIL-v1)), the comparison
<figure><img src="image_14.png" /><figcaption>**Fig. 7** Effect of increased mineralization rates on soil P pools simulated for grassland at the global scale. In panel (a), $k_{m1} = 2.7e^{-5} \text{day}^{-1}$ and $k_{m2} = 2.7e^{-4} \text{day}^{-1}$ for respectively $P_{o-sta}$ and $P_{o-lab}$ (i.e. residence time of 100 yr and 10 yr, respectively) following (model = GPASOIL-v0). In panel (b), $k_{m1} = 1.8e^{-4} \text{day}^{-1}$ and $k_{m2} = 1.4e^{-3} \text{day}^{-1}$ (i.e. a residence time of 15 yr and 2 yr, respectively). In panel (c), same mineralization rates as in panel (b) are used but in addition, losses by erosion are set to zero.</figcaption></figure>

between the two models can nevertheless focus on $P_{i-lab}$ as the maximum value for the ratio $P_{i-sol}/P_{i-lab}$ reaches only 0.2% in (model = GPASOIL-v1). We found that (model = GPASOIL-v1) allows soil P pools simulated to be more consistent with the P plant uptake prescribed (last column of Figs. 5, 6), especially for grassland (Fig. 6).

Soil exchanges with $P_{x-occ}$ play a key-role in the difference between (model = GPASOIL-v0) and (model = GPASOIL-v1). This in particular true for grassland for which, the model compensates negative soil P budget by P taken from $P_{x-occ}$ (3rd column of Fig. 6) while de-occlusion is not allowed in (model = GPASOIL-v0). The amount of P transferred from $P_{x-occ}$ to other pools ($P_{i-sec}$ then $P_{i-lab}$) in (model = GPASOIL-v1) depends on the mineralization rates of organic pools ($P_{o-lab}$ and $P_{o-sta}$) used in the model. Sensitivity tests showed that the $P_{x-occ}$ for global grassland decreased much less in time if we increased the mineralization rates (i.e. decrease the organic residence times) (Fig. 7). Our soil P dynamic model for organic P is too simplistic (no biochemical mineralization, no link with carbon, no microbial pools<sup>27</sup>) and building a more realistic model based on measurements would require a specific study. As it is likely that the temporal variations of global grassland soil P were small during the 1st half of the 20th century, we chose mineralization rates that allow to approach this global pattern in our simulations, i.e. $k^{o-sta \to i-lab} = 1.8e^{-4} \text{day}^{-1}$ and $k^{o-lab \to i-lab} = 1.4e^{-3} \text{day}^{-1}$ (i.e. a residence time of 15 yr and 2 yr, respectively) (Fig. 7). Such configuration is called GPASOIL-v1.1. Even with increased mineralization, a decrease of $P_{x-occ}$ over time remained, partly attributed to losses through erosion and to transfers to sustain the prescribed P uptake (Fig. 7c). The effect on the global average of the change in mineralization rate is large (from 149 kgP ha<sup>-1</sup> for cropland and 38 for grassland with GPASOIL-v1.0 to 187 for cropland and 91 for grassland with GPASOIL-v1.1).

We also found that the global averaged soil P pools were sensitive to the value of $P_{c,\infty}$ chosen. The choice of $P_{c,\infty}$ has a small effect on the global average for common values of $P_{c,\infty}$ (global averages of cropland $P_{i-lab}$ of 149 kgP ha<sup>-1</sup> for default value of $P_{c,\infty}$ (0.1 mgP L<sup>-1</sup>) vs 174 and 148 for kgP ha<sup>-1</sup> for $P_{c,\infty}$ equal to the 1st and 3rd quartile of $P_{c,\infty}$ found in Helfenstein *et al.*<sup>25</sup>, respectively; 2nd and 3rd line of Fig. S10) but a large effect for lowest values of $P_{c,\infty}$ such as values that could be found in Andosols or Ferralsols (global averaged $P_{i-lab}$ of 239 kgP ha<sup>-1</sup> for $P_{c,\infty} = 0.005 \text{ mgP L}^{-1}$, 1st line of Fig. S10). This suggest that more work is needed to constrain further spatial variation of $P_{c,\infty}$ (see the “Usages Notes” section).

**Best simulations, uncertainty associated and comparison to previous estimates.** In the following, we focused on the GPASOIL-v1.1 simulations as we considered them as our best estimates. Figure 8 shows the spatial distribution of the soil P budget (soil P input - output) for $P_{i-lab}$. Figure 9 shows the input/output for $P_{i-lab}$ (including soil P input/output and fluxes resulting from soil P dynamics) and soil P pools simulated for GPASOIL-v1.1 for both cropland and grassland for few countries and at the global scale. It underlines the transfer from $P_{i-lab}$ to $P_{i-sec}$ (through $P_{i-sol}$) for countries as France (in 1970s) or China (in the current time period) in cropland, when soil P input are larger than soil P output. Flux from $P_{i-sec}$ to $P_{i-lab}$ was simulated for cropland in Russia or for many countries in grassland. Such flux also occurred for France in the current time-period in cropland but with a smaller magnitude (average of 0.9 kgP ha<sup>-1</sup> yr<sup>-1</sup> over 2010–2018 in cropland for France), which corresponds to transfer from P inherited from past fertilizer applications to current plant uptake.

Figure 10 shows the comparison, for the top 0–0.3 m soil layer, between GPASOIL-v1.1 and other published modelling products available: GPASOIL-v0, Zhang *et al.*<sup>12</sup>, and McDowell *et al.*<sup>58</sup>. In parallel to the studies based on process-based model of soil P dynamics mentioned in the introduction<sup>11,12</sup>, we also included a comparison against a global dataset of soil plant available P based on a statistical model trained on a database of measured
<figure><img src="image_15.png" /><figcaption>**Fig. 8** Spatial distribution of the soil P budget (soil P input - output) for $P_{i-lab}$ for both cropland and grassland. Only plant P uptake allowed by the soil P pools simulated was considered in this budget and thus the soil P budget plotted is representative to GPASOIL-v1.1.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 9** Temporal evolution of soil $P_{i-lab}$ input/output and simulated soil P pools for several (group of) countries and for both cropland and grassland. Simulations used correspond to GPASOIL-v1.1. In left columns, only flux corresponding to soil P dynamics (weathering, mineralization, net flux from $P_{i-sec}$) are simulated while others are prescribed through (data = GPASOIL-v1).</figcaption></figure>

soil P concentration⁵⁸. McDowell et al.⁵⁸ did not focus on agricultural soils but their dataset encompassed these soils, and we filtered them in the scope of our comparison. GPASOIL-v1.1, GPASOIL-v0 and Zhang et al.¹² provide soil P pools for top 0–0.3 m in kgP ha⁻¹ for cropland (and grassland) while McDowell et al.⁵⁸ provide soil P concentration for top 0–0.2 m in mgP (kg of soil)⁻¹ for any land-use categories at 1 km² resolution. We used the land-use information used in McDowell et al.⁵⁸ to filter their estimate and computed soil P for cropland and grassland only, before regrading them to half-degree regular resolution. McDowell et al.⁵⁸ provide soil P pools for cropland, rangeland, improved grassland, forest and non-productive areas. According to our definition of grassland, we gathered ‘improved grassland’ and ‘rangeland’ classes of McDowell et al.⁵⁸ to compare to our grassland soil P pools estimates. Then, we convert soil P concentration in kgP ha⁻¹ for top 0–0.3 m by using the bulk density used in our study and by assuming that the P concentration provided by McDowell et al.⁵⁸ is also representative to 0–0.3 m.

Global average of $P_{i-lab}$ decreased for both cropland and grassland from GPASOIL-v0 to GPASOIL-v1.1. The relative decrease was even larger for grassland (from 208 to 91 kgP ha⁻¹) than for cropland (from 253 to 187 kgP ha⁻¹). $P_{i-lab}$ reached 238 kgP ha⁻¹ in Zhang et al.¹² for cropland, and 62 and 32 kgP ha⁻¹ in McDowell et al.⁵⁸ for cropland and grassland, respectively. $P_{i-lab}$ provided by McDowell et al.⁵⁸ is difficult to compare to other estimates used here as it is based on Olsen P extraction vs. pool design based on Hedley fractionation used in other studies. For cropland, we simulated larger P in South America in GPASOIL-v1.1 than in GPASOIL-v0, a lower P
<figure><img src="image_17.png" /><figcaption>**Fig. 10** Comparison between previous estimates (McDowell *et al.*⁵⁸, Zhang *et al.*, GPASOIL-v0) and our best estimates (GPASOIL-v1.1) for Pᵢ-lab and Pₓ-tot for cropland and grassland when available. Grassland were not simulated in Zhang *et al.*¹² (and while Sattari *et al.*⁸², based on the same approach, focused on grassland, they did not provide nor discuss the distribution in grasslands). McDowell *et al.*⁵⁸ provided Olsen P concentration for top 0–0.2 m at 1 km² resolution for any land-use categories. We used the land-use information used in their study to filter their estimates and computed soil P for cropland and grassland only before changing the projection and regidding to half-degree resolution. Soil P concentration was converted in kgP ha⁻¹ for top 0–0.3 m by using the bulk density used in our study and by assuming that the P concentration is representative to 0–0.3 m top soil layer. Note that Pᵢ-lab of McDowell *et al.*⁵⁸ is based on Olsen P extraction while pools design in other estimates plotted in this Figure is based on Hedley fractionation. McDowell *et al.*⁵⁸ did not provide Pₓ-tot.</figcaption></figure>

content in Africa in GPASOIL-v1.1 than in GPASOIL-v0 while some patterns remain stable in the two versions (high soil P in India and China, low P content in West Europe and Russia).

Differences in the global averaged soil P between GPASOIL-v1.1 and McDowell *et al.*⁵⁸ are likely to be caused by the shallower depth used in McDowell *et al.*⁵⁸ and the different P fractions measured (Olsen P) or estimated (Pᵢ-lab). The effect of difference in soil depth considered is not straightforward, depending on the farming practices (tillage, fertilization, etc.⁶¹). However, we do know the likely effect of the Olsen-P vs Pᵢ-lab. Labile inorganic P from Hedley is close to the Colwell P test used in Australia; both use a 16 hours soil extraction with 0.5 M HCO₃⁻. There are equations to estimate Colwell P from Olsen-P while these equations differ according to soil P sorption capacity (e.g., Colwell-P is 2.869 × Olsen-P - 2.93 for non-calcareous soils and 1.376 × Olsen-P + 8.80 for calcareous soils)⁶². After converting the values derived from McDowell *et al.*⁵⁸ for (more common) non-calcareous soils, we get 165 and 81 kgP ha⁻¹, for cropland and grassland respectively; meaning that the global values calculated in McDowell *et al.*⁵⁸ are within the same range as Pᵢ-lab simulated here.

As the biogeochemical background was found to be a major driver of the spatial distribution of Pₓ-tot¹¹, and because the same biogeochemical background was considered for both cropland and grassland, Pₓ-tot was similar between cropland and grassland for GPASOIL-v1.1. The same reasoning could be applied to GPASOIL-v0. Pₓ-tot for cropland shows also some similarities between Zhang *et al.*¹² and GPASOIL-v0 as both relied on the same dataset to approach the initial soil P content⁴. Differences in Pₓ-tot between GPASOIL-v0 and GPASOIL-v1.1 were linked to differences in Pₓ-tot of unmanaged soil used in each study (Yang *et al.*⁴ for v0 vs. He *et al.*¹⁴ for v1).

The coefficients of variation for GPASOIL-v1.1 for Pᵢ-lab and Pₓ-tot were plotted in Fig. 11 for both cropland and grassland. These coefficients were due to uncertainty in drivers estimates only. Uncertainty was larger for Pᵢ-lab than for Pₓ-tot as the uncertainty in BIOG explained exclusively the uncertainty in Pₓ-tot while all drivers contributed to the uncertainty in Pᵢ-lab. Uncertainty of Pₓ-tot was similar for both cropland and grassland (global averaged coefficient of variation of 0.17) while much larger for grassland (0.54) than for cropland (0.22) for Pᵢ-lab. This was explained by the large uncertainty in plant uptake due to uncertainty in both grassland NPP and P concentration of grass. Using mean values for variables related to FARM (instead of random ones) makes the global averaged coefficients of variation of Pᵢ-lab to decrease up to 0.14 for cropland and 0.24 for grassland, underlying the large contribution of FARM uncertainty in uncertainty of simulated Pᵢ-lab.

**Evaluation of the simulated soil P pools at different scales (country, continental/watershed, global scale).** The evaluation of predicted soil P is challenging due to lack of data compiling Hedley P measurements from cropland/grassland soils at the global scale. Different chemical extraction protocols are used worldwide with difficulties to find correspondence among different extractions. In addition, some chemical extraction in cropland/grassland soils at long-term field sites can be found in the literature but these sites usually encompass different P treatments (varying within their amount/rate/nature of fertilizer) with difficulties to link with farming practices occurring in the region grid-cell including the sites. Here we relied on 3 regional datasets based on soil monitoring networks that carried out measurements on soil samples: RMQS over France⁶³, LUCAS
<figure><img src="image_18.png" /><figcaption>**Fig. 11** Coefficients of variation (CV) for P<sub>1-lab</sub> and P<sub>x-tot</sub> for both cropland and grassland. CV was computed by using mean and standard-deviation among 100 simulations performed to assess how the uncertainty in the driver estimates has an effect on the soil P pools simulated. In the 100 simulations used for this plot, the uncertainty of all drivers was considered.</figcaption></figure>

<table><caption>**Table 12.** Information about the regional dataset used to evaluate the soil P distribution simulated in our study (GPASOIL-v1.1).</caption><thead><tr><th>Name</th><th>RMQS</th><th>LUCAS</th><th>STS</th></tr></thead><tbody><tr><td>Reference</td><td>Jolivet et al.<sup>63</sup></td><td>Jones et al.<sup>19</sup></td><td>TFI<sup>64</sup></td></tr><tr><td>Website</td><td>Data can be downloaded from: https://doi.org/10.15454/QSXKGA</td><td>Data can be downloaded after a request at ESDA website<sup>60</sup>; https://esdac.jrc.ec.europa.eu/content/lucas2015-topsoil-data</td><td>http://soiltest.tfi.org.</td></tr><tr><td>Spatial extent</td><td>France</td><td>Europe</td><td>Canada and USA</td></tr><tr><td>Years of sampling</td><td>1<sup>st</sup> compaign of measurement (2002–2009)</td><td>2015</td><td>Different years were available in STS but we used the last year in common to GPASOIL-v1.1 and STS (i.e. 2015).</td></tr><tr><td>Soil layer</td><td>Value of the 1<sup>st</sup> soil layer considered was used (this soil layer is not necessary exactly equal to 0–30 cm but very close to it)</td><td>0–0.2 m</td><td>Value for the top 0–0.3 m was computed by using the P concentration available for each soil layer within this horizon in STS.</td></tr><tr><td>Chemical extraction for P</td><td>Olsen P</td><td>Olsen P (information given in Töth et al.<sup>81</sup>)</td><td>All measures were done using different P extraction method and then all data were converted to Bray and Kurtz P1 Equivalent.</td></tr><tr><td>Land-cover used to approach cropland and grassland</td><td>The categories “cropland” and “fruits” for the land-cover “occupation 1” were used for our land-cover cropland. The category “grassland” was used for our land-cover “grassland”.</td><td>Categories “cropland” and “grassland” for the “main group of Land Cover” (LC0) were considered</td><td>Values are provided at county/province without distinction between cropland and grassland</td></tr><tr><td>Unit</td><td>gP<sub>2</sub>O<sub>5</sub> (kg of soil)<sup>-1</sup></td><td>mgP (kg of soil)<sup>-1</sup></td><td>mgP (kg of soil)<sup>-1</sup></td></tr><tr><td>Treatments performed for comparison to our dataset</td><td>- change in projection<br/>- average of all sites within 0.5 grid-cell (with average done separately for cropland and grassland)<br/>- conversion in kgP ha<sup>-1</sup> by using the 0.5 resolution bulk density and volumetric fraction of coarse fragments used in our approach</td><td>- average of all sites within 0.5 grid-cell (with average done separately for cropland and grassland)<br/>- conversion in kgP ha<sup>-1</sup> by using the 0.5 resolution bulk density and volumetric fraction of coarse fragments used in our approach</td><td>- depth average to compute 0–0.3 m value<br/>- conversion in kgP ha<sup>-1</sup> by using the country-scale value of bulk density that we computed from the 0.5 resolution bulk density and volumetric fraction of coarse fragments used in our approach</td></tr></tbody></table>

over Europe<sup>19</sup>, and STS over USA/Canada<sup>64</sup>. They differ in the chemical extraction used, the spatial extent and resolution, and the distinction (or not) between cropland and grassland. Characteristics of each dataset is given in Table 12.

For both RMQS and LUCAS, we averaged all sites within each 0.5 grid-cell separately for cropland and grassland, by considering land-uses types of each dataset (Table 12). For both RMQS and LUCAS, a detection limit for P-Olsen measurement was given in the dataset description (5 mgP (kg of soil)<sup>-1</sup> for RMQS, 10 mgP (kg of soil)<sup>-1</sup> for LUCAS) but we did not exclude the sites below this threshold to prevent bias in our grid-cell averages. All dataset (RMQS, LUCAS, STS) provided soil P in terms of concentration and we converted them in kgP ha<sup>-1</sup> by using the half-degree grid-cell bulk density and volumetric fraction of coarse fragments used previously in our study and by considering the P concentration given in each dataset representative to the top 0–0.3 m. The soil P content at half-degree resolution in kgP ha<sup>-1</sup> from the 3 datasets were provided in Figs. S11–S13.

As the datasets and our simulation did not rely on the same chemical extraction (pools representative to Hedley extraction in our simulation vs Olsen P for RMQS and LUCAS vs Bray and Kurtz P1 Equivalent for STS, Table 12), we focused our evaluation on the comparison of the relative spatial distribution. To do so, we first mapped our simulation and each dataset by using deciles (Figs. 12–14). As the number of grid-cells considered vary between our GPASOIL-v0 and GPASOIL-v1.1, we masked the grid-cells in the dataset that are excluded from our simulations (Figs. 12–14). Spearman correlation was computed on the absolute value (i.e. not in decile but in kgP ha<sup>-1</sup>) and provided in Table 13. In the scope to compare our simulation to STS, we computed first state/province scale values from our simulation before computing relative distribution with decile. We excluded from the comparison the states/provinces for which our simulation does not provide 75% of the land in farm for the considered state/province. The treatment of our simulation to compare it to STS was plotted in Fig. S14. The year used in our simulation for the comparison to the different datasets was the exact year of the dataset for GPASOIL-v1.1 (2002–2009 average for RMQS, 2015 for LUCAS, 2015 for STS) and 2005 for GPASOIL-v0 (as it is the only one year available in this latter). The comparison distinguished cropland and grassland when possible (RMQS and LUCAS).
Pᵢ-lab category, cropland,
McDowell et al. (2023)

LUCAS category, cropland,
same mask as McDowell (2023)

Pᵢ-lab category, grassland,
McDowell et al. (2023)

LUCAS category, grassland,
same mask as McDowell (2023)

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

Pᵢ-lab category, cropland,
Zhang et al. (2017)

RMQS category, cropland,
same mask as Zhang et al. (2017)

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /></figure>

Pᵢ-lab category, grassland,
GPASOIL-v0

RMQS category, grassland,
same mask as GPASOIL-v0

Pᵢ-lab category, cropland,
GPASOIL-v0

RMQS category, cropland,
same mask as GPASOIL-v0

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

Pᵢ-lab category, cropland,
GPASOIL-v1.1

RMQS category, cropland,
same mask as GPASOIL-v1.1

Pᵢ-lab category, grassland,
GPASOIL-v1.1

RMQS category, grassland,
same mask as GPASOIL-v1.1

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

RMQS category, cropland,
initial mask

RMQS category, grassland,
initial mask

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

**Fig. 12** Comparison between RMQS and either McDowell *et al.*<sup>58</sup>, Zhang *et al.*<sup>12</sup>, GPASOIL-v0 or GPASOIL-v1.1. RMQS compiled Olsen P with a cropland vs grassland distinction over France. The comparison focuses on the distribution in deciles as the simulations and RMQS were not based on same chemical extraction. All RMQS sites within half-degree grid-cell were averaged, with a distinction between cropland and grassland. Years used for the modelling approach to perform the comparison are given in Table 13. As the number of grid-cells vary between the modelling approach of McDowell *et al.*<sup>58</sup>, Zhang *et al.*<sup>12</sup>, GPASOIL-v0 or GPASOIL-v1.1, we masked the RMQS dataset according to the mask of each modelling approach output. Note that both RMQS and McDowell *et al.*<sup>58</sup> are based on Olsen P.

Overall, we found that GPASOIL-v1.1 had higher correlation than GPASOIL-v0 for grassland over France when using RMQS (+0.32 for GPASOIL-v1.1 vs non-significant correlation for GPASOIL-v0), for cropland over Europe (+0.18 for GPASOIL-v1.1 vs non-significant correlation for GPASOIL-v0, both by using LUCAS) and for Canada/USA (no distinction between cropland and grassland) (+0.32 for GPASOIL-v1.1 vs non-significant correlation for GPASOIL-v0). In particular, GPASOIL-v1.1 was able to catch the higher soil P content in grassland of West of France (Fig. 12). However, a similar agreement was not found between GPASOIL-v1.1 and LUCAS after extracting French sites from LUCAS (Table 13 and Fig. S15). This can be explained by a relatively low spearman correlation between RMQS and LUCAS over France for grassland (0.47; p-value = 0.00). GPASOIL-v1.1 was also able to catch the high soil P content in states/provinces of the north centre and south centre of the USA and Canada.

The update from GPASOIL-v0 to GPASOIL-v1.1 did not improve the correlation with datasets neither for cropland in France nor for grassland in Europe (Table 13). In particular, GPASOIL-v1.1 simulated a low soil
Pᵢ-lab category, cropland,
McDowell et al. (2023)

<figure><img src="image_36.png" /></figure>

Pᵢ-lab category, cropland,
Zhang et al. (2017)

<figure><img src="image_37.png" /></figure>

Pᵢ-lab category, cropland,
GPASOIL-v0

<figure><img src="image_38.png" /></figure>

Pᵢ-lab category, cropland,
GPASOIL-v1.1

<figure><img src="image_39.png" /></figure>

LUCAS category, cropland,
same mask as McDowell (2023)

<figure><img src="image_40.png" /></figure>

LUCAS category, cropland,
same mask as Zhang et al. (2017)

<figure><img src="image_41.png" /></figure>

LUCAS category, cropland,
same mask as GPASOIL-v0

<figure><img src="image_42.png" /></figure>

LUCAS category, cropland,
same mask as GPASOIL-v1.1

<figure><img src="image_43.png" /></figure>

LUCAS category, cropland,
initial mask

<figure><img src="image_44.png" /></figure>

Pᵢ-lab category, grassland,
McDowell et al. (2023)

<figure><img src="image_45.png" /></figure>

LUCAS category, grassland,
same mask as McDowell (2023)

<figure><img src="image_46.png" /></figure>

Pᵢ-lab category, grassland,
GPASOIL-v0

<figure><img src="image_47.png" /></figure>

Pᵢ-lab category, grassland,
GPASOIL-v1.1

<figure><img src="image_48.png" /></figure>

LUCAS category, grassland,
same mask as GPASOIL-v0

<figure><img src="image_49.png" /></figure>

LUCAS category, grassland,
same mask as GPASOIL-v1.1

<figure><img src="image_50.png" /></figure>

LUCAS category, grassland,
initial mask

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

**Fig. 13** Comparison between LUCAS and either McDowell *et al.*⁵⁸, Zhang *et al.*¹², GPASOIL-v0 or GPASOIL-v1.1. LUCAS compiled Olsen P with a cropland vs grassland distinction over Europe. The comparison focuses on the distribution in deciles as the simulations and LUCAS were not based on same chemical extraction. All LUCAS sites within half-degree grid-cell were averaged, with a distinction between cropland and grassland. Years used for the modelling approach to perform the comparison are given in Table 13. As the number of grid-cells vary between the modelling approach of McDowell *et al.*⁵⁸, Zhang *et al.*¹², GPASOIL-v0 or GPASOIL-v1.1, we masked LUCAS according to the mask of each modelling approach output. Note that LUCAS was part of the data source used in McDowell *et al.*⁵⁸ to train their statistical model.

P content in cropland of north of France contrary to what was found in RMQS. A low soil P content was also found in the north of France in GPASOIL-v0 (Fig. 12). In GPASOIL-v1.1, this pattern was explained by a large P uptake of root crops (sugar beet, potatoes) in that region. Overall, the soil P pools simulated in north of Europe are sensitive to root crop parameters related to uptake and residues. GPASOIL-v1.1 has some difficulties to catch the increasing soil P along a south-east to north-west gradient found in both cropland and grassland at the European scale in LUCAS (Fig. 13). In particular, the spearman correlation between LUCAS and GPASOIL-v1.1 was significantly negative for grassland (−0.13, p-value = 0.00, Table 13). GPASOIL-v1.1 had also difficulties to reproduce low soil P in East and West of Canada (Fig. 14).

Correlations between McDowell *et al.*⁵⁸ and regional soil P datasets are higher than correlations found with GPASOIL-v1.1 for both RMQS and LUCAS, but not for STS. The high correlations between McDowell *et al.*⁵⁸ and LUCAS is explained by the fact that LUCAS was one the data source used to train the statistical model in McDowell *et al.*⁵⁸. Also, the fact that both RMQS and McDowell *et al.*⁵⁸ are based on Olsen P can contribute to explain the high correlation between them. Difficulties of McDowell *et al.*⁵⁸ to match STS (that did not provide
<figure><img src="image_53.png" /><figcaption>**Fig. 14** Comparison between STS and either McDowell *et al.*<sup>58</sup>, GPASOIL-v0 or GPASOIL-v1.1. State/province averages were performed on simulations to allow the comparison to STS. The comparison focuses on the distribution in deciles as the simulations and STS were not based on the same chemical extraction. Note that McDowell *et al.*<sup>58</sup>, GPASOIL-v0, GPASOIL-v1.1 are not representative to the same year (Table 13). We excluded from the comparison the states/provinces for which our simulation does not provide 75% of the land in farm for the considered states/province.</figcaption></figure>

<table><caption>**Table 13.** Spearman correlations computed between the different dataset (RMQS, LUCAS, STS) and P<sub>i-lab</sub> simulated by McDowell *et al.*<sup>58</sup>, Zhang *et al.*<sup>12</sup>, GPASOIL-v0 or GPASOIL-v1.1. Bold front is used for p-value ≤ 0.05.</caption><thead><tr><th>Dataset used for the comparison</th><th>Spatial extent and years of sampling</th><th>Modeling approach</th><th>Years used in the modeling approach to compare with the dataset</th><th>Spearman correlation over cropland</th><th>Spearman correlation over grassland</th><th>Spearman correlation over cropland + grassland</th></tr></thead><tbody><tr><td rowspan="4">RMQS</td><td rowspan="4">France 2002–2009</td><td>McDowell et al.<sup>58</sup></td><td>2009</td><td>+0.33 (p-value = 0.00)</td><td>+0.36 (p-value = 0.00)</td><td></td></tr><tr><td>Zhang et al.<sup>12</sup></td><td>2010</td><td>-0.09 (p-value = 0.27)</td><td></td><td></td></tr><tr><td>GPASOIL-v0</td><td>2005</td><td>+0.03 (p-value = 0.77)</td><td>-0.11 (p-value = 0.27)</td><td></td></tr><tr><td>GPASOIL-v1.1</td><td>2002–2009</td><td>+0.01 (p-value = 0.94)</td><td>+0.32 (p-value = 0.00)</td><td></td></tr><tr><td rowspan="4">LUCAS</td><td rowspan="4">France 2015</td><td>McDowell et al.<sup>58</sup></td><td>2009</td><td>+0.39 (p-value = 0.00)</td><td>+0.56 (p-value = 0.00)</td><td></td></tr><tr><td>Zhang et al.<sup>12</sup></td><td>2010</td><td>+0.01 (p-value = 0.92)</td><td></td><td></td></tr><tr><td>GPASOIL-v0</td><td>2005</td><td>+0.04 (p-value = 0.68)</td><td>-0.22 (p-value = 0.02)</td><td></td></tr><tr><td>GPASOIL-v1.1</td><td>2015</td><td>-0.01 (p-value = 0.89)</td><td>+0.07 (p-value = 0.34)</td><td></td></tr><tr><td rowspan="4">LUCAS</td><td rowspan="4">Europe 2015</td><td>McDowell et al.<sup>58</sup></td><td>2009</td><td>+0.56 (p-value = 0.00)</td><td>+0.58 (p-value = 0.00)</td><td></td></tr><tr><td>Zhang et al.<sup>12</sup></td><td>2010</td><td>-0.07 (p-value = 0.03)</td><td></td><td></td></tr><tr><td>GPASOIL-v0</td><td>2005</td><td>+0.05 (p-value = 0.11)</td><td>-0.10 (p-value = 0.01)</td><td></td></tr><tr><td>GPASOIL-v1.1</td><td>2015</td><td>+0.18 (p-value = 0.00)</td><td>-0.13 (p-value = 0.00)</td><td></td></tr><tr><td rowspan="3">STS</td><td rowspan="3">Canada and USA 2015</td><td>McDowell et al.<sup>58</sup></td><td>2009</td><td></td><td></td><td>-0.20 (p-value = 0.35)</td></tr><tr><td>GPASOIL-v0</td><td>2005</td><td></td><td></td><td>-0.23 (p-value = 0.35)</td></tr><tr><td>GPASOIL-v1.1</td><td>2015</td><td></td><td></td><td>+0.32 (p-value = 0.03)</td></tr></tbody></table>
separate values for cropland vs grassland) could be partly explained by the fact that we considered all “range-
land” pixels in McDowell et al.⁵⁸ as grassland, contributing to unrealistic contribution of rangeland soil P pools
to agricultural pools at states/provinces scale.

At that stage, it is difficult to understand if the mismatch between GPASOIL-v1.1 and the different soil P data-
set was explained by misrepresentation of some drivers (in particular FARM or SPRO), incorrect soil P dynamics
modeling or uncertainties related to the use of the different dataset for the comparison (see next Section).

Usages Notes

The soil P pools available in the GPASOIL-v1.1 dataset are model estimates and not free from error due to imperfect modeling (e.g. soil organic P dynamics), inaccurate/uncertainties in inputs (e.g. use of dataset about N for manure to derive P in manure, bulk density not representative to agricultural soils, etc.), limitation in the estimate of some drivers (e.g. difficulty to represent P uptake in grassland), misreporting in agricultural census statistics at the basis of few input (e.g. crop P uptake), and constant in space information (e.g. soil plough layer). Only few sources of uncertainty are considered in our approach (uncertainty related to the datasets describing the drivers to provide the coefficients of variation, uncertainty related to the mineralization rates for GPASOIL-v1.0 and GPASOIL-v1.1).

Caution is particularly required if the user wants to focus on a restricted geographical area (e.g. ≤ continent) of the estimates.

The variables used as input of the soil P dynamic model (the so-called drivers) provided with the soil P pools
simulated thanks to our approach should not be used separately to the soil P pools and the users should instead
use original dataset if he/she wants to focus on a given driver.

It is worth noting that GPASOIL-v1.1 does not exactly match the spatial distribution given by regional soil P dataset. This mismatch can be explained by misrepresentation of some drivers (in particular FARM or SPRO), incorrect soil P dynamics modeling or uncertainties related to the use of the different dataset for the comparison. About the misrepresentation of some drivers, we made use of the fact that both RMQS and LUCAS provide, in addition to soil P used in the evaluation of our simulations, information about pH, soil organic carbon, and soil texture, i.e. soil properties that are involved in our soil P dynamics parameterization (SPRO driver). Figure S16 shows the comparison for these soil properties between Soilgrids 2.0 (used in our approach), RMQS and LUCAS. Some spatial patterns are not consistent between Soilgrids 2.0, RMQS and LUCAS (e.g. pH). Also, the comparison lets suggest that we would need to distinguish cropland vs grassland for soil parameters involved in the soil P dynamics model.

To improve our understanding of the mismatch between GPASOIL-v1.1 and soil P dataset, it could also be
particularly interesting to investigate if the correlation between GPASOIL-v1.1 and LUCAS would be improved
in case of simulations over Europe where both soil P input/output and soil dynamics parameterization were
improved. European simulation of GPASOIL-v1.1 could be forced by soil P input/output of Panagos et al.⁶⁵. How
the correlation with LUCAS evolved in sensitivity tests where some crop-specific parameters related to uptake
and residues (in particular for root crops) are made varying should also be studied. Soil P dynamics could be
improved for European simulations as LUCAS provide information about soil oxides, which can allow to use
initial parameterizations of Wang et al.¹⁸ (instead of modified ones as descried in Table 10) and in particular
spatially varying concentration of P in soil solution at steady-state (P<sub>c,∞</sub>). Also, it could be interesting to evaluate
if the soil P dynamic model, even without being fed by oxide concentrations, can match the known global distri-
bution of soil phosphorus retention potential⁶⁶. Finally, the extension of dataset such as LUCAS or RMQS with
Hedley fractionation measurements would also allow a more straightforward comparison to modelling studies.

Code availability

The scripts of the modelling approach described in this paper was made available at recherche.data.gouv.fr⁵⁷.
Our code is available within the GPASOILScripts.tgz file along the tar file that contain the GPASOIL-v1.0 and
GPASOIL-v1.1 soil P pools.

GPASOILScripts.tgz contains the whole file tree and associated scripts that were used to generate GPASOIL-v1_output.tgz. The directory GPASOIL/ contains few sub-directories:

- GrabData_and_PrepDriver/ that contains all scripts (bash and python) used to generate the input of the soil P dynamic model. These input correspond to the different drivers described in the current manuscript. The input files generated (in netcdf format) are within the directory output_prepDriver/ and are here provided. The procedures to download the original dataset (if this dataset is available on the web) used to generate the input files and information are given for each driver in specific README.txt files.

- main/ contains the soil P dynamic model. Different version are provided: (model = v0), (model = v1.0) and (model = v1.1). The 2nd one (i.e. (model = v1.0)) was used to generate GPASOILv1.0 while the latter (i.e. (model = v1.1)) was used to generate GPASOILv1.1

- evaluation/ contains the script to compare the model output to regional databases. The databases are not provided.
General information about the procedure to get the regional databases are provided in a specific README.txt file.

- output_main/ is the directory that receives the different soil P maps simulated.

A container generated with Singularity 3 was also provided and allows users to run the different scripts on
a server or local computer without issue of incompatibility about Ubuntu distribution or Python packages.
Modeling and analysis were performed in using Python (Python Software Foundation. Python Language
Reference, version 3.6.10., available at http://www.python.org, last access: January 2020).

Received: 13 February 2023; Accepted: 14 November 2023;

Published online: 02 January 2024
## References

1. Hou, E. et al. Global meta-analysis shows pervasive phosphorus limitation of aboveground plant production in natural terrestrial ecosystems. *Nat. Commun.* **11**, (2020).

2. Kvakić, M. et al. Quantifying the Limitation to World Cereal Production Due To Soil Phosphorus Status. *Glob. Biogeochem. Cycles* https://doi.org/10.1002/2017GB005754 (2018).

3. Carpenter, S. R. Eutrophication of aquatic ecosystems: Bistability and soil phosphorus. *Proc. Natl. Acad. Sci.* **102**, 10002–10005 (2005).

4. Yang, X., Post, W. M., Thornton, P. E. & Jain, A. The distribution of soil phosphorus for global biogeochemical modeling. *Biogeosciences* **10**, 2525–2537 (2013).

5. Achat, D. L., Pousse, N., Nicolas, M., Brédoire, F. & Augusto, L. Soil properties controlling inorganic phosphorus availability: general results from a national forest network and a global compilation of the literature. *Biogeochemistry* **127**, 255–272 (2016).

6. He, X. et al. Global patterns and drivers of soil total phosphorus concentration. *Earth Syst. Sci. Data* **13**, 5831–5846 (2021).

7. Elser, J. & Bennett, E. A broken biogeochemical cycle. *Nature* **478**, 29–31 (2011).

8. MacDonald, G. K., Bennett, E. M., Potter, P. A. & Ramankutty, N. Agronomic phosphorus imbalances across the world's croplands. *Proc. Natl. Acad. Sci. USA* **108**, 3086–3091 (2011).

9. Demay, J., Ringeval, B., Pellerin, S. & Nesme, T. Half of global agricultural soil phosphorus fertility derived from anthropogenic sources. *Nat. Geosci.* **16**, 69–74 (2023).

10. Sattari, S. Z., Bouwman, A. F., Giller, K. E. & van Ittersum, M. K. Residual soil phosphorus as the missing piece in the global phosphorus crisis puzzle. *Proc. Natl. Acad. Sci. USA* **109**, 6348–6353 (2012).

11. Ringeval, B. et al. Phosphorus in agricultural soils: drivers of its distribution at the global scale. *Glob. Change Biol.* **23**, 3418–3432 (2017).

12. Zhang, J. et al. Spatiotemporal dynamics of soil phosphorus and crop uptake in global cropland during the 20th century. *Biogeosciences* **14**, 2055–2068 (2017).

13. Bouwman, L. et al. Exploring global changes in nitrogen and phosphorus cycles in agriculture induced by livestock production over the 1900–2050 period. *Proc. Natl. Acad. Sci. USA* **52**, https://doi.org/10.1073/pnas.1012878108 (2013).

14. He, X. et al. Global patterns and drivers of phosphorus fractions in natural soils. *Biogeosciences* **20**, 4147–4163 (2023).

15. Wang, Y. P., Law, R. M. & Pak, B. A global model of carbon, nitrogen and phosphorus cycles for the terrestrial biosphere. *Biogeosciences* **7**, 2261–2282 (2010).

16. Helfenstein, J. et al. Combining spectroscopic and isotopic techniques gives a dynamic view of phosphorus cycling in soil. *Nat. Commun.* **9**, 3226 (2018).

17. Helfenstein, J. et al. Estimates of mean residence times of phosphorus in commonly considered inorganic soil phosphorus pools. *Biogeosciences* **17**, 441–454 (2020).

18. Wang, Y. et al. Toward a Global Model for Soil Inorganic Phosphorus Dynamics: Dependence of Exchange Kinetics and Soil Bioavailability on Soil Physicochemical Properties. *Glob. Biogeochem. Cycles* **36**, (2022).

19. Jones, A., Fernandez-Ugalde, O. & Scarpa, S. LUCAS 2015 topsoil survey: presentation of dataset and results. vol. https://doi.org/10.2760/616084, JRC121325 (EUR 30332 EN, Publications Office of the European Union: Luxembourg, 2020).

20. Chini, L. et al. Land-use harmonization datasets for annual global carbon budgets. *Earth Syst. Sci. Data* **13**, 4175–4189 (2021).

21. Poggio, L. et al. SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. *SOIL* **7**, 217–240 (2021).

22. Bronick, C. J. & Lal, R. Soil structure and management: a review. *Geoderma* **124**, 3–22 (2005).

23. Klein Goldewijk, K., Beusen, A., Doelman, J. & Stehfest, E. Anthropogenic land use estimates for the Holocene – HYDE 3.2. *Earth Syst. Sci. Data* **9**, 927–953 (2017).

24. Krinner, G. et al. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. *Glob. Biogeochem. Cycles* **19**, n/a–n/a (2005).

25. Helfenstein, J., Jegminnat, J., McLaren, T. L. & Frossard Soil solution phosphorus turnover: derivation, interpretation, and insights from a global compilation of isotope exchange kinetic studies. *Biogeosciences* **15**, 105–114 (2018).

26. Buendía, C., Kleidon, A. & Porporato, A. The role of tectonic uplift, climate, and vegetation in the long-term terrestrial phosphorous cycle. *Biogeosciences* **7**, 2025–2038 (2010).

27. McConnell, C. A., Kaye, J. P. & Kemanian, A. R. Reviews and syntheses: Ironing out wrinkles in the soil phosphorus cycling paradigm. *Biogeosciences* **17**, 5309–5333 (2020).

28. Beck, M. A. & Sanchez, P. A. Soil phosphorus movement and budget after 13 years of fertilized cultivation in the Amazon basin. *Plant Soil* **184**, 23–31 (1996).

29. Sun, Y. et al. Global evaluation of the nutrient-enabled version of the land surface model ORCHIDEE-CNP v1.2 (r5986). *Geosci. Model Dev.* **14**, 1987–2010 (2021).

30. Zhang, B. et al. Global manure nitrogen production and application in cropland during 1860–2014: a 5 arcmin gridded global dataset for Earth system modeling. *Earth Syst. Sci. Data* **9**, 667–678 (2017).

31. Xu, R. et al. Increased nitrogen enrichment and shifted patterns in the world's grassland: 1860–2016, 13, *Earth Syst. Sci. Data*, **11**, 175–187, https://doi.org/10.5194/essd-11-175-2019 (2019).

32. Lun, F. et al. Global and regional phosphorus budgets in agricultural systems and their implications for phosphorus-use efficiency. *Earth Syst. Sci. Data* **10**, 1–18 (2018).

33. Velthof, G. L., Bannink, A. & Oenema, O. Relationships between animal nutrition and manure quality.

34. Lu, C. & Tian, H. Global nitrogen and phosphorus fertilizer use for agriculture production in the past half century: shifted hot spots and nutrient imbalance. *Earth Syst. Sci. Data* **9**, 181–192 (2017).

35. Sun, Y. et al. Diagnosing phosphorus limitations in natural terrestrial ecosystems in carbon cycle models. *Earths Future* **5**, 730–749 (2017).

36. Wang, J. et al. Vegetation type controls root turnover in global grasslands. *Glob. Ecol. Biogeogr.* **28**, 442–455 (2019).

37. Zhang, X. & Wang, W. The decomposition of fine and coarse roots: their global patterns and controlling factors. *Sci. Rep.* **5**, 9940 (2015).

38. Smith, B. et al. Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model. *Biogeosciences* **11**, 2027–2054 (2014).

39. Monfreda, C., Ramankutty, N. & Foley, J. A. Farming the planet: 2. Geographic distribution of crop areas, yields, physiological types, and net primary production in the year 2000: GLOBAL CROP AREAS AND YIELDS IN 2000. *Glob. Biogeochem. Cycles* **22**, n/a–n/a (2008).

40. Sun, Y. et al. Field-Based Estimation of Net Primary Productivity and Its Above- and Belowground Partitioning in Global Grasslands. *J. Geophys. Res. Biogeosciences* **126**, (2021).

41. Kastner, T. et al. Land use intensification increasingly drives the spatiotemporal patterns of the global human appropriation of net primary production in the last century. *Glob. Change Biol.* gcb.15932, https://doi.org/10.1111/gcb.15932 (2021).

42. Wang, Y. et al. GOLUM-CNP v1.0: a data-driven modeling of carbon, nitrogen and phosphorus cycles in major terrestrial biomes. *Geosci. Model Dev.* **11**, 3903–3928 (2018).

43. Lun, F. et al. Influences of international agricultural trade on the global phosphorus cycle and its associated issues. *Glob. Environ. Change* **69**, 102282 (2021).

44. Smil, V. Phosphorus in the environment: natural flows and human interferences. *Annu. Rev. Energy Environ.* **25**, 53–88 (2000).
45. Zechmeister-Boltenstern, S. et al. The application of ecological stoichiometry to plant-microbial-soil organic matter transformations. *Ecol. Monogr.* **85**, 133–155 (2015).

46. Spohn, M. Increasing the organic carbon stocks in mineral soils sequesters large amounts of phosphorus. *Glob. Change Biol.* **26**, 4169–4177 (2020).

47. Bentsen, N. S., Felby, C. & Thorsen, B. J. Agricultural residue production and potentials for energy and materials services. *Prog. Energy Combust. Sci.* **40**, 59–73 (2014).

48. Ye, Y. et al. Carbon, Nitrogen and Phosphorus Accumulation and Partitioning, and C:N:P Stoichiometry in Late-Season Rice under Different Water and Nitrogen Managements. *PLoS ONE* **9**, e101776 (2014).

49. Wang, R. et al. Significant contribution of combustion-related emissions to the atmospheric phosphorus budget. *Nat. Geosci.* **8**, 48–54 (2015).

50. Mahowald, N. et al. Global distribution of atmospheric phosphorus sources, concentrations and deposition rates, and anthropogenic impacts. *Glob. Biogeochem. Cycles* **22**, n/a–n/a (2008).

51. Wang, R. et al. Global forest carbon uptake due to nitrogen and phosphorus deposition from 1850 to 2100. *Glob. Change Biol.* **23**, 4854–4872 (2017).

52. van Puijenbroek, P. J. T. M., Beusen, A. H. W. & Bouwman, A. F. Global nitrogen and phosphorus in urban waste water based on the Shared Socio-economic pathways. *J. Environ. Manage.* **231**, 446–456 (2019).

53. Donatello, S. & Cheeseman, C. R. Recycling and recovery routes for incinerated sewage sludge ash (ISSA): A review. *Waste Manag.* **33**, 2328–2340 (2013).

54. Borrelli, P. et al. An assessment of the global impact of 21st century land use change on soil erosion. *Nat. Commun.* **8**, (2017).

55. van den Hurk, B. et al. LS3MIP (v1.0) contribution to CMIP6: the Land Surface, Snow and Soilmoisture Model Intercomparison Project – aims, setup and expected outcome. *Geosci. Model Dev.* **9**, 2809–2832 (2016).

56. Valade, A. et al. Modeling sugarcane yield with a process-based model from site to continental scale: uncertainties arising from model structure and parameter values. *Geosci. Model Dev.* **7**, 1225–1245 (2014).

57. Ringeval, B. GPASOIL-v1, a global dataset on phosphorus in agricultural soils. *Recherche Data Gouv*. https://doi.org/10.57745/XZTW7Z (2023).

58. McDowell, R. W., Noble, A., Pletnyakov, P. & Haygarth, P. M. A Global Database of Soil Plant Available Phosphorus. *Sci. Data* **10**, 125 (2023).

59. Iizumi, T. & Sakai, T. The global dataset of historical yields for major crops 1981–2016. *Sci. Data* **7**, 97 (2020).

60. Senthilkumar, K., Nesme, T., Mollier, A. & Pellerin, S. Conceptual design and quantification of phosphorus flows and balances at the country scale: The case of France. *Glob. Biogeochem. Cycles* **26**, GB2008 (2012).

61. Tian, J. et al. Accumulation and distribution of phosphorus in the soil profile under fertilized grazed pasture. *Agric. Ecosyst. Environ.* **239**, 228–235 (2017).

62. Moody, P. W., Speirs, S. D., Scott, B. J. & Mason, S. D. Soil phosphorus tests I: What soil phosphorus pools and processes do they measure? *Crop Pasture Sci.* **64**, 461 (2013).

63. Jolivet, C., Boulonne, L. & Ratié, C. Manuel du réseau de Mesures de la Qualité des Sols. *Unité InfoSol - INRA Orléans - Fr.* (2006).

64. TFI. Soil Test Levels in North America, 2020 Summary Update. (2020).

65. Panagos, P. et al. Improving the phosphorus budget of European agricultural soils. *Sci. Total Environ.* **853**, 158706 (2022).

66. Batjes, N. H. Global distribution of soil phosphorus retention potential.

67. Bouwman, A. F. & Kram. *Integrated modelling of global environmental change: an overview of IMAGE 2.4.* (Netherlands Environmental Assessment Agency, 2006).

68. Xu, R. et al. Half-degree gridded manure and fertilizer nitrogen inputs in global grassland systems during 1860–2016, PANGAEA, https://doi.org/10.1594/PANGAEA.892940 (2018).

69. Demay, J. *Supplementary Information - Half of global agricultural soil phosphorus fertility derived from anthropogenic sources*. https://doi.org/10.57745/LEPJCS (2023).

70. Kastner, T. et al. Data supplement for ‘Land use intensification increasingly drives the spatiotemporal patterns of the global human appropriation of net primary production in the last century’, https://doi.org/10.5281/ZENODOV.5519104 (2022).

71. Lu, C. & Tian, H. Half-degree gridded nitrogen and phosphorus fertilizer use for global agriculture production during 1900-2013, https://doi.org/10.1594/PANGAEA.863323 (2016).

72. Zhang, B. et al. Manure nitrogen production and application in cropland and rangeland during 1860 - 2014: A 5-minute gridded global data set for Earth system modeling, https://doi.org/10.1594/PANGAEA.871980 (2017).

73. Klein Goldewijk, Dr. ir. C.G.M. Anthropogenic land-use estimates for the Holocene; HYDE 3.2. DANS. https://doi.org/10.17026/dans-25g-gez3 (Utrecht University, 2017).

74. Van Oost, K. et al. The Impact of Agricultural Soil Erosion on the Global Carbon Cycle. *Science* **318**, 626–629 (2007).

75. Hurtt, G. C. et al. Harmonization of land-use scenarios for the period 1500–2100: 600 years of global gridded annual land-use transitions, wood harvest, and resulting secondary lands. *Clim. Change* **109**, 117–161 (2011).

76. Chini, L. et al. *Land-Use Harmonization Datasets for Annual Global Carbon Budgets*. https://essd.copernicus.org/preprints/essd-2020-388/essd-2020-388.pdf, 10.5194/essd-2020-388 (2021).

77. Decharme, B., Martin, E. & Faroux, S. Reconciling soil thermal and hydrological lower boundary conditions in land surface models: LOWER BOUNDARY CONDITIONS OF SOIL IN LSM. *J. Geophys. Res. Atmospheres* **118**, 7819–7834 (2013).

78. He, X. et al. *Global patterns and drivers of phosphorus pools in natural soils version 2.0*. 24653728 Bytes https://doi.org/10.6084/M9.FIGSHARE.16988029.V2 (2023).

79. Yang, X. & Post, W. M. Phosphorus transformations as a function of pedogenesis: A synthesis of soil phosphorus data using Hedley fractionation method. *Biogeosciences* **8**, 2907–2916 (2011).

80. Panagos, P. et al. European Soil Data Centre 2.0: Soil data and knowledge in support of the EU policies. *Eur. J. Soil Sci.* **73** (2022).

81. Tóth, G., Guicharnaud, R.-A., Tóth, B. & Hermann, T. Phosphorus levels in croplands of the European Union with implications for P fertilizer use. *Eur. J. Agron.* **55**, 42–52 (2014).

82. Sattari, S. Z., Bouwman, A. F., Martinez Rodriguez, R., Beusen, A. H. W. & van Ittersum, M. K. Negative global phosphorus budgets challenge sustainable intensification of grasslands. *Nat. Commun.* **7**, 10696 (2016).

## Acknowledgements

This research was financially supported by INRAE (Institut national de recherche pour l'agriculture, l'alimentation et l'environnement), the Agroecosystems division of INRAE and by the French funding account CASDAR (Compte d'Affectation Spécial “Développement Agricole et Rural”) under grant agreement n° 20AIP1626879 (PhosphoBio project). We thank NASA as the source of the material used for the LUCC driver. For the driver CLIM, we acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF. We thank Mark Irvine for his help with the computing aspects, Tovo Rabemanantsoa for his help with the datasets, Bertrand Decharme for his help
with the CMIP6 output, Christelle Aluome for her help with the licence of the different datasets, and Francesca Degan (Arvalis) for helpful discussion. This work has been realized with the support of MESO@LR-Platform at the University of Montpellier. DSG and XH acknowledge support from the EJP SOIL ICONICA project, DSG acknowledges support from the ANR CLAND Convergence Institute.

**Author contributions**

B.R. designed the study, wrote the different scripts, and performed the first analysis of the results. J.D., X.H., E.H., S.M., K.E., R.W., F.L., P.B., R.W.M., P.P. provided help about their respective dataset, D.S.G., Y-P.W., L.A. provided help with the soil P dynamic model. All co-authors participated in the analysis of the results and the writing.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41597-023-02751-6.

**Correspondence** and requests for materials should be addressed to B.R.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_54.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2024