Combining livestock production information in a process-based vegetation model to reconstruct the history of grassland management

. Grassland management type (grazed or mown) and intensity (intensive or extensive) play a crucial role in the greenhouse gas balance and surface energy budget of this biome, both at ﬁeld scale and at large spatial scale. However, global gridded historical information on grassland management intensity is not available. Combining modelled grass-biomass productivity with statistics of the grass-biomass demand by livestock, we reconstruct gridded maps of grassland management intensity from 1901 to 2012. These maps include the minimum area of managed vs. maximum area of unmanaged grasslands and the fraction of mown vs. grazed area at a resolution of 0.5 ◦ by 0.5 ◦ . The grass-biomass demand is derived from a livestock dataset for 2000, extended to cover the period 1901–2012. The grass-biomass supply (i.e. forage grass from mown grassland and biomass grazed) is simulated by the process-based model ORCHIDEE-GM driven by historical climate change, rising CO 2 concentration, and changes in nitrogen fertilization. The global area of managed grassland obtained in this study increases from 6.1 × 10 6 km 2 in 1901 to 12.3 × 10 6 km 2 in 2000, although the expansion pathway varies between different regions. ORCHIDEE-GM also simulated augmentation in global mean productivity and herbage-use efﬁciency over managed grassland during the 20th century, indicating a gen-eral intensiﬁcation of grassland management at global scale but with regional differences. The gridded grassland management intensity maps are model dependent because they depend on modelled productivity. Thus speciﬁc attention was given to the evaluation of modelled productivity against a series of observations from site-level net primary productivity (NPP) measurements to two global satellite products of gross primary productivity (GPP) (MODIS-GPP and SIF data). Generally, ORCHIDEE-GM captures the spatial pattern, seasonal cycle, and interannual variability of grassland productivity at global scale well and thus is appropriate for global applications presented here.

Abstract. Grassland management type (grazed or mown) and intensity (intensive or extensive) play a crucial role in the greenhouse gas balance and surface energy budget of this biome, both at field scale and at large spatial scale. However, global gridded historical information on grassland management intensity is not available. Combining modelled grass-biomass productivity with statistics of the grassbiomass demand by livestock, we reconstruct gridded maps of grassland management intensity from 1901 to 2012. These maps include the minimum area of managed vs. maximum area of unmanaged grasslands and the fraction of mown vs. grazed area at a resolution of 0.5 • by 0.5 • . The grassbiomass demand is derived from a livestock dataset for 2000, extended to cover the period 1901-2012. The grassbiomass supply (i.e. forage grass from mown grassland and biomass grazed) is simulated by the process-based model ORCHIDEE-GM driven by historical climate change, rising CO 2 concentration, and changes in nitrogen fertilization. The global area of managed grassland obtained in this study increases from 6.1 × 10 6 km 2 in 1901 to 12.3 × 10 6 km 2 in 2000, although the expansion pathway varies between different regions. ORCHIDEE-GM also simulated augmentation in global mean productivity and herbage-use efficiency over managed grassland during the 20th century, indicating a general intensification of grassland management at global scale but with regional differences. The gridded grassland management intensity maps are model dependent because they depend on modelled productivity. Thus specific attention was given to the evaluation of modelled productivity against a series of observations from site-level net primary productivity (NPP) measurements to two global satellite products of gross primary productivity (GPP) (MODIS-GPP and SIF data). Generally, ORCHIDEE-GM captures the spatial pattern, seasonal cycle, and interannual variability of grassland productivity at global scale well and thus is appropriate for global applications presented here.

Introduction
The rising concentrations of greenhouse gases (GHGs), such as carbon dioxide (CO 2 ), methane (CH 4 ), and nitrous oxide (N 2 O), are driving climate change through increased radiative forcing (IPCC, 2013). It is estimated that, globally, livestock production (including crop-based and pasture-based) currently accounts for 37 and 65 % of the anthropogenic CH 4 and N 2 O emissions respectively (Martin et al., 2010;FAO, 2006). Grassland ecosystems support most of the world's livestock production, thus contributing indirectly a significant share of global CH 4 and N 2 O emissions. For CO 2 fluxes, however, grassland can be either a sink or a source with respect to the atmosphere. The annual changes in carbon storage of managed grassland ecosystems in Europe (hereafter referred to as net biome productivity, NBP) was found to be correlated with carbon removed by grazing and/or mowing (Soussana et al., 2007). Thus, knowledge of management type (grazed or mown) and intensity (intensive or extensive) is crucial for simulating the carbon stocks and GHG fluxes of grasslands.
For European grasslands, Chang et al. (2015a) constructed management intensity maps over the period 1961-2010 based on (i) national-scale livestock numbers from statistics (FAOSTAT, 2014), (ii) static sub-continental grass-fed fractions for each animal type (Bouwman et al., 2005), and (iii) the grass-fed livestock numbers supported by the net primary productivity (NPP) of the ORCHIDEE-GM (ORganizing Carbon and Hydrology In Dynamic Ecosystems grassland management) model. That study estimated an increasing NBP (i.e. acceleration of soil carbon accumulation) over the period 1991-2010. The increasing NBP was attributed to climate change, CO 2 trends, nitrogen (N) addition, and landcover and management intensity changes. The observationdriven trends of management intensity were found to be the dominant driver explaining the positive trend of NBP across Europe (36-43 % of the total trend with all drivers; Chang et al., 2016). That study confirmed the importance of management intensity in drawing up a grassland carbon balance. However, the national-scale management intensity and the identical history maps between 1901 and 1960 in that study carried several sources of uncertainty (Chang et al., 2015a). It implies that long-term history of large-scale gridded information on grassland management intensity is needed. The HYDE 3.1 land-use dataset (Klein Goldewijk et al., 2011) provides reconstructed gridded changes of pasture area over the past 12 000 years. Here, "pasture" represents managed grassland providing grass biomass to livestock. This reconstruction is based on population density data and countrylevel per capita use of pasture land derived from FAO statistics (FAOSTAT, 2008) for the post-1961 period and assumed by those authors for the pre-1960 period. It defines land used as pasture but does not provide information about management intensity. To our knowledge, global maps of grassland management intensity history are not available.
Recently, Herrero et al. (2013) garnered global livestock data to create a dataset with gridded grass-biomass-use information for year 2000. In this dataset, grass used for grazing or silage is separated from grain feed, occasional feed, and stover (fibrous crop residues). A variety of constraints have been taken into account in creating this global dataset, including the specific metabolisable energy (ME) requirements for each animal species and regional differences in animal diet composition, feed quality, and feed availability. This grass-biomass-use dataset provides a starting point for constraining the amount of carbon removed by grazing and mowing (i.e. the target of grass-biomass use) and is suitable for adoption by global vegetation models to account for livestock-related fluxes.
The major objective of this study is to produce global gridded maps of grassland management intensity since 1901 for  (Chang et al., 2013(Chang et al., , 2015b with gridded grass-biomass use extrapolated from Herrero et al. (2013). First, ORCHIDEE-GM is calibrated to simulate the distribution of "potential" (maximal) harvested and grazed biomass from mown and grazed grasslands respectively. In a second step, the modelled productivity maps are used in combination with livestock data to reconstruct annual maps of grassland management intensity, at a spatial resolution of 0.5 • by 0.5 • . This is done for each country since 1961 and for 18 large regions of the globe for 1901-1960. The reconstructed management intensity defines the fraction of mown, grazed, and unmanaged grasslands in each grid cell. The gridded grassland management intensity maps are model dependent because they rely on simulated NPP. Thus, in this study we also give a specific attention to the evaluation of modelled productivity against both a new set of site-level NPP measurements and satellite-based models of gross primary productivity (GPP). In Sect. 2, we describe the ORCHIDEE-GM model, the adjustment of its parameters for the C4 grassland biome, model input, the method proposed to reconstruct grassland management intensity, and the data used for evaluation. The derived management intensity maps and the comparison between modelled and observed productivity are presented in Sect. 3 and discussed in Sect. 4. Concluding remarks are made in Sect. 5.
2 Material and methods

Model description
ORCHIDEE is a process-based ecosystem model developed for simulating carbon fluxes, and water and energy fluxes in ecosystems, from site level to global scale (Krinner et al., 2005;Ciais et al., 2005;Piao et al., 2007). ORCHIDEE-GM (Chang et al., 2013) is a version of ORCHIDEE that includes the grassland management module from PaSim (Riedo et al., 1998;Vuichard et al., 2007a, b;Graux et al., 2011), a grassland model for field-level to continental-scale applications. Accounting for the management practices such as mowing, livestock grazing and organic fertilizer application on a daily basis, ORCHIDEE-GM proved capable of simulating the dynamics of leaf area index, biomass, and C fluxes of managed grasslands. ORCHIDEE-GM version 1 was evaluated and some of its parameters calibrated, at 11 European grassland sites representative of a range of management practices, with eddy-covariance net ecosystem exchange and biomass measurements. The model successfully simulated the NBP of these managed grasslands (Chang et al., 2013). Chang et al. (2015b) then added a parameterization of adaptive management through which farmers react to a climatedriven change of previous-year productivity. Though a full N cycle is not included in ORCHIDEE-GM, the positive effect of nitrogen fertilizers on grass photosynthesis rates, and thus on subsequent ecosystem productivity and carbon storage, is parameterized with an empirical function calibrated from literature estimates (version 2.1; Chang et al., 2015b). ORCHIDEE-GM v2.1 was applied over Europe to calculate the spatial pattern, interannual variability (IAV), and the trends of potential productivity, i.e. the productivity that maximizes simulated livestock densities assuming an optimal management system in each grid cell (Chang et al., 2015b). This version was further used to simulate NBP and NBP trends over European grasslands during the last 5 decades at a spatial resolution of 25 km and a 30 min time step (Chang et al., 2015a). ORCHIDEE-GM v1 and v2.1 were developed based on ORCHIDEE v1.9.6. To benefit from recent developments and bug corrections in the ORCHIDEE model, ORCHIDEE-GM is updated in this study with ORCHIDEE Trunk.rev2425 (available at https://forge.ipsl.jussieu.fr/orchidee/browser/ trunk#ORCHIDEE). We further made the adjustment of its parameters for the C4 grassland biome (Sect. 2.2) and implemented a specific strategy for wild herbivores grazing (Sect. 2.3; also see Supplement Sect. S1). The updated model is referred to hereafter as ORCHIDEE-GM v3.1.

Model parameter settings
ORCHIDEE-GM was applied to simulate GHG budgets and ecosystem carbon stocks under climate, CO 2 , and management changes for Europe. However, an extension of model application to regions outside Europe requires first a calibration of key productivity-related parameters. Two sensitive parameters representing photosynthetic capacity (the maximum rate of Rubisco carboxylase activity at a reference temperature of 25 • C; Vc max 25) and the morphological plant traits (the maximum specific leaf area; SLA max ) were reported by Chang et al. (2015a) for simulating grassland NPP. The Vc max 25 = 55 µmol m −2 s −1 and SLA max = 0.048 m 2 per g C in ORCHIDEE-GM were previously defined from observations and indirectly evaluated against eddy-flux tower measurements of GPP for temperate C3 grasslands in Europe (Chang et al., 2013(Chang et al., , 2015b. The global TRY database gives SLA values for C4 grasses, of 0.0192 m 2 g −1 dry matter (DM) (0.0403 m 2 per g C with a mean leaf carbon content per DM of 47.61 %; Kattge et al., 2011). Thus, we have set the value of SLA max = 0.044 m 2 per g C for C4 grasses in ORCHIDEE-GM to fit the mean value from the TRY estimate, as we did previously for C3 grasses (Chang et al., 2013). The parameter Vc max 25 cannot be directly measured, but it is usually derived from A/C i curves in C3 or C4 photosynthesis models (C3: Farquhar et al., 1980;C4: Collatz et al., 1992), where A is the leaf-scale net CO 2 assimilation rate and C i the partial pressure of CO 2 in leaf intercellular spaces. Several researches provide observation-based estimates of Vc max 25 (Feng and Dietze, 2013;Verheijen et al., 2013; range of 24-131 µmol m −2 s −1 for C3 grasses and of 15-46 µmol m −2 s −1 for C4 grasses). Based on these estimates, we keep the value of Vc max 25 = 55 µmol m −2 s −1 previously calibrated in Europe for all C3 grasses and set Vc max 25 = 25 µmol m −2 s −1 for C4 grasses. These values may reflect neither differences in nitrogen and phosphorus availability between locations nor adaptation or species changes within a C3 or C4 grassland, but they are within the range of observations made under different conditions and consistent with values used by other terrestrial ecosystem models (Table S1 in the Supplement). All other parameters of ORCHIDEE model are kept the same as in Trunk.rev2425. The parameter settings for grassland management module are in consistent with that in ORCHIDEE-GM v1 (Chang et al., 2013) and v2.1 (Chang et al., 2015a, b).

Model input
ORCHIDEE-GM v3.1 was run on a global grid over the globe using the 6-hourly CRU+NCEP reconstructed climate data at 0.5 • × 0.5 • spatial resolution for the period 1901-2012 (Viovy, 2013). The fields used as input of the model are temperature, precipitation, specific humidity, solar radiation, wind speed, pressure, and long-wave radiation. Other input data are (1) yearly domestic grazing-ruminant stocking density maps, (2) wild-herbivores population density maps, (3) N fertilizer application maps including manure-N and mineral-N fertilizers, and (4) atmospheric-N deposition maps. These input maps all cover the period from 1901 to 2012 and are briefly described below (also see Supplement Sects. S2-S5). Table 1 lists all variables shown in this section, including their abbreviations, units, related equations, and data sources.
Grazing-ruminant stocking density maps: spatial statistical information on grazing-ruminant stocking density is not available at global scale. In this study, we combined the domestic ruminant stocking density maps (Supplement Sect. S2) and historic land-cover change maps (Supplement Sect. S3) to construct gridded grazing-ruminant stocking density.
Assuming that all the ruminants in each grid cell were grazing on the grassland within the same grid, we defined the grazing-ruminant stocking density in grid cell k in year m (D grazing,m,k , livestock unit (LU) per ha of grassland area) as where D m,k is the total domestic ruminant stocking density (unit: LU per ha of land area; Supplement Sect. S2) and f grass,m,k is the grassland fraction in grid cell k in year m from a set of historic land-cover-change maps (Supplement Sect. S3). To avoid unrealistic densities of ruminant grazing over grassland (which might cause grasses to die during the growing season), a maximum value of 5 LU ha −1 was set for the density map. In addition, a minimum grazing-ruminant density of 0.2 LU ha −1 was set to avoid economically implausible stocking rates. Figure S1 in the Supplement shows the example maps of domestic ruminant stocking density (D) and the corresponding grazing-ruminant stocking density (D grazing ) for reference year 2006. Wild herbivore density maps: gridded maps of wild herbivore density are not available; therefore the gridded population density of wild herbivores (D wild ; unit: LU per ha of grassland area) is derived from the literature data and from Bouwman et al. (1997) (see Table S2 for detail). The population of these herbivores from literature was first converted to LU according to the ME requirement calculated from their mean weight (Table S2) and then distributed to suitable grasslands based on grassland aboveground (consumable) NPP simulated from ORCHIDEE-GM v3.1 (Supplement Sect. S4; Fig. S2). The wild herbivores density was assumed to remain constant during the period of 1901-2012, because no worldwide historical wild-animal population information was available. A specific grazing strategy for wild herbivores is incorporated in the model (Supplement Sect. S1). We assumed wild herbivores eat fresh grass biomass during the growing season and eat dead grass during the non-growing season.
Nitrogen application rates from mineral fertilizers and manure: grassland is fertilized with organic N fertilizer (e.g. manure, slurry) and/or even mineral-N fertilizer, though this is not as common as for cropland. Gridded fertilizer application rates on grassland are not available worldwide. The only exception that we are aware of is for European grasslands (Leip et al., 2008(Leip et al., , 2011(Leip et al., , 2014; data available for EU-27 as used in Chang et al., 2015a). For countries/regions other than EU-27, the following data were used. The amount of manure-N fertilizer for 17 world regions at 1995 was derived from various sources (e.g. IFA, 1999;FAO/IFA/IFDC, 1999;FAO/IFA, 2001) and synthesized by Bouwman et al. (2002a, b; Table S3). For mineral-N fertilizers on grassland, countryscale data of fertilized area and mean fertilization rate for 1999/2000 are available in FAO/IFA/IFDC/IPI/PPI (2002) with grassland/pasture been fertilized in 13 non-EU countries. The regional/country-scale data were downscaled to a 0.5 • × 0.5 • grid and extended to cover the period 1901-2012 (see Supplement Sect. S5 for detail).
Atmospheric-nitrogen deposition maps: the historical atmospheric-N deposition maps were simulated by the LMDz-INCA-ORCHIDEE global chemistry-aerosolclimate model (Hauglustaine et al., 2014). Hindcast simulations for the years 1850, 1960, 1970, 1980, 1990, and 2000 have been performed using anthropogenic emissions from Lamarque et al. (2010). The total nitrogen deposition fields (wet and dry; NH x and NO y ) of all nitrogen-containing gas-phase and aerosol species have been simulated at a spatial resolution of 1.9 • in latitude and 3.75 • in longitude. Linear interpolation was performed between the hindcast snapshot years to produce temporally variable  (2) this study Y graze Annual potential biomass consumption over grazed grasslands this study f unmanaged Maximum fraction of unmanaged grassland , (11) this study a The subscripts of these variables in this study: i is ruminant category; j is country; k is grid cell; m is year; q is region. b When not specified, the ha −1 (or m −2 ) in the units indicate per ha (or per m 2 ) of grassland area.
atmospheric-N deposition maps (N deposition , unit: kg N per ha of grassland area per year).

Simulation set-up
Considering different photosynthetic pathways and management types, six grassland plant functional types (PFTs) are defined: C3 natural (unmanaged) grassland, C3 mown grass-land, C3 grazed grassland, C4 natural (unmanaged) grassland, C4 mown grassland, and C4 grazed grassland. In the simulation, we ideally consider that grassland PFTs are distributed all over the world. Post-processing will incorporate the information of grassland distribution in the real world (Supplement Sect. S3  Figure 1. Illustration of the procedures for reconstructing management intensity maps. Italic texts indicate the major steps of the reconstruction. The meanings, units, related equations, and data sources of the variables (i.e. gridded maps) are shown in Table 1. D grazing is grazing-ruminant stocking density; D wild is wild herbivore density; N manure is organic (manure) nitrogen fertilizer application rate; N mineral is mineral-nitrogen fertilizer application rate; N deposition is atmospheric-nitrogen deposition rate; Y mown is annual potential harvested biomass from mown grasslands; Y graze is annual potential grazed biomass from grazed grasslands; GBU is grass-biomass use; f mown is minimum fraction of mown grassland; f grazed is minimum fraction of grazed grassland; f unmanaged is maximum fraction of unmanaged grassland.
(N deposition ). For each grassland PFT, specific forcing and management strategies are used (summarized in Fig. 1). Unmanaged grasslands are forced by wild herbivore density maps (D wild ). Both mown and grazed grasslands are forced by the historical N fertilizer maps described above, which include manure (N manure ) and mineral fertilizers (N mineral ). Grazed grassland is additionally forced by the historical gridded grazing-ruminant stocking density (D grazing ).

Grassland management intensity and historical changes
Figure 1 briefly illustrates the procedures of combining model output, grass-biomass-use data, and grassland area data to reconstruct grassland management intensity maps. This section presents the procedures of the reconstruction in detail. Table 1 lists all variables shown in this section, including their abbreviations, units, related equations, and data sources. Herrero et al. (2013) established a global livestock production dataset containing a high-resolution (8 km × 8 km) gridded map of grass-biomass use for the year 2000. In this study, this dataset is extrapolated annually over 1901-2012 to constrain the grass-biomass consumption in ORCHIDEE-GM v3.1. Assuming that grass-biomass use for grid cell k in country j and year m (GBU m,j,k , unit: kg DM per year) varies proportionally with the total ME requirement of domestic ruminants in each country, GBU m,j,k can be calculated from its value of the year 2000 given by Herrero et al. (2013), according to where D m.k and D 2000,k are the total ruminant stocking density for grid cell k in year m and in year 2000 calculated by Eqs. (S4) and (S5) in Supplement Sect. S2, which take into account the changes in category-specific ME requirement at country scale (1961-2012) or regional scale . ORCHIDEE-GM v3.1 simulates the annual potential (maximal) harvested biomass from mown grasslands (Y mown , unit: kg DM m −2 yr −1 from mown grassland) and the annual potential biomass consumption per unit area of grazed grassland (Y grazed , unit: kg DM m −2 yr −1 from grazed grassland) in each grid cell. Under mowing, the frequency and magnitude of forage harvests in each grid cell is a function of grown biomass (Vuichard et al., 2007a). The effective yield on grazed grassland (i.e. Y grazed ) depends on the grazing stocking rate (here, D grazing ) and on the environmental conditions of the grid cell (Chang et al., 2015a); it is calculated as where IC is the daily intake capacity for 1 LU (∼ 18 kg DM per day calculated in Supplement Sect. S1 of Chang et al., 2015b), and T grazing,m,k is the number of grazing days in grid cell k at year m. Due to the impact of livestock on grass growth through trampling, defoliation (i.e. biomass intake), etc., and because grassland cannot be continuously grazed during the vegetation period, thresholds of shoot biomass are set for starting, stopping, and resuming grazing (Vuichard et al., 2007a). The "recovery" time required under grazing is obtained in the model using threshold (Vuichard et al., 2007a;Chang et al., 2015a), which determines when grazing stops (dry biomass remaining lower than 300 kg DM ha −1 ) or when grazing can start again (dry biomass recovered to a value above 300 kg DM ha −1 for at least 15 days). Y grazed is usually lower than Y mown in temperate grasslands due to the lower herbage-use efficiency of grazing simulated by ORCHIDEE-GM (Chang et al., 2015b). However, in some arid regions the grass biomass does not grow enough during the season to trigger harvest; i.e. it does not reach the threshold in the model at which farmers are assumed to decide to cut grass for feeding forage to animals (see Chang et al., 2015b), so that Y grazed can become larger than Y mown (Fig. S3). The following set of rules was used to reconstruct historical changes in grassland management intensity, based on NPP simulated by ORCHIDEE-GM v3.1.
-Rule 1: for each grid cell and year, the total biomass removed by either grazing and cutting must be equal to the grass-biomass use, GBU m,k .
Thus, for grid cell k in year m, the minimum fraction of grazed (f grazed,m,k ), the minimum fraction of mown (f mown,m,k ), and the maximum fraction of unmanaged grassland (f unmanaged,m,k ) are calculated with the following equations (definitions of minimum and maximum in this context are given below). If where A grass,m,k (unit: m 2 ) is the grassland area for grid cell k in year m of the series of historic land-cover change maps If GBU m,k cannot be fulfilled by any combination of modelled Y grazed and Y mown , we diagnose a modelled grassbiomass production deficit and apply the following equations: if Y grazed >Y mown , then f grazed,m,k = 1, f mown,m,k = 0, and f unmanaged,m,k = 0, if Y grazed <Y mown , then f mown,m,k = 1, f grazed,m,k = 0, and f unmanaged,m,k = 0.
This set of equations is valid for a mosaic of different types of grasslands in each grid cell, some managed (grazed and/or mown) and some remaining unmanaged. In reality, (1) farm owners could increase the mown fraction to produce more forage, which corresponds approximately to the mixed and landless systems of Bouwman et al. (2005); and (2) animals could migrate a long way across grazed and unmanaged fractions (as they do in real rangelands) and only select the most digestible grass in pastoral systems, which corresponds to extensively grazed grasslands. Yet, given the approximations made in this study, f grazed,m,k and f mown,m,k represent the minimum fractions of grazed/mown grasslands rather than the actual fractions, and f unmanaged,m,k corresponds to a maximum fraction of unmanaged grasslands since both mixed and landless and extensive grazing are not modelled.
Herbage-use efficiency (Hodgson, 1979) is defined as the forage removed expressed as a proportion of herbage growth. It can be an indicator of management intensity over managed grassland, in addition to the fraction of managed area obtained above. In this study, the forage removed is modelled annual grass-biomass use including Y grazed and Y mown , and herbage growth is modelled annual grass GPP.

Model evaluation: datasets and model-data agreement metrics
The gridded grassland management intensity maps are model dependent because they depend on modelled productivity. Thus the evaluation of modelled productivity becomes necessary. In this study, modelled productivity (NPP and GPP) is compared to a new set of site-level NPP measurements (Sect. 2.7.1) and two satellite-based models of GPP (MODIS-GPP, from the Moderate Resolution Imaging Spectroradiometer, Sect. 2.7.2; sun-induced chlorophyll fluorescence (SIF) data, Sect. 2.7.3). Modelled NPP (or GPP) combines grassland productivity of all PFTs (Sect. 2.4), accounting for the variable fractions of grazed, mown, and unmanaged grassland in each grid cell calculated by Eqs. (4-11), and hereafter is referred to as NPP model (or GPP model ).
Model-data agreement of NPP and GPP was assessed using Pearson's product-moment correlation coefficients (r) and root mean squared errors (RMSEs).

Grassland NPP observation database
NPP is a crucial variable in vegetation models and it is essential that this variable is properly validated. High-quality measurements of grassland NPP are scarce, partly due to the difficulty of measuring some NPP components such as fineroot production (Scurlock et al., 1999(Scurlock et al., , 2002. An updated version of the Luyssaert et al. (2007) database comprising non-forest biomes (Campioli et al., 2015) was used here. This database contains a flag indicating managed or unmanaged to each site and provides mean annual temperature, annual precipitation, and downwelling solar radiation based on site measurements from the literature, CRU database (Mitchell and Jones, 2005), MARS database (http://mars.jrc.ec. europa.eu/mars/About-us/AGRI4CAST/Data-distribution/ AGRI4CAST-Interpolated-Meteorological-Data), or World-Clim database (Hijmans et al., 2005). Three additional datasets used in this study present NPP measurements from 30 sites across China (Zeng et al., 2015;Y. Bai, personal communication, 2015) and 16 sites across western Siberia (Peregon et al., 2008; with data updated to 2012). Data from China include NPP observations at fenced (i.e. unmanaged) and unfenced (i.e. managed) grassland for each site, and data of western Siberia are observations from natural wetland.
In total, we have 305 NPP observations (NPP obs ) with separated aboveground and belowground NPP from 129 sites all over the world (including grassland, wetland, and savanna; Fig. S4). Duplicate observations from the same site year were averaged and considered as a single entry. NPP measurements with different management (managed or unmanaged) at the same site were considered as identical observations. In total, 270 grassland NPP measurements were compared to the simulation of ORCHIDEE-GM v3.1 for the grid cell, corresponding to each site and for the time period of observation. Depending on the status of measured grassland (unmanaged or managed), modelled NPP from unmanaged or managed grassland is used for comparison. Modelled NPP over managed grassland accounts for the NPP from mown and grazed grassland and their corresponding fractions.

Grassland GPP from MODIS products
The MOD17A3 dataset (version 55; Zhao et al., 2005;Zhao and Running, 2010) -a MODIS product on vegetation production -provides the seasonal and annual GPP data at a spatial resolution of 1 km from 2000 to 2013.
To obtain the grassland GPP from the MOD17 dataset, we first extract the MOD17 GPP at 1 km resolution over grassland grids in the MOD12Q1 dataset. Here, the grassland in the MOD12Q1 dataset includes the "open shrubland", "savanna", and "grassland" in the Boston University's UMD classification scheme. The extracted annual and seasonal MODIS GPP was then averaged and aggregated to 0.5 • × 0.5 • spatial resolution to be comparable to model output.

Sun-induced chlorophyll fluorescence data
Space-based observations of SIF provide a time-resolved measurement of a proxy of photosynthesis . Similar to the MPI-BGC data-driven GPP product (Jung et al., 2011), SIF values exhibit a linear relationship (r 2 = 0.79) with monthly tower GPP at grassland sites in western Europe . Compared to MODIS EVI (MOD13C2 products), SIF observations drop to zero during the non-growing season, thus providing a clearer signal of photosynthetic activity  than other vegetation indices based on visible and near-infrared reflectances. SIF also provides a better seasonal agreement with GPP from flux towers as compared to vegetation indices .
In this study, we used monthly GOME-2 (version 26, level 3) SIF data products with the spatial resolution of 0.5 • × 0.5 • (available from 2007 to 2012). SIF-GPP is calculated by a SIF-GPP linear model adjusted from Guanter et al. (2014) (SIF-GPP = −0.1 + 4.65 × SIF (V26); see Supplement Sect. S6 for detail). To reduce the contamination of SIF by non-grassland PFTs, we restrict the model-data comparison to grassland-dominated grid cells, defined as those with grassland cover in the MOD12Q1 dataset (Sect. 2.5.2) is larger than 50 %. Figure 2 shows the minimum fractions of mown and grazed grasslands and the maximum fraction of unmanaged out of total grassland (f mown , f grazed , and f unmanaged respectively; Sect. 2.4) in the year 2000. Grazed grasslands comprise most of the managed grasslands in the maps (Fig. 2b). Significant fractions of mown grasslands are only found in regions with high ruminant stocking density such as eastern China, India, eastern and northern Europe, and eastern United States, where Y grazed cannot fulfil the grass-biomass demand (Fig. 2a). Using the FAO-defined regions (see caption to Table 2), the largest fractions of managed grasslands are modelled in regions of high ruminant stocking density (Fig. S1) such as in eastern Europe with a mean fraction of 90 ± 17 % (the mean is the average fraction of mown and grazed grasslands over all the grid cells in this region, and the standard deviation is taken from differences between grid cells), South Asia (59 ± 46 %), and western Europe (55 ± 36 %). The lowest managed grasslands fractions are modelled in the Russian Federation (17 ± 34 %).

Maps of grassland management intensity
In some grid cells, the simulated grassland productivity is not sufficient to fulfil the grass-biomass use given by Herrero et al. (2013; Fig. 2d). Of the 2.4 billion tonnes of grassbiomass use (in dry matter for the reference year 2000) given by Herrero et al., 16 % cannot be fulfilled by the productivity simulated by ORCHIDEE-GM v3.1. This translates into a modelled grass-biomass production deficit of 0.38 billion tonnes ( Table 2). Out of all regions, the largest modelled production deficit (f global in Table 2) is found in South Asia (49 %). This South Asian deficit is predominantly in India (35 % of the modelled global total deficit) and Pakistan (10 % of the modelled global total deficit). Other regions with a biomass production deficit are the Near East and North Africa (18 %) and sub-Saharan Africa (13 %). Overall, 32 % of the global production deficit comes from regions with dry climate and low NPP (less than 50 g C m −2 yr −1 ), and 34 % of it comes from regions with low grassland cover (less than 10 % of total land cover). The causes of this grass-biomass production deficit diagnosed by ORCHIDEE-GM are discussed in Sect. 4.2.
Modelled herbage-use efficiency over managed grassland during the 2000s (grazed plus mown; Fig. 3) ranges between 2 and 20 % in most regions and generally follows the spatial pattern of grazing-ruminant density (Fig. S1). High herbageuse efficiency (over 20 %) is found in regions with significant mown grassland (f mown ) simulated due to the larger fraction of biomass removed over mown grassland than that over grazed grassland in the same grid cell (Fig. S3). Figure 4 displays the NPP per unit area and the production (Prod = NPP × grassland area) of each type of grassland for 10 FAO-defined regions and the globe. Even when grassland management is included, the production of unman-  aged grassland (Prod unmanaged ) still comprises 63 % of the total production (Prod total ) in the 1990s. The production of grazed grasslands (Prod grazed ) accounts for 34 % of Prod total , while the production of mown grasslands (Prod mown ) is only 3 %, given the small area under this management practice (Fig. 4). Mown grasslands only contribute to production in the regions where climate conditions and fertilizers maintain a high NPP, and Y grazed is not enough to fulfil the animal requirement, which triggers the harvest practice in Eqs. (7-11). Over unmanaged grassland (Fig. S2), ORCHIDEE-GM v3.1 simulated a total annual consumption by wild herbivores of 147-654 million tonnes DM of the 5778 million tonnes DM in aboveground NPP (consumable NPP) over suitable grassland (Table S5), which comprises 3-11 % of the consumable NPP, similar to the range given by Warneck (1988). The fraction of consumption in consumable NPP varied from 1 % in the former USSR to 9 % in Scandinavia, indicating the different significance of wild herbivores on grassland.

Historical changes in the area and productivity of managed grassland
The global minimum area of managed grassland (A managed-gm ) is of 6.1 × 10 6 km 2 in 1901 and increased to 12.3 × 10 6 km 2 in 2000 (Table 3; Fig. 5) -an increase of The potential harvested biomass from mown grassland (Ycut) and the potential biomass consumption over grazed grassland (Ygraze) are 10-year averages for the period 1901-1910 (1900s) and 1991-2000 (1990s), representing the productivity at the beginning and at the end of the 20th century respectively. c Ruminant numbers (in units of livestock unit, LU) are calculated based on the total metabolisable energy (ME) requirement by all ruminant. The ME requirement by all ruminants is based on ruminant numbers from statistics (for 1961-2012; data derived from FAOSTAT, 2014) and literature estimates (for 1901-1960data derived from Mitchell (1993data derived from Mitchell ( , 1998a and available in HYDE database at http://themasites.pbl.nl/tridion/en/themasites/hyde/landusedata/livestock/index-2.html), using the calculation method given in the Supplement Sect. S1 of Chang et al. (2015a).   (Hodgson, 1979) is defined as the forage removed expressed as a proportion of herbage growth. In this study, the forage removed is modelled annual grass-biomass use including Y grazed and Y mown , and herbage growth is modelled annual grass GPP.
102 % during the 20th century. This expansion of managed grasslands is mainly explained by the increase in the area of grazed lands (+5.7 × 10 6 km 2 ), while mown grassland increased only marginally (+0.5 × 10 6 km 2 ). The largest extension of A managed-gm is found in sub-Saharan Africa (+1.8 × 10 6 km 2 ) and Latin America and the Caribbean (+1.7 × 10 6 km 2 ; Fig. 5). The regions with the largest relative expansion of managed grasslands (as a percentage of 1901 areas) are sub-Saharan Africa (+219 %), East and Southeast Asia (+204 %), nd Latin America and the Caribbean (+175 %), and the regions where the number of domestic ruminants (N ruminant ) increased by nearly or over a factor of 3. Only small increases of A managed-gm were modelled in western Europe (+41 × 10 3 km 2 ; i.e. 8 %) and eastern Europe (+27 × 10 3 km 2 ; i.e. 8 %), despite an increase of N ruminant by a factor of 1.5 in western Europe (+27×10 6 LU) and of 1.4 in eastern Europe (+5 × 10 6 LU). This means that livestock production intensified in those two regions, first by giving crop feedstock given to animals (Bouwman et al., 2005) and second through the optimization of forage harvesting and grazing to feed higher animalstocking densities. Note that the animal density in eastern and western Europe peaked at 123 × 10 6 LU near 1990 and has declined by 29 % since then. Besides the extension of managed grassland area, modelled herbage-use efficiency over managed grassland increased from 6.2 to 6.6 % during the 20th century, indicating the intensification of grassland management. Large increase in herbage-use efficiency is modelled in South Asia (+3.6 %) and eastern Europe (+2.7 %), while marginal decrease of herbage-use efficiency is found in the Near East and North Africa (−0.1 %) and Oceania (−0.2 %; Table 3).
The global mean potential productivity of mown grassland (Y mown ) increased by 62 % from 0.29 kg DM m −2 yr −1 for 1900s to 0.48 kg DM m −2 yr −1 for the 1990s, while that of grazed grassland Y grazed increased by 40 %, from 0.10 kg DM m −2 yr −1 for the 1900s to 0.14 kg DM m −2 yr −1 for the 1990s (Table 3). During the last century, Y mown increased by more than 40 % in most regions except in Latin America and the Caribbean (14 %), while the increase of Y grazed ranged from 25 % in sub-Saharan Africa and 80 % in eastern Europe (Table 3). Figure 6 shows the grassland productivity (NPP model ; Fig. 6a) and the NPP differences between NPP model and NPP from unmanaged grassland (Fig. 6b). The effect of including management does not produce a big difference in simulated NPP, which has similar patterns in most regions (Fig. 6b) (Tables S3 and S4) cause a higher productivity (Fig. 6b). Figure 7a shows the comparison between site-scale NPP observations (NPP obs ) and the model results at the corresponding grid cells (NPP model ). The NPP model is positively correlated with NPP obs across 129 sites but with the low cor- relation coefficient of r = 0.35 (p<0.01) and the RMSE of 380 g C m −2 yr −1 . Figure 7b presents box-and-whisker plot of the observed and modelled annual whole-plant NPP, aboveground NPP, and belowground NPP. The mean value and range of modelled whole-plant NPP are both higher than those of NPP obs . The NPP overestimation by the model is mainly due to a too-high aboveground NPP, while belowground NPP is only little higher for its mean or even lower for its median than belowground NPP obs .

Evaluation of modelled GPP against MODIS-GPP for annual mean and interannual variability
At global scale, MODIS-GPP gives a mean grassland GPP of 537 g C m −2 yr −1 and ORCHIDEE-GM v3.1 simulates a mean value of 796 g C m −2 yr −1 , ≈ 50 % higher than MODIS-GPP. A higher modelled GPP (GPP model ) than MODIS is found for all latitude bands especially in boreal (50-80 • N) and tropical regions (20 • S-20 • N; Fig. 8).
The linear regression between gridded MODIS-GPP and GPP model suggests a similar spatial pattern (slope = 1.05, and the correlation coefficient r spatial = 0.84; Fig. S5). The temporal correlation coefficient between the detrended time series of global GPP model and MODIS-GPP was found to be high (r IAV-global 0.88, p<0.01). Within the grid  Peel et al. (2007) using climate data from WorldClim (http://www. worldclim.org/). In subplot (b), the "whisker" indicates the crossmeasurement (total 270 measurements) uncertainty. cells covered by grass over more than 20 % of total land in MOD12Q1, significant positive interannual correlations between GPP model and MODIS-GPP were found for 39 % of the grid cells (i.e. 40 % of the grassland area), except in some tundra areas of Siberia and North America, grassland on the Qinghai-Tibet Plateau, and savannah in sub-Saharan Africa (Fig. 9).

Evaluation of modelled seasonal cycle of GPP
against MODIS-GPP and GOME-2 SIF products Figure 10 compares the normalized seasonal variation of GPP model , MODIS-GPP, and SIF-GPP for five latitude bands and the globe. Similar mean seasonal variations of grassland productivity are found between modelled GPP, MODIS-GPP, and SIF (r seasonal range from 0.55 to 0.89; Table 4). Compared to both MODIS-GPP and SIF data, ORCHIDEE-GM v3.1 captures the seasonal variation of productivity in bo-Biogeosciences, 13, 3757-3776, 2016 www.biogeosciences.net/13/3757/2016/ Table 4. Mean ± standard deviation of r seasonal comparing the seasonal cycle of modelled GPP (GPP model ), MODIS-GPP, and SIF data for the five latitude bands and global scale. r seasonal is expressed as mean ± standard deviation of grid level correlation coefficient within each latitude band and global. To avoid the strong impact of other land-cover types (e.g. crop and forest) to the seasonal cycle, we only consider r seasonal for grid cells with grassland covering more than 50 % of total land in the MOD12Q1 dataset.  real and temperate regions of the Northern Hemisphere well (r seasonal >0.8; Table 4). In the band from 60 • S to 30 • N, relatively low average r seasonal correlations are found both with MODIS-GPP and SIF (ranging from 0.55 to 0.71). However, note that the r seasonal between the two remote sensing GPP related products is relatively low for grassland between 60 • S and 30 • N, particularly between 0 and 60 • S (Table 4).

Managed area of grassland and management intensity: comparison with previous estimates
The area of managed grasslands obtained in this study is lower than the pasture area of HYDE 3.  -STAT, 2008). In this study, A managed-gm is constrained by grass-biomass-use data (i.e. requirement of biomass for animals) and the simulated grassland productivity (i.e. supply of biomass to animals). In fact, the actual (real-world) managed grassland area could be larger than A managed-gm in regions where grasslands are not strictly unmanaged, i.e. not fully occupied by A managed-gm in the management intensity maps (i.e. f unmanaged >0; Fig. 2c). In pastoral systems such as open rangeland and mountain areas, animals keep moving to search for the most digestible grass. Tracts of grasslands can be grazed for a short period, with only a small part of the annual grass productivity being digested (i.e. very low herbageuse efficiency). This type of grassland could be recognized as extensively grazed grassland, whereas it is considered as unmanaged in this study. For example, lower herbage-use efficiency than that simulated in this study (Fig. 3)  Reclassifying these areas would result in a larger area of extensively managed grassland. Few studies reported the herbage-use efficiency of managed grassland. One exception is the network of European eddy-covariance flux sites. For these sites the average herbage-use efficiency (expressed as forage defoliated as a proportion of GPP) is 7.1 % ± 6.1 % for grazed sites, and 13.3 % ± 6.4 % for mown sites (J.-F. Soussana, personal communication, 2015); a similar range, between 2 and 20 %, is simulated in this study (Fig. 3). The time evolution of A managed-gm since 1901 in this study is arguably more realistic than HYDE because it considers changes in animal stocking density from statistics and the evolution in per-head use of pasture. A managed-gm takes into account (1) changes in grass-biomass requirement, considering both ruminant numbers and meat/milk productivity (Supplement Sect. S2; N ruminant in Table 3); (2) changes in grassland productivity driven by climate change, rising CO 2 concentration, and changes in N fertilization (Y mown and Y grazed in Table 3); and (3) changes in management types (mown and grazed grassland areas in Table 3 and Fig. 5). For example in intensively managed grasslands, an increase in ruminant stocking density causes a shift from grazed to mown grassland (globally and regionally, except in western Europe; Table 3 and Fig. 5), because mown grassland provides more grass biomass than grazed grassland per unit of area (Fig. S3).
A pasture-hyde is consistent with country-specific pasture area censuses and thus may be suitable for reconstructing land cover, but it does not provide information about management intensity. A managed-gm and its split between mown, grazed, and unmanaged fractions provide specific global distributions of pasture management intensity and its historical changes. However, there are several limitations, which may cause uncertainties in our maps of management intensity: (1) the grass fraction in ruminant diet has likely been changing during the last century while, due to a lack of information, we assumed that it was static in each region up to the year 2000; (2) technical developments (such as ruminant breeding) are not considered but may affect the feeding efficiency (meat/milk production per amount of feed) and thus feedback on the grass-biomass requirement; (3) the spatial distribution of ruminants was kept constant in our estimate, whereas it could have changed, depending on geographic changes in human population distribution; and (4) the results depend on the accuracy of NPP modelling in ORCHIDEE-GM. Despite these limitations, the maps of grassland management intensity provide new information for drawing up global estimates of management impact on biomass production and yields (Campioli et al., 2015) and for global vegetation models like ORCHIDEE-GM to enable simulations of carbon stocks and GHG budgets beyond simple tuning of grassland productivities (e.g. like in LPJmL; Bondeau et al., 2007) to account for management. These maps can also be tested in other vegetation models, or the same algorithm can be implemented in other models to give the management intensity consistent with simulated NPP.

Causes of regional grass-biomass production deficits
Grass-biomass production is constrained by the gridded biomass consumption for the year 2000 (Herrero et al., 2013). In some grid cells, the gridded biomass consumption by year 2000 cannot be fulfilled by the potential grass production simulated by ORCHIDEE-GM v3.1 (Fig. 2d). These modelled grass-biomass production deficits could be due to several reasons.
-Land-cover maps used as input to ORCHIDEE-GM v3.1 do not represent grasslands well in the mixed and landless systems and grasslands providing occasional feed to ruminant (e.g. roadside, forest understory grazing land, and small patches). This failing could cause the model to miss a significant part of grass productivity in this study. For example, the largest modelled grass-biomass production deficit is found in India because the simulated grassland productivity is far from agreeing with the grass-biomass-use data. In this country, occasional feed may constitute an important fraction of ruminant diet (30 or 50 % in mixed and landless or pastoral systems of South Asia from Bouwman et al., 2005), which is not represented by the land-cover maps used as input to ORCHIDEE-GM v3.1 and thus is not modelled.
-In arid regions such as Pakistan, Sudan, Iran, Egypt, and northwestern China, grass can grow in places where the water table is near to the surface and groundwater resources are available (e.g. oases, riparian zones, lakes). However, ORCHIDEE-GM v3.1 is driven by gridded climate data and does not taken into account local topography-dependent water resources such as rivers and lakes and thus is not being able to simulate local grass growing areas in arid regions.
-Grassland irrigation, though it is not as common as in cropland, is applied in arid regions such as Saudi Arabia but is not considered by ORCHIDEE-GM v3.1.
-In some semi-arid open rangeland, ruminants may walk long distances to acquire enough grass. For example, in semi-arid sub-Saharan Africa, Uzbekistan, and central Australia, animals usually keep moving in order to search for grass. This displacement of grazing animals from grass sources is not considered in the model.
-The grass fraction in ruminant diet is defined per region according to specific production systems. However, the grass fraction can differ within a region depending on local fodder crop production and grassland use. For example, the large numbers of ruminants in eastern China are mostly fed by grain and stovers (fibrous crop residues) instead of grass, because little grassland exists in that region.

Model performance: comparison of modelled and observed grassland productivity
In Sect. 3.3, the spatial patterns of NPP model or GPP model were compared with observations (NPP obs or MODIS-GPP). ORCHIDEE-GM v3.1 captured well the spatial pattern of grassland productivity, with (i) high r spatial between GPP model and MODIS-GPP (Sect. 3.3.2) and (ii) NPP model extracted from global simulation showing significant correlation with site-level NPP observation from 129 sites all over the world (Sect. 3.3.1). However, GPP model is higher than MODIS-GPP in all latitude bands (Fig. 8). It should be kept in mind that MODIS-GPP had a calculated 18 % uncertainty due to climate forcing (Zhao et al., 2006). Besides, a low bias of MODIS-GPP for grasslands has been reported in a tallgrass prairie in the United States (Turner et al., 2006) and in an alpine meadow on the Tibetan Plateau (Zhang et al., 2008) when compared to the GPP from flux-tower measurements. The underestimate of MODIS-GPP is mostly related to the low value of the maximum light-use efficiency parameters used in the MODIS-GPP algorithm (Turner et al., 2006;Zhang et al., 2008). The relatively low r value between NPP model and site-level NPP obs (r = 0.35, p<0.01; Sect. 3.3.1) could be related to the fact that local climate, soil properties, and topographic features are not considered in the model. For example, the r between the site-level climate and that from the CRU+NCEP climate forcing data (0.5 • × 0.5 • resolution) is 0.96 for annual mean temperature but only 0.86 for annual total precipitation and 0.86 for solar radiation. The relatively low correlation for annual total precipitation may cause inaccuracy in the model simulations of productivity, because water availability could be a major factor limiting grass growth (e.g. in temperate regions; Le Houerou et al., 1988;Silvertown et al., 1994;Briggs and Knapp 1995;Knapp et al., 2001;Nippert et al., 2006;Harpole et al., 2007). Further, a similar mean belowground NPP and an overestimation of mean aboveground NPP by ORCHIDEE-GM v3.1 is found in Sect. 3.3.1, which suggests that (1) the model tends to overestimate aboveground NPP possibly due to overestimation of GPP (compared to MODIS-GPP) and (2) the model tends to overestimate the ratio of aboveground and belowground biomass allocation (R above/below ) compared to observation. This overestimation could be the result of nitrogen limitation on the carbon allocation scheme for grassland. For example, a large nitrogen supply has been observed to increase R above/below (Aerts et al., 1991;Cotrufo and Gorissen, 1997), while nitrogen limitation might cause it to decrease. However, nitrogen limitation in grassland is not accounted for in ORCHIDEE-GM v3.1, which possibly leads to the model's overestimation of R above/below . The model could be improved by incorporating the full nitrogen cycle.
For the seasonal cycle, we compared modelled GPP seasonality to both MODIS-GPP and GOME-2 SIF data. ORCHIDEE-GM v3.1 captures the seasonal variation of productivity in most regions where grassland is the dominant ecosystem (coverage > 50 %), as shown by the high r seasonal between GPP model and MODIS-GPP (Fig. S6a) or SIF data (Fig. S6b). However, the model does not capture the seasonal amplitude of grassland productivity in some arid/semi-arid regions (e.g. southwestern United States and central Australia; Fig. S6a and b). In arid/semi-arid regions, grass productivity is triggered by discrete precipitation events and depends on the timing and magnitude of these pulses (Sala et al., 1982;Schwinning and Sala, 2004;Huxman et al., 2004).
These precipitation pulses are infrequent, discrete, and not represented in a global climate re-analysis dataset such as CRU+NCEP used in our simulation. In particular, NCEP, like all climate models tends to produce "general circulation model drizzle" (Berg et al., 2010), i.e. too many frequent small rainfall events. This forcing uncertainty could be a major obstacle for our model to capture the seasonality of productivity in these regions. In dry grasslands, the dominant species could change during the season, but the resultant changes in SLA and Vc max 25 by different dominant species cannot be reflected in ORCHIDEE-GM v3.1. This withinseason variability could be another reason for the modeldata discrepancy in arid/semi-arid grassland seasonality. For the savanna of sub-Saharan Africa, eastern Africa, and South America (Fig. S6), the relatively low r seasonal could be a result of the fact that the frequent fires are not simulated in the current version of the model used here.
ORCHIDEE-GM v3.1 captures the IAV of grassland GPP at global scale and in many regions of the world (40 % of global grassland area) compared to the MODIS-GPP. One exception where IAV is not in phase with MODIS-GPP is sub-Saharan Africa (Fig. 9). Possible causes of this discrepancy are (1) the frequent fires which affect the IAV of GPP, which are not simulated in this study; (2) model biases in the IAV of soil moisture, which could affect the model performances for the productivity of semi-arid Africa, given its two-layer bucket hydrology; (3) the problems with MODIS-GPP dry areas, which may degrade the model-data agreement. The cold Qinghai-Tibet plateau and boreal tundra are the other regions where the model does not capture the GPP IAV (Fig. 9). The low model-data agreement in IAV could be due to shortcomings in (1) the specific characteristics, functioning traits, and nutrient availability of the tundra/alpine-grassland ecosystem that are not well parameterized or accounted for in our model (e.g. Tan et al., 2010, for Qinghai-Tibet plateau) and (2) the snow scheme. The timing of snowmelt will impact the grass phenology, while early spring soil moisture impacted by snow water storage may affect the grassland productivity. The single-bucket snowpack scheme (Chalita and Le Treut, 1994) in the current version of ORCHIDEE-GM may not represent the snow processes sufficiently accurately. The mechanistic intermediatecomplexity snow scheme (ISBA-ES; Boone and Etchevers, 2001) implemented into ORCHIDEE-ES  may improve the model performance in simulating grassland productivity.

Concluding remarks
In this study, we have derived the global gridded maps of grassland management intensity, including the minimum area of managed grassland with fraction of mown/grazed part, the grazing-ruminant stocking density, and the density of the wild animal population at a resolution of 0.5 • by 0.5 • . The management intensity maps are built based on the assumption that grass-biomass production from managed grassland (simulated by ORCHIDEE-GM v3.1) in each grid cell is just enough to satisfy the grass-biomass requirement by ruminants in the same grid (data derived from Herrero et al., 2013). Furthermore, the maps are extended to cover the period 1901-2012, taking into account both the changes in grass-biomass requirement and supply. The evolution in grass-biomass requirement is determined by the ME-based ruminant numbers calculated in this study, while the changes in grass-biomass supply are simulated by ORCHIDEE-GM v3.1, considering variable drivers such as climate, CO 2 concentration, and N fertilization. Despite the multiple sources of uncertainty, these maps, to our knowledge for the first time, provide global, time-dependent information on grassland management intensity. Global vegetation models such as ORCHIDEE-GM, containing an explicit representation of grassland management, are now able to use these maps to make a more accurate estimate of global carbon and GHG budgets.
The gridded grassland management intensity maps are model dependent because they depend on NPP. Thus in this study we also give a specific attention to the evaluation of modelled productivity against both a new set of sitelevel NPP measurements and global satellite-based products (MODIS-GPP and GOME2-SIF). Generally, ORCHIDEE-GM v3.1 captures the spatial pattern, seasonal cycle, and IAV of grassland productivity at global scale, except in regions with either arid or cold climates (tundra) and high-altitude mountains/plateaus. Because the major purpose of a global vegetation model like ORCHIDEE-GM is to simulate carbon, water, and energy fluxes at a large scale, it uses a limited number of plant functional types and generic equations. The model is not expected to accurately capture productivity variations everywhere. Thus we conclude that its current version, ORCHIDEE-GM v3.1, is suitable to simulate global grassland productivity.
The Supplement related to this article is available online at doi:10.5194/bg-13-3757-2016-supplement.