Overview

We coupled two models in our study, that is, the soil P budget model (Fig. 1) and the DPPS model, to reproduce the soil P budget for the 1970–2005 period and to estimate future P requirements in grasslands. The soil P budget model considers two system boundaries: the grassland system boundary and the grassland soil boundary. Phosphorus is leaving the grassland systems’ boundary via animal products (mainly meat and milk), livestock manure and runoff or erosion. On the other hand there is P import from cropland into grassland through animal feed.

Within the grassland soil boundary, we distinguish different sources of soil P inputs, including manure generated inside the grasslands and spread in the grasslands (grassland-based, internal, livestock manure; Fig. 1), manure produced by pigs and poultry outside the grassland boundaries (non-grassland-based, external, livestock manure; Fig. 1), mineral P fertilizer (only in intensive livestock production systems, see section below) and P from weathering and atmospheric deposition. Grass P uptake, runoff and erosion are regarded as P outflows from the grassland soil (Fig. 1). Manure generated within the grassland boundaries is assigned to three uses: spreading in grassland soil, spreading in croplands and other uses of manure. The net transfer of P between croplands and grasslands over the period 1970–2005 is the difference between P imported to grassland from cropland through livestock feed, and P export from grassland to cropland in the form of manure (Fig. 1 and Supplementary Fig. 1). The DPPS model is presented in detail below.

All abbreviations used in the soil P budget model are explained in Table 3, and Table 4 describes the P flows, the animal categories and the data sources. The data used in this paper as well as the DPPS model are available on: http://models.pps.wur.nl/content/DPPS-Grassland.

Table 3 Abbreviations used in the soil P budget model, description and units. Full size table

Table 4 Description of phosphorus flows, the animal categories involved and the data source. Full size table

Livestock production systems and animal groups

In their global classification, Seré and Steinfeld32 distinguished livestock production systems (LPSs) and mixed farming production systems (MPS). LPS includes systems in which more than 90% of dry matter fed to animals comes from rangelands, pastures, annual forages and purchased feed. MPS are systems in which more than 10% of the dry matter fed to animals comes from crop by-products or stubble, or more than 10% of the total value of production comes from non-livestock farming activities32. They further distinguished four subgroups: landless LPS; grassland LPS; rainfed MPS and irrigated MPS (Supplementary Fig. 1a). They collected data on buffaloes, cattle, goats, pigs, poultry and sheep population numbers and meat and milk production in each system. Later, Bouwman et al.23 modified the classification from Seré and Steinfeld32 for their global analysis of ruminant production systems. Landless LPS and MPS were merged into mixed-landless systems (or intensive systems) while grassland LPS was renamed into pastoral systems (Supplementary Fig. 1b,c).

Following Bouwman et al.23, in the present study also two production systems were distinguished, that is, grasslands in mixed and landless (referred to as intensive hereafter) and extensive pastoral LPSs. Within each system, two groups of animals were considered: ‘grassland-based livestock’ including asses, buffaloes, camels, dairy cattle, horses, mules, non-dairy cattle, sheep and goats, and ‘non-grassland-based livestock’ including pigs and poultry (Fig. 3). Owing to lack of data, it was not possible to include all animal categories for all the calculations (Table 4). Furthermore, non-grassland-based livestock categories were not included in products’ P flow calculations since in our definition they are not located within the grassland system boundaries (Fig. 1).

Data used

The P flows considered for different animal categories and their data sources are listed in Table 4. Phosphorus inflows to the grasslands comprise feed, P fertilizer, external manure P and atmospheric deposition; P outflows from grasslands boundaries include P in livestock products, exported manure to croplands, other use of manure P and P in runoff or erosion.

Long-term livestock production data (meat and milk) were obtained from the Integrated Model to Assess the Global Environment (IMAGE)33, which is based on FAOSTAT22 data; the average P content of meat and milk and various animal parts was obtained from the literature23.

The animal stocks within pastoral and intensive livestock production systems and their P excretion rates were used to calculate the total P in manure within the two systems (see section on livestock manure P).

The amount of feed for each country, year and animal category was estimated based on IMAGE33 output data and on the use of grass and feed crops and energy requirements of the different animal categories).

The total mineral P applications to grasslands in mixed and landless livestock production systems for 1970–2005 were from IMAGE at the country scale, based on FAOSTAT and the International Fertilizer Association34,35).

Weathering and atmospheric deposition as P inflows to grasslands and soil P erosion and runoff as P outflows were estimated from literature review and IMAGE output data. The sections below provide more details on the methods and data used for the various model components.

Livestock products

Livestock production data include milk, meat and livestock by-products. Livestock by-products include adipose tissue, skeleton, viscera, blood, skin, hair and digestive content. Country-scale meat production data for non-dairy cattle and sheep and goats, as well as milk production were obtained from FAOSTAT22. Meat and milk production for buffaloes, horses, asses, mules and camels was not included due to lack of data. The total amount of P in livestock products is the sum over P in milk, meat and livestock by-products.

The carcass weight is the most common way of expressing livestock meat production. Carcass weight is defined as the weight left after slaughter and removal of head, skin, genitourinary organs and offals. The ratio between the carcass weight and the live weight is called the dressing percentage (DP). Live weight consists of the following fractions: muscles; adipose tissue; skeleton; viscera; blood; skin; hair; and digestive content. Live weight partitioning allows for a more accurate P accounting due to large differences in P concentration in different fractions such as bones and blood.

To allocate the total national production to mixed-landless and pastoral systems, national, annual values of the fraction allocated to mixed-landless systems (Fraction intensive) were obtained from IMAGE33. A list of countries and regions as used in the present study based on IMAGE has been shown in Supplementary Table 1.

In addition, data on DPs were obtained from Kempster et al.36, who reported 53 and 50% DP for, respectively, non-dairy cattle and for sheep and goats. Live weight partitioning fraction (LWF) data and P content of those fractions were obtained from different sources (Supplementary Table 2).

Equation (1) was used to calculate P in meat and livestock by-products:

where y denotes year, c country, s production system, a animal category and i, fraction of meat and by-products.

Phosphorus in meat and livestock by-products is denoted by PLMe (kg P). LPD and DP refer to the livestock production data (kg carcass weight) and DP (kg carcass per kg live weight), respectively. LWF (kg meat or livestock by-products per kg live weight) represents live weight partitioning and FPC refers to P content (kg P per kg meat or livestock by-products).

The total amount of P in milk was calculated by multiplying the total amount of milk production with the P content of milk, as shown in equation (2). The P content was considered to be 8.4·10−4 kg P per kg milk30.

where y denotes year, c country and s production systems.

Phosphorus in milk is expressed by PLMi (kg P). MPD and MPC refer to the milk production data (kg milk) and milk phosphorus content (kg P per kg milk), respectively.

Livestock manure phosphorus

Total manure production within pastoral and intensive systems was computed from animal stocks and P excretion rates. We used P excretion rates per head for dairy and non-dairy cattle, buffaloes, sheep and goats, pigs, poultry, horses, asses, mules and camels based on various sources37,38,39,40,41 (Supplementary Table 3). We used constant excretion rates per head; in case of increasing production per head, fewer animals are needed to produce the same amount of meat or milk; total P excretion will thus decrease, and the excretion per unit of product decreases reflecting increasing conversion efficiency due to improved feeding and management. For each country, animal stocks and P in the manure for each animal category were spatially allocated across intensive and pastoral systems. For the time period 2005–2050, the distribution over these systems is provided by the Rio+20 study21.

The methodology and equations described in this section were applied to all animal categories.

The manure allocation comprises five steps: a first calculation of total excreted manure and its P content and a subsequent fractioning into four different flows (Fig. 3).

The total amount of manure P was calculated according to equation (3).

where y denotes year, c country, s production systems and a animal category).

‘Manure’ refers to the total P in manure excreted by livestock (kg P); LPN is livestock population (number of heads); ‘Excretion rate’ is the annual amount of P excreted for each animal category (kg P per head per year). The sum of manure for grassland-based animal categories (cattle, buffaloes, sheep, goats, camels, horses and asses) corresponds to the ‘Livestock excretion flow’ depicted in Fig. 1.

The amount of P (kg) in manure used for non-agricultural purposes (other uses) was calculated according to equation (4).

‘FrOthers’ expresses the fraction of manure used in any way so that it is effectively removed from agricultural systems, as it could be used as fuel or building material, or digested for generating energy.

From the remaining manure, the amount of P (kg) in the manure allocated to grazing was estimated by equation (5).

‘FrGrazing’ refers to the fraction of manure that is deposited in grasslands by grassland-based grazing animals. In the case of non-grassland-based species, it represents their deposition as they scavenge in grasslands.

The rest of the manure P was considered to be excreted in animal houses and stored, and available for use as organic fertilizer either in grasslands or croplands. Two fractions were calculated, ‘Grasslands’ and ‘Croplands’ as shown in equations (6) and (7), respectively.

‘FrGrass’ refers to the fraction of stored manure that is applied in grasslands as organic fertilizer.

The amount of P (kg) in the manure that is transferred from animal houses or stored manure to cropland as fertilizer was estimated by equation (7).

The sum of ‘Grazing’ and ‘Grasslands’ for grassland-based and non-grassland-based species corresponds to the internal and external manure input flows, respectively, as depicted in Fig. 1.

The values for the manure allocation fractions (other uses, grazing and application to grasslands) are listed in Supplementary Table 4.

Livestock feed

The IMAGE model33 used FAOSTAT data on crops for feed and energy requirements of animals to estimate the amount of feed for each country, year and animal category. Only the animal categories of non-dairy cattle, dairy cattle and sheep and goats were considered for estimating feed crop use. Thus, a 100% grass ration was assumed for the rest of grassland-based animal categories (buffaloes, horses, asses, mules and camels).

IMAGE provides feed data for 11 feed items: oil crops; maize; pulses; rice; root and tuber crops; temperate cereals; tropical cereals; crop residues; and other by-products such as residues from breweries, and grass and scavenging (such as road-side grazing, food waste and so on). FAOSTAT provides data on the amount of feed crop use per country, but does not specify feed use by animal category. In IMAGE the specification of feed crop use for individual animal categories is based on feed rations that are calibrated to match FAOSTAT data on total feed use at the scale of world regions. Therefore, disaggregation of world-region feed crop use by animal category to the country scale was needed. To do so, it was assumed for each animal category that the fraction of feed crops in the ration is proportional to livestock productivity. The total amount of feed crops available for a specific animal category within a world region was thus allocated on the basis of the productivity of the animal category considered in a specific country relative to the regional productivity. Country and regional productivity were calculated for each animal category as the country and regional amount of product divided by the country and regional animal numbers, respectively. Since this calculation was done for each animal category, productivity numbers were expressed in terms of kg of carcass weight per head (for non-dairy and sheep and goats) or kg milk per head (for dairy cattle). Then, for each country the ratio between its productivity over the regional productivity was computed, which yielded a weighing factor (first term of equation (8)). With the information on the regional feed amount and the animal numbers, the regional average feed intake was calculated for each feed item and animal category. The total amount of P for each feed item for a specific animal category within a certain country was thus the product of the weighting factor times the country’s animal numbers times the regional average feed intake times the P content of the feed item. By adding up all the feed items and animal categories, the country total P amount in livestock feed crops was calculated.

where y denotes year, c country, s production systems, a animal category, f feed item and R regional data).

PFE (kg P) is total amount of P in livestock feed. LPR refers to the livestock productivity (kg per head). It is calculated as the total amount of products associated with the animal category (carcass weight or milk) over the total number of animals. LPN refers to the animal numbers and the total amount of feed item used in a certain region as presented by FEED (kg). PFI (kg P per kg feed item) is the P content for a given feed item.

It is clear that part of the P in feed used in grassland-based systems is imported from other countries. Owing to lack of data we ignore feed trade, but in regions such as Western Europe, P in feed produced in croplands in other world regions may be a significant part of total feed use. However, this does not affect the global budget of the grassland-based systems, but may lead to overestimation of the net transfer from cropland to grassland within feed importing countries.

Livestock grass intake and grass P uptake

The estimation of grass uptake (and grass intake by animals) is based on a mass balance approach. At the animal level, P inputs (feed plus grass intake) equal the P outputs (manure plus products). Phosphorus in the livestock grass intake was assumed to be equal to grass P uptake from the soil. Thus, this approach ignores any P losses during mowing, transporting or stall-feeding of grass. The grass P uptake was calculated according to equation (9). The amount of feed (PFE) and products (PLMe and PLMi) were assumed to be zero for buffaloes, horses, asses, mules and camels. For these animal categories the amount of P in grass intake equals the excretion of P.

where y denotes year, c country and s production systems.

PGU (kg P) is the total grass phosphorus uptake and PFE (kg P) refers to the total amount of phosphorus (kg) in livestock feed.

Mineral phosphorus fertilization

The IMAGE data include the total use of mineral phosphorus fertilizers in grasslands for 1950, 1970, 1980, 1990, 1995, 2000 and 2005 at country level. The values are based on FAOSTAT22 and the International Fertilizer Industry Association34,35 for the period 1970–2000, and extrapolations for earlier and later years4,21. Those numbers were exclusively allocated to the mixed and landless production systems, as we assumed no use of mineral fertilizer in pastoral systems.

Soil phosphorus erosion

Soil P loss estimates were based on a recent modelling approach, which distinguishes two nutrient loss pathways to estimate P runoff and erosion42, that is, losses from recent nutrient applications (P rec ) in the form of fertilizer, manure or organic matter43, and a residual or ‘memory’ (P mem ) effect related to long-term historical changes in soil nutrient stocks for the top 30 cm (refs 44, 45). The approach presented by Cerdan et al.46 based on a large database of measurements was used as a basis for calculating P mem based on slope, soil texture and land cover type. This approach yields an average soil loss of 40 tonnes of soil per km2 year for Europe46.

P rec was calculated from P inputs and depends on slope (using the approach of Bogena et al.47) and was further modified by land use and soil texture, that is, those factors that reduce surface runoff according to Velthof et al.48,49

The initial P stock in the top 30 cm was taken from Yang et al.50 for the year 1900, and residual soil P was calculated from the soil P budget up till the year 2005 (ref. 42). All inputs and outputs of the soil balance were assumed to occur in the top 30 cm; the soil budget model replaces P enriched or depleted soil material lost at the surface by erosion with fresh soil material (with the initial soil P content) at the bottom.

Simulated cumulative P loss amounted to 33 Tg P over the 1970–2005 period, which is about 1 Tg P per year or on average 30 kg P per km2 per year. For a mean P content of soils of 0.5 kg per tonne of soil51,52 this means a soil loss of 60 ton km−2, which exceeds the European estimates by 50% due to larger erosivity of grasslands in especially tropical and (semi-)arid climates.

Estimates of erosion rates are uncertain. Erosion rates in overgrazed grasslands have been estimated to amount to 1,500 t of soil per km2 per year52,53. Such rates may locally be possible, but as a global average they seem unrealistic as they exceed the European average erosion rates of 40 t of soil per km2 per year46 by a factor of 40.

Weathering and atmospheric deposition

On the basis of the ratio of cropland area to world total land areas, Liu et al.52 calculated the amount of weathering and atmospheric deposition on croplands. Employing the same approach, we have estimated about 4.0 Tg P per year as weathering and atmospheric deposition in grasslands.

Scenario for the period 2006–2050

There are different global scenarios, such as from the Organization for Economic Cooperation and Development54, MEA28, IPCC-SRES55 and Rio+20 (ref. 21). For this study, we needed scenarios that focus on nutrients and included future food demand, production and the nutrient uptake. Particularly future production and consequently P uptake was needed as a target for selecting the appropriate scenarios. A few recent examples of scenarios that included projections of nutrient uptake are the MEA28 and more recently the Rio+20 (ref. 21) scenarios.

Future demand of P fertilizer in global croplands was calculated based on the target crop yields, given by four different MEA scenarios for 2050 (ref. 15).

Recently, van Vuuren et al.21 developed new pathways (Rio+20) to achieve global sustainability goals for food, land and biodiversity, as well as for energy and climate by 2050. The Rio+20 study describes four scenarios, that is, the Trend scenario and three challenge pathways. The Trend scenario describes possible trends in the absence of climate and sustainability policies. The three challenge pathways were designed to assess the potential to achieve sustainability goals.

We used the Rio+20 (ref. 21) Trend scenario for simulating future P requirement in grasslands. Baseline scenarios represent a continuation of current trends, with no marked changes or shifts in production and management systems, attitude towards environmental problems and so on. The Rio+20 Trend scenario is a baseline or business-as-usual scenario, and comparable with the baseline scenario of the Environmental Outlook of the Organization for Economic Cooperation and Development54, the Global Orchestration scenario of the MEA28 and the A1 scenario of IPCC-SRES55.

DPPS model description

DPPS is a simple two-pool P-model56 including labile and stable soil pools and long-term P input and output data15. The DPPS model reproduces historical grass uptake as a function of P inputs (fertilizer and manure).

P inputs are allocated to two dynamic P pools, namely the stable (P S ; 20%) and the labile P pools (P L ; 80%). The model simulates the P transfers between the pools, the uptake of P by the grass and the size of both pools. To calculate the dynamics of P in these two pools, two differential equations are used:

The rates of P transfer from P L to P S and vice versa are denoted by μ LS and μ SL , respectively (per year). The coefficient ρ refers to the total P input (mineral fertilizer and manure). The coefficients f and 1−f refer to the fraction of ρ that transfers to P L and P S , respectively. Coefficient α represents the grass P uptake fraction from P L . Parameters ω and δ are weathering and deposition inputs to P L and P S , respectively. The parameter stands for soil P erosion-runoff outflow from P L . A large μ LS makes P L less available for plant uptake and a large μ SL indicates that the stable pool acts as a buffer that replenishes the labile pool.

This model considers the essential P fluxes between grass and soil with a yearly time step. It can calculate P transfer between the labile and stable soil P pools, the P uptake by grass and the pool sizes. The model can also be formulated in a target-oriented approach57, in which the (future) P uptake is a model input and the required P application a result assuming no change in grassland area.

The DPPS model successfully simulated the historical patterns of crop P uptake as a response to the application rates in all continents and the entire globe as shown by Sattari et al.15 In the present paper, the model was used to calculate future P fertilizer and manure application rates in grasslands based on target grass productions in 2050 derived from the Rio+20 Trend (baseline) scenarios21. For more details of the model see Sattari et al.15 and Wolf et al.56

Comparing the results with other studies

The estimates of global P budgets of grassland soils and P transfers between grasslands and croplands are compared in Supplementary Table 5, with those reported by other authors. Although it is not easy to compare all results obtained in this study with those in other studies, due to differences in approaches, system boundaries and the scope, a comparison of selected flows calculated here at country level with other studies is possible (Supplementary Table 5). Frequently, in literature larger values than the ones reported here are found, primarily due to inclusion of pigs and poultry, while in our study only ruminants are considered.

Sensitivity analysis

The sensitivity of the results to 26-model parameters was investigated for 6 output variables for the soil P budget model, and the sensitivity of two output variables of the DPPS model to variation of 13 model parameters was investigated. These output variables represent global results for P uptake by grass (Supplementary Table 6). To limit computational load in the sensitivity analysis, the Latin Hypercube Sampling (LHS) technique58 was used. LHS offers a stratified sampling method for the separate input parameters, based on subdividing the range of each of the k parameters into disjunct equiprobable intervals based on a uniform distribution. The intervals were selected on the basis of earlier analysis of the livestock system. By sampling one value in each of the N intervals according to the associated distribution in this interval, we obtained N sampled values for each parameter. The number of runs N was 500 for the soil budget model and for the DPPS model.

The sampled values for the first model parameter are randomly paired to the samples of the second parameter, and these pairs are subsequently randomly combined with the samples of the third source and so on. This results in an LHS consisting of N combinations of k parameters. The parameter space is thus representatively sampled with a limited number of samples.

LHS can be used in combination with linear regression to quantify the uncertainty contributions of the input parameters to the model outputs58,59. The output Y considered (see columns in Supplementary Tables 7 and 8) is approximated by a linear function of the parameters X i expressed by

where β i is the so-called ordinary regression coefficient and e is the error of the approximation. The quality of the regression model is expressed by the coefficient of determination (R2), representing the amount of variation Y explained by Y−e. Since β i depends on the scale and dimension of X i , we used the standardized regression coefficient (SRC), which is a relative sensitivity measure obtained by rescaling the regression equation on the basis of the s.d.’s σ Y and :

SRC i can take values in the interval [−1, 1]. SRC is the relative change ΔY/σ Y of Y due to the relative change ΔX i / of the parameter X i considered (both with respect to their s.d. σ). Hence, SRC i is independent of the units, scale and size of the parameters. A positive SRC i value indicates that increasing a parameter value will cause an increase in the calculated model output, while a negative value indicates a decrease in the output considered caused by a parameter increase.

The sum of squares of SRC i values of all parameters equals the coefficient of determination (R2), which for a perfect fit equals 1. Hence, SRC i 2/R2 yields the contribution of parameter X i to Y. For example, a parameter X i with SRC i =0.1 adds 0.01 or 1% to Y in case R2 equals 1.

Here we discuss the SRC for 26 parameters in the soil P budget model (Supplementary Table 7) and 13 parameters in the DPPS model (Supplementary Table 8). We focus on significant SRC values exceeding 0.2, which we consider to be an important influence on the model sensitivity or an important parameter (that is, a contribution of 0.22=0.04 or about 4% to the variation of global results). Results for a similar analysis on the regional or smaller scale or different year would yield different results, depending on the P balances and local production system.

Soil budget model

The fraction of the livestock production in mixed and industrial systems (FrProdint) is very important for both P uptake in 2005 and cumulative P uptake for 1970 and 2005 in mixed respectively. pastoral systems, but not in the aggregated livestock system (Supplementary Table 7). Excretion rates and the number of animals (LPN) are important for uptake in 2005 and cumulative uptake, in each system (mixed and pastoral) and in the aggregated livestock system.

The manure deposited in grasslands, the amount that is diverted outside agricultural systems and the manure used as organic fertilizer in grasslands (represented by ‘FrGrazing’, ‘FrOthers’ and ‘FrGrass’, respectively) are not significant and not important for the sensitivity of uptake and cumulative uptake.

The kg meat or by-products per kg live weight (LWF), the DP—the ratio live weight to carcass—and the P contents of the various body parts (FPC) are all important and significant for uptake in 2005 and cumulative uptake in mixed systems and the aggregated livestock system, and for uptake in 2005 in pastoral systems. The parameter LWF is the fraction of each animal ‘product’: meat; bones; blood; and so on. Therefore, when LWF increases, P in products and thus the extraction will be higher (SRC positive). Higher values of DP lead to a lower amount of P found in the (by) products and hence less uptake. Livestock production is not only used in equation (1) but also in the calculation of feed (equation (8)). The higher the production induces an increase of feed use, since it is linked to productivity. Therefore, the increase in feed can exceed the increase in the P in products leading to a decrease in uptake (since uptake=P in products+P in manure–P in feed).

Dynamic Phosphorus Pool Simulator

Median simulated P uptake of global grasslands in 2005 is 15 Tg, with a variation of 11.0–17.6 Tg (2.5 and 97.5 percentile) as a result of variation of all the parameters in Supplementary Table 8. The cumulative P uptake (1970–2005) is 374–598 Tg around the median of 510 Tg. U U0 (initial P uptake from unfertilized soil) and U F0 (initial P uptake from fertilized soil) are the most important parameters for the sensitivity of all 13 parameters tested. U U0 , U F0 , α and β are all fitting parameters used to achieve the values for the size of LP and SP in year 0 of the simulation. U U0 (SRC=0.44; indicating sensitivity of 20%) and U F0 (SRC=0.64, indicating sensitivity of 41%) have a strong influence on the cumulative uptake of the total livestock system, and a smaller influence on the sensitivity of the uptake in 2005. The coefficient C u , the maximum fraction of crop P uptake from the labile pool is an important factor for the sensitivity of the cumulative uptake but has a smaller influence on the instantaneous uptake in 2005.

P input (ρ; manure plus mineral) has a stronger influence on the sensitivity of the instantaneous uptake in 2005 than on that of the cumulative uptake. This is due to depletion of the P in the two pools, which results in a strong and growing instantaneous effect of fertilizer, which is mostly flowing into the labile pool from which uptake occurs.