# Modeled Changes in Potential Grassland Productivity and in Grass-Fed Ruminant Livestock Density in Europe over 1961–2010

Jinfeng Chang, Nicolas N. Viovy, Nicolas N. Vuichard, Philippe Ciais, Matteo

Campioli, Katja Klumpp, Raphaël Martin, Adrian Leip, Jean-François

Soussana

► To cite this version:

Jinfeng Chang, Nicolas N. Viovy, Nicolas N. Vuichard, Philippe Ciais, Matteo Campioli, et al. Modeled Changes in Potential Grassland Productivity and in Grass-Fed Ruminant Livestock Density in Europe over 1961–2010. PLoS ONE, 2015, 10 (5), e0127554. doi:10.1371/journal.pone.0127554

HAL Id: hal-01806202

https://hal.science/hal-01806202

Submitted on 27 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

# Modeled Changes in Potential Grassland Productivity and in Grass-Fed Ruminant Livestock Density in Europe over 1961–2010

Jinfeng Chang¹*, Nicolas Viovy¹, Nicolas Vuichard¹, Philippe Ciais¹, Matteo Campioli², Katja Klumpp³, Raphaël Martin³, Adrian Leip⁴, Jean-François Soussana³

1. Laboratoire des Sciences du Climat et de l'Environnement, UMR8212, CEA-CNRS-UVSQ, Gif-sur-Yvette, France,
2. Research Group of Plant and Vegetation Ecology, Department of Biology, University of Antwerp, Wilrijk, Belgium,
3. INRA, Grassland Ecosystem Research Unit, UREP, Clermont-Ferrand, France,
4. European Commission, Joint Research Centre (JRC), Institute for Environment and Sustainability, Ispra (VA), Italy

* jinfeng.chang@lsce.ipsl.fr

## Abstract

About 25% of European livestock intake is based on permanent and sown grasslands. To fulfill rising demand for animal products, an intensification of livestock production may lead to an increased consumption of crop and compound feeds. In order to preserve an economically and environmentally sustainable agriculture, a more forage based livestock alimentation may be an advantage. However, besides management, grassland productivity is highly vulnerable to climate (i.e., temperature, precipitation, CO₂ concentration), and spatial information about European grassland productivity in response to climate change is scarce. The process-based vegetation model ORCHIDEE-GM, containing an explicit representation of grassland management (i.e., herbage mowing and grazing), is used here to estimate changes in potential productivity and potential grass-fed ruminant livestock density across European grasslands over the period 1961–2010. Here “potential grass-fed ruminant livestock density” denotes the maximum density of livestock that can be supported by grassland productivity in each 25 km × 25 km grid cell. In reality, livestock density could be higher than potential (e.g., if additional feed is supplied to animals) or lower (e.g., in response to economic factors, pedo-climatic and biotic conditions ignored by the model, or policy decisions that can for instance reduce livestock numbers). When compared to agricultural statistics (Eurostat and FAOstat), ORCHIDEE-GM gave a good reproduction of the regional gradients of annual grassland productivity and ruminant livestock density. The model however tends to systematically overestimate the absolute values of productivity in most regions, suggesting that most grid cells remain below their potential grassland productivity due to possible nutrient and biotic limitations on plant growth. When ORCHIDEE-GM was run for the period 1961–2010 with variable climate and rising CO₂, an increase of potential annual production (over 3%) per decade was found: 97% of this increase was attributed to the rise in CO₂, -3% to climate trends and 15% to trends in nitrogen fertilization and deposition. When compared with statistical data, ORCHIDEE-GM captures well the observed phase of climate-driven interannual variability in grassland production well, whereas the

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

### OPEN ACCESS

**Citation:** Chang J, Viovy N, Vuichard N, Ciais P, Campioli M, Klumpp K, et al. (2015) Modeled Changes in Potential Grassland Productivity and in Grass-Fed Ruminant Livestock Density in Europe over 1961–2010. PLoS ONE 10(5): e0127554. doi:10.1371/journal.pone.0127554

**Received:** November 24, 2014

**Accepted:** April 16, 2015

**Published:** May 27, 2015

**Copyright:** © 2015 Chang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This work was supported by the AnimalChange project (grant no. 266018). The authors' acknowledge the PhD funding for JC by the French l'Institut national de la recherche agronomique (INRA) and the French Commissariat à l'énergie atomique (CEA).

**Competing Interests:** The authors have declared that no competing interests exist.
magnitude of the interannual variability in modeled productivity is larger than the statistical data. Regional grass-fed livestock numbers can be reproduced by ORCHIDEE-GM based on its simple assumptions and parameterization about productivity being the only limiting factor to define the sustainable number of animals per unit area. Causes for regional model-data misfits are discussed, including uncertainties in farming practices (e.g., nitrogen fertilizer application, and mowing and grazing intensity) and in ruminant diet composition, as well as uncertainties in the statistical data and in model parameter values.

## Introduction

Consumption of animal-product is increasing in the global food diet, and this trend is projected to continue with large increases in demand for both dairy products and ruminant meat—mainly beef but also lamb and goat [1]. In turn, this rising demand for animal products will impact the demand for animal feed (grown as crops and pasture) and is expected to lead to an intensification of livestock production. Currently, the proportion of pasture-based feed is decreasing, mirrored by increased use of crop-based feeds [2]. Yet grass forage still accounts for around 70% of the global dry matter intake of ruminants (world mean in 1993 [3]). For Europe, Leip et al. [4] estimate that herbage intake provides about 25% of livestock's protein in the feed, while about 60% are from crops and compound feeds. Of the approximately 13 billion hectares of ice-free land on Earth, 3.4 billion hectares (about 25%) are grasslands, an area 2.3 times larger than that used for arable crops ([5], data of 2010).

Improving grassland productivity and the management of livestock systems is especially important if we are to create economically, socially and environmentally sustainable agriculture [6]. The European livestock sector produces 15 and 21% of the global totals of meat, and of milk and dairy products, respectively, amounting to an economic value of $185 billion ([5], data of 2010). Ruminant livestock is largely bred using grassland as a key resource. Grassland provides fresh herbage for grazing and is also mown to create hay and silage that are predominantly used during the winter period. Most European grassland is infrequently (permanent grasslands) or never (rough pastures) re-sown, but these semi-natural grasslands, which are characterized by a moderate to high plant species diversity, have been increasingly fertilized with both inorganic and organic nitrogen over recent decades [7]. Apart from these so-called rough pastures, which are usually located in mountainous areas with steep slopes and limited accessibility, fulfilling the forage requirement of ruminant livestock requires permanent and sown grasslands in Europe to be used (at farm scale) for both grazing and mowing.

Grassland productivity is not only affected by management, but also responds to rising atmospheric CO₂ concentration [8] and climate change. Elevated CO₂ concentration has the dual effect of increasing leaf photosynthesis and reducing stomatal conductivity. These effects tend to increase water-use efficiency [9] and often reduce the consumption of soil moisture [8], but decreased conductance could be offset by increased leaf area index (LAI) in the plant transpired water flux to the atmosphere. Elevated CO₂ effects on photosynthesis are larger with the C3 species that dominate temperate pastures than with C4 species [10, 11]. Temperate C3 grasses show a typical 30% increase in canopy photosynthesis under doubled atmospheric CO₂ concentration [12, 13]. Water availability plays a major role in the response of grasslands to climate change, with marked declines in productivity under increased soil water deficits—although there are differences in species response [14]. In the future, elevated atmospheric CO₂ concentration is expected to reduce the sensitivity of grassland productivity to low
precipitation in grassland ecosystems [12, 15]. However, trends and variability in temperature and precipitation, as well as nitrogen limitations, will all interact with the effects of elevated CO₂ to determine actual changes in grassland productivity [8, 16].

Despite the importance of grasslands in sustaining ruminant livestock production, spatially explicit information about European grassland productivity and its response to climate change is scarce. A few studies have attempted to estimate the potential productivity for individual regions [17–19]. Smit et al. [20] presented a European diagnostic map of grassland productivity derived from census statistics across regions in Europe, integrated with MODIS satellite NDVI observations. However, the mechanisms that control productivity were not investigated in that study. Climate change impacts, including the altered frequency and/or severity of extreme events, and nitrogen deposition impacts on grasslands are better understood at the site or small regional level than at larger scales (e.g., [21, 22]; reviewed by Tubiello et al. [23]).

Projections of climate change impacts on grassland productivity are mostly based on local modeling studies that require many local variables, which limits up-scaling to regional or continental scale (e.g., [24]). Model-based impact assessment studies use relationships to simulate grassland productivity that have usually been derived from field experiments (e.g., [25–27]). However, the applicability of these relationships to larger scales and future conditions is highly uncertain [20]. Process-based vegetation models with equations representing biogeochemical and biophysical mechanisms and the combined response to multiple drivers have the advantage of being applicable from plot to continental scale. Process-based models can also be coupled with regional or global climate models to quantify land-atmosphere feedbacks. This type of mechanistic model is increasingly used for global impact studies focusing on agricultural productivity and terrestrial carbon fluxes (e.g. ISI-MIP, http://www.isi-mip.org).

In this study, a process-based land surface model specifically calibrated and parameterized to simulate the carbon cycle of European managed grasslands ORCHIDEE-GM [28–29] is used to address two questions:

1. What is the geographic distribution of potential (maximal) grassland productivity and potential grass-fed livestock density of managed grasslands over the EU27?

2. How did potential grassland productivity evolve during the past 50 years in response to climate change, atmospheric CO₂ concentration and nitrogen fertilizer supply?

# Material and Methods

## Model description

ORCHIDEE is a process-based ecosystem model built for simulating global carbon, water and energy fluxes [30]. ORCHIDEE-GM is a version of this model that has been developed to explicitly represent grassland management such as mowing and livestock density. Grassland management was developed by implementing the management module from PaSim [28, 31–33], a grassland model developed initially for site applications, into ORCHIDEE. The equations that differentiate ORCHIDEE-GM from the standard version of ORCHIDEE are described by Chang et al. [29], together with a comparison of both model versions at 11 sites equipped with eddy-covariance and biometric measurements. The model assumes the simple principle that farmers will produce the maximum number of animals that can be fed locally by above ground productivity, optimally combining grazing and mowing (local forage production) [28] as well as nitrogen fertilizer application, although the model does not explicitly simulate the nitrogen cycle [29].
## Effects of nitrogen addition on productivity in ORCHIDEE-GM

Application of nitrogen to pasture leads to higher nitrogen concentration in plant leaves and enhanced leaf photosynthetic capacity [34, 35]. In a meta-analysis gathering many nitrogen addition studies, Xia & Wan [36] reported a 46.5% increase in grass biomass under nitrogen addition, close to another estimate (53% increase of aboveground net primary production (ANPP); [37]). In addition, the effect of N addition to productivity depends on some other factors such as original soil nitrogen availability, nitrogen-fixing plants in grassland community, and management intensity (e.g., herbivores can modify soil N cycling; see a review by Bardgett and Wardle [38]). But the data about these factors are not available over European grassland. Ideally, the relationships among photosynthesis (productivity), leaf nitrogen content, and soil nitrogen availability (can be altered by nitrogen addition) can be modeled with coupled nitrogen/carbon cycling (e.g., O-CN [39]). However, the nitrogen cycle, and thus the interactions between nitrogen and carbon cycles, is not included in ORCHIDEE-GM.

In this study, we tried to apply a simplified function to reproduce empirically the effect of nitrogen addition on grass growth, based on the following fact: i) generally, increased soil nitrogen availability for increased photosynthesis is the result of nitrogen addition, ii) the experimental phenomena [40] showed that grassland productivity was increased along a nitrogen addition gradient, but with a limit after which the productivity will not increase. A variable, representing nitrogen response $N_{add}$ was added to ORCHIDEE-GM. $N_{add}$ has a directly additive impact on two photosynthetic parameters, $Vc_{max}opt$ (the maximum rate of Rubisco carboxylase activity), and $J_{max}opt$ (the maximum rate of photosynthetic electron transport), as given by:

$$Vc_{max}opt = (1 + N_{add}) \times Vc_{max}opt^* \quad (1)$$

$$J_{max}opt = (1 + N_{add}) \times J_{max}opt^* \quad (2)$$

where $Vc_{max}opt^*$ and $J_{max}opt^*$ are the parameter values determining the photosynthetic rates and indirectly, plant productivity and growth in ORCHIDEE. $N_{add}$ is assumed to be a function of the amount of nitrogen added ($N_{amount}$ expressed as kg N ha$^{-1}$ yr$^{-1}$) with a saturating curve limited by $N_{addmax}$ representing the maximum asymptotic response to nitrogen addition [40]. A saturating exponential function was chosen to describe the N addition effect:

$$N_{add} = N_{addmax} \times (1 - a^{N_{amount}/30}) \quad (3)$$

where a is a constant (set to 0.75 here) and $N_{addmax}$ is set to 0.6. With those parameter values, leaf photosynthesis is increased by about 60% under very high nitrogen supply, which is in agreement with the positive effects of nitrogen addition on grassland productivity [36, 37].

## Biological potential productivity and grass-fed livestock density

ORCHIDEE-GM simulates two options, through which managed grassland provides forage for livestock: grazing and mowing. Under grazing conditions, thresholds of shoot biomass are set for starting, stopping and resuming grazing [28], while under mowing, the frequency and magnitude of forage harvests in each grid cell is a function of grown biomass [28], and the model carbon dynamics equations then calculate the potential productivity. Here, we define grassland productivity as the annual production of forage from cut grassland. Simulated productivity can be validated against the variable names yields in the European Eurostat statistics at the scale of administrative units (Eurostat Agriculture, forestry and fisheries; referred as Eurostat-AFF [41] hereafter). Spatial statistical information on stocking rates (grazing-animal density) and on the ratio between cut and grazed areas in each grid cell is not available at European scale. Thus a
set of rules describing an idealized self-sufficient herbage-based ruminant livestock farm [28] has been introduced into ORCHIDEE-GM. These rules are based upon three assumptions, applied for each grid cell, 1) livestock is only fed by herbage (i.e. arable crop-feed products are not considered), 2) grazing occurs throughout the vegetation period, and that mown biomass must satisfy the feed requirements of the animals outside the growing season, and that no lateral herbage (e.g. hay) goes in or out of the grid cell, and 3) the use of grassland production is maximized by determining a combination of grazing and mowing, that maximizes the number of animals. Our model thus calculates how much livestock may potentially be fed by grass biomass in each grid cell. Accordingly, we call the results of the model built upon these rules, potential or optimal productivity and animal stocking rate.

Under these assumptions, the optimal animal stocking rate, $S_{opt}$ (number of livestock units (LSU) per hectare) and the optimal proportion of grazed versus cut grasslands, $F_{opt}$ (within [0,1]) are calculated for each grid cell. The annual herbage production, $Y_{cut}$, (kg of dry matter (DM) per hectare) of cut grasslands that occupy a fraction $(1-F_{opt})$ of a grid cell should be equal to the herbage dry matter required by herbivores under cover (Eq 4), while the production under grazing ($Y_{graze}$) needs to meet dry matter requirements during the vegetation period (Eq 5). Due to the impact of livestock on grass growth through trampling, defoliation etc., and because grassland cannot be continuously grazed during the vegetation period, the effective yield on grazed grassland depends on the stocking rate and on the environmental conditions of the grid cell (Eq 6). Given these interdependencies, $S_{opt}$ and $F_{opt}$ are calculated using the optimization algorithm of Vuichard et al. [28]. The ‘recovery’ time required under grazing is obtained in the model using threshold, which determine when grazing stops (dry biomass remaining lower than 300 kg DM ha⁻¹), or when grazing can start again (dry biomass recovered to a value above 300 kg DM ha⁻¹ for at least 15 days).

$$IC \times T_{farm} \times S_{opt} \times F_{opt} = Y_{cut} \times (1 - F_{opt}) \quad (4)$$

$$IC \times (1 - T_{farm}) \times S_{opt} \times F_{opt} = Y_{graze} \times F_{opt} \quad (5)$$

$$Y_{graze} = f\{\text{grid cell}, S_{opt}\} \quad (6)$$

where IC is the daily intake capacity (with a mean value of 13 kg DM LSU⁻¹ day⁻¹, assuming 1 LSU is the grazing equivalent of one adult dairy cow producing 3000 kg of milk annually, without additional concentrated foodstuffs [42]; using typical digestible energy expressed as a percentage of gross energy (DE%: 45%- 55%) suggested by IPCC guidelines [43], IC ranges from 12–14 kg DM LSU⁻¹ day⁻¹ calculated based on Eq 10.18b in IPCC guidelines [43]) of animals during grazing and under cover, $T_{farm}$ is the number of non-growing season days during which herbivores need to be fed with cut herbage. The potential animal density ($D_{opt}$) is given by:

$$D_{opt} = S_{opt} \times F_{opt} \quad (7)$$

$D_{opt}$ is calculated in the iterative algorithm by increasing the input animal stocking rate ($S_{opt}$) until convergence is reached and the stocking rate reaches its potential value and cannot increase further. $D_{opt}$ must be interpreted as a potential herbage-only limited (i.e., grass-fed) livestock density of a given model grid cell (similar to the term “livestock carrying capacity” [44]). When $D_{opt}$ is reached, the grassland herbage production is fully used by livestock and the herbage intake capacity of the livestock is reached. The model assumes that the grassland herbage is highly digestible, i.e., dominated by leafy herbage with a low proportion of stem and dead material, which is usually not the case in rough pastures [7]. Furthermore, other limiting factors, such as low soil pH, low soil nitrogen mineralization and low soil phosphorus, the density of
unpalatable plants that animals refuse to eat, soil degradation and slope are not taken into account. Therefore, the actual grass-fed livestock density tends to be lower than $D_{opt}$ because those factors reduce grassland productivity. On the other hand, the use of additional arable crop-feed diet supplements justified by economic considerations (e.g., estimated in [3]), for instance in the production of high economic value dairy products, may locally increase the actual livestock density above $D_{opt}$.

## Modeling adaptive management in response to climate variability and CO₂

The algorithm developed by Vuichard et al. [32] is designed to define the optimal number of livestock ($D_{opt}$) in equilibrium with a given climate condition. We then improved the algorithm to be able to simulate the progressive adaptation of livestock density to climate change. New rules were defined and incorporated into ORCHIDEE-GM to account for how grassland management might change in response to a climate-driven change in productivity. We start from animal stocking rate ($S_{opt}$) and proportion of grazed grasslands ($F_{opt}$) optimized for a productivity in equilibrium with a given climate condition. In the case productivity changes (e.g., because of climate change), the profit-driven farmers will tend to fully use the herbage production from the farm, and in practically, take adaptive actions to restore the biological potential grass-fed livestock density by changing livestock numbers (Fig 1). Specifically, at the end of a given year, $i$, we model that the farm manager in each grid cell takes the decision to add (in case of forage surplus) or to remove (in case of forage deficit) animals to or from the grid cell for the next year, thus changing the previous optimal animal density ($D_{opt,i}$) by a step $\Delta D_{opt,i}$ as function of the cut forage surplus or deficit ($\Delta M_i$). The new optimal animal density for the next year $i+1$ ($D_{opt,i+1}$) is given by:

$$D_{opt,i+1} = D_{opt,i} + \Delta D_{opt,i} \qquad (8)$$

<figure><img src="image_4.png" /><figcaption>Fig 1. Strategy for modeling adaptive management in ORCHIDEE-GM. See text for symbol definitions.</figcaption></figure>

doi:10.1371/journal.pone.0127554.g001
where $D_{opt,i+1}$ is adapted in the model by changing $F_{opt}$ (in Eq 7), and the step density change ($\Delta D_{opt,i}$) is given by:

$$ \Delta D_{opt,i} = \Delta D_{opt} \times \alpha = \Delta M_i / IC / T_{year} \times \alpha \quad (9) $$

where IC is the daily intake capacity, $T_{year} = 365$ days. $\Delta D_{opt,i}$ is the maximum animal density change resulting from the forage surplus (or deficit) of the previous year ($\Delta M_i$), and $\alpha$ is the fraction of $\Delta D_{opt}$ that is realized by the farmer. Given the highly variable productivity caused by the year-to-year variability of climate, the farmers may face a very high risk by adjusting each year their livestock numbers to productivity variations (the situation of $\alpha$ is 100%). For example, in a situation of an extreme drought year $i+1$ follows a very productive year $i$, farmers will suffer a heavy loss in year $i+1$ due to their acquisition of too many animal based on the productivity of year $i$. Therefore, the farmers were assumed to be able to search a balanced response to productivity change between the profit (e.g., get more profit from productivity increase) and the risk (year-to-year variability of productivity). We assumed they tend to choose a reasonable “moderate risk”, which is the response to the productivity change during a relatively long period (e.g., previous 5-years’ mean productivity, corresponding to the $\alpha$ of 20% set in this study). $\Delta M_i$ is given by:

$$ \Delta M_i = X'_{i} - X_i \quad (10) $$

$$ X'_{i} = Y_{cut_i} \times (1 - F_{opt_i}) \quad (11) $$

$$ X_i = IC \times T_{farm_i} \times S_{opt_i} \times F_{opt_i} \quad (12) $$

where, in each grid cell, $X'_{i}$ is the annual total forage grass production from cut grasslands for the year $i$, calculated from Eq 11, and $X_i$ is the herbage dry matter required by herbivores under cover during the non-growing season for the year $i$, calculated from Eq 12.

## Simulation set-up

ORCHIDEE-GM is integrated on a grid over Europe using the harmonized climate forcing data from the ERA-WATCH reanalysis for the period 1901–2010 and at a spatial resolution of 25' by 25' [45]. Mean and standard deviation of the ERA-Interim time series [46] were adjusted according to the WATCH time series [47, 48] by using the overlapping period 1989–2001. The harmonized dataset was spatially downscaled to 25' by overlapping CRU CL2.0 [49] monthly means to the spatial anomaly of the harmonized datasets for each single climatic variable. An altitude-based correction was applied for downscaling surface pressure according to a digital elevation map from CRU CL2.0. This resolution (25' by 25') is sufficient to represent regional meteorological regimes accurately in low lying regions, but not in mountainous areas.

Gridded mineral fertilizer and manure nitrogen application rate for European grasslands in the European Union (EU27) was estimated by CAPRI model (see [50, 51]) based on combined information from official and harmonized data sources such as Eurostat, FAOstat and OECD, and spatially dis-aggregated using the methodology described in Leip et al. [52]. The data are estimated at a spatial resolution of clusters of 1 km by 1 km and were re-aggregated here to a spatial resolution of 25' by 25' (S1 Fig). For French regions, we use data from the French national statistics [53]. A set of rules was used to rebuild the temporal evolution of gridded nitrogen fertilization from 1901 to 2010: 1) organic fertilizer is assumed to have remained constant over time; 2), the application rate of mineral fertilizer evolved with time following the total mineral nitrogen fertilizer consumption of European Union [54]; 3) mineral fertilizers are set to be applied since 1951, and application rates linearly increased from 0 to the observed level of 1961 during the period 1951–1960. Besides nitrogen fertilizer application, nitrogen deposition
<figure><img src="image_5.png" /><figcaption>**Fig 2.** Illustration of the simulation protocol, forcing data and initial state for various simulations.</figcaption></figure>

doi:10.1371/journal.pone.0127554.g002

was considered as nitrogen addition as well. Gridded nitrogen deposition rates for Europe was originated from European Monitoring & Evaluation Programme (EMEP) dataset, and was a product of EU-PF7 project (GHG-Europe; data are available at http://gaia.agraria.unitus.it/ghg-europe/data/others-data) where the decadal means were linear interpolated to annual values. It is noteworthy that the N fixation by legumes was not accounted for in this study.

A series of simulations was carried out, and the illustration of the simulation protocol is shown in Fig 2. ORCHIDEE-GM was first run for a spin-up (simulation E1) without management using the first 10 years of the period 1901–1910 recycled in a loop, and atmospheric CO₂ concentration for 1900 (296 ppm) until all carbon pools reached equilibrium (long term Net Ecosystem Exchange, NEE = 0 at each grid point). This first spin-up usually takes 10,000 years. Starting from the end of this first spin-up, two separate transient simulations were performed. In the first one (simulation E2), optimal animal stocking rates ($S_{opt}$) and fractions of grazed grassland ($F_{opt}$) for the period 1901–1910 were defined by running ORCHIDEE-GM with its optimization algorithm that maximizes stocking rates in each grid cell (see above). In the simulation E3, ORCHIDEE-GM was run during the historical period (1911–1960) with increasing atmospheric CO₂, variable climate, and with the new climate-adaptive management change
algorithm. This transient E3 simulation started with the reference distributions of $S_{opt}$ and $F_{opt}$ over Europe (obtained in simulation E2), as well as with soil carbon pools for the year 1910 (end of the spin-up simulation). Starting from the end of simulation E3, a transient simulation E4 was carried out for the period 1961–2010 with increasing atmospheric CO₂, variable climate, and with the adaptive management change algorithm.

Three further simulations (E5, E6 and E7) investigated the relative contribution of atmospheric CO₂, climate change and nitrogen fertilization trends on the estimated trend in productivity. The simulations E5, E6 and E7 are the same as simulation E4, but with atmospheric CO₂ concentration fixed at the level of 1961 (E5), using the first five years of climatology data (1961–1965) recycled (E6), or with nitrogen fertilization and nitrogen deposition fixed to the level of 1961 (E7), respectively. The difference in productivity trend between simulations E4 and E5 reflects the effects of increased CO₂. The effects of climate variation and nitrogen fertilization were derived as the difference between simulations E4 and E6 and between simulations E4 and E7, respectively.

## Regional census statistics (grassland area, livestock numbers)

The European statistical database, Eurostat [41], divides the European Union (EU) into regions at three different levels abbreviated as NUTS (Nomenclature des Unités Territoriales Statisques [55]). Above NUTS-1 (referred as *NUTS-0* hereafter) is the country level of the Member State. Eurostat [41] provides data on grassland area and livestock numbers for NUTS regions [55]. A full scope Farm Structure Survey (FSS; Council Regulation (EEC) No 571/88, and Regulation 1166/2008) is carried out every 10 years since 1990, with sample surveys every 2 or 3 years. In this study, we compiled data from 30 countries (EU28 plus Norway, and Switzerland) and 279 regions (NUTS-2 level except for Germany, where only NUTS1-level data are available) of Europe. Eurostat statistics divide grassland into classes: *temporary grassland* (productive grassland less than five years after sowing), *permanent grassland* (productive grassland, five years or more after sowing), and *rough grazing* (low productivity pasture, extensively grazed, often in fragile areas). Ruminant livestock of various species and ages are aggregated into livestock units (LSU), by applying specific coefficients initially established on the basis of the nutritional or feed requirement for each type of animal. Livestock unit coefficients are derived from Eurostat [42] (e.g., 0.1, 0.8 and 1 LSU per head for sheep and goats, and heifers and dairy cows, respectively).

## Data for evaluating the simulated productivity of European grasslands

Smit et al. [20] constructed a map of Europe showing the spatial distribution of grassland productivity by integrating census statistics, literature, and expert judgment using the NUTS classification [55]. The biological potential of grassland productivity from ORCHIDEE-GM (on a spatial resolution of 25') was aggregated to the NUTS-2 level weighted by the corresponding grassland area in each grid cell (from CORINE Land Cover map, CLC2000 [56]). We then compared the spatial pattern of modeled potential productivity with the map of Smit et al. [20] averaged over the period 1995–2004.

Net primary productivity (NPP), including above- and below-ground plant organs, represents the net flux of carbon from the atmosphere into plants per unit time (one year in this study). NPP is a crucial variable in vegetation models. However, high quality measurements of grassland NPP are sparse due to the difficulty of measuring some NPP components such as fine-root production [57]. A new version of the Luyssaert et al. [58] database comprising non-forest biomes (Campioli, unpublished) was obtained. NPP data from seven European
temperate grassland sites of this database were compared with NPP of cut grassland simulated with ORCHIDEE-GM at their corresponding locations on the grid.

The temporal evolution of grassland productivity is not available for most regions. However for some countries, the Eurostat statistical database provides long-term grassland yield data from 1973 onwards, which are directly calculated as harvested production per unit area [59]. NDVI has been shown to be positively correlated with temperate grassland productivity [60]. We therefore used the NOAA/AVHRR NDVI composites at a spatial resolution of 8-km and 15-day intervals produced by the Global Inventory Modeling and Mapping Studies (GIMMS) [61] from 1982 onwards. Summer NDVI (June, July and August, JJA) of grid cells dominated by grassland (cells with more than 50% of the area are attributed to grassland in CLC2000 land cover map) was extracted from GIMMS and then aggregated to NUTS-0 level for the countries with long-term grassland productivity data from Eurostat. The time series of grassland productivity statistics (aggregated over the NUTS regions) were compared with the modeled productivity of cut grassland (see above), and grassland NDVI was correlated with aggregated modeled productivity in the 25 km grid cells corresponding to those used for NDVI.

## Livestock distribution and its ratio to productivity

FAO [62] provides a 5' by 5' global livestock distribution maps for major animal species (cows, pigs, poultry, sheep, goats, and buffalo), which are consistent with regional statistics. Cattle, sheep and goat density is expressed in this study as the number of animals (head) per square kilometer of land suitable for livestock production [62]. A ruminant livestock density in LSU has¹, was calculated from these data by using head-to-LSU conversion factors of 0.8, 0.1, and 0.1 for cattle, sheep and goats, respectively. This LSU density distribution was then aggregated to the 25' by 25' grid used by ORCHIDEE-GM. We also aggregated the livestock density distribution from FAO to the NUTS-2 level, for comparison with grassland productivity statistics from Eurostat summarized by Smit et al. [20]. Then, the average ratio of FAO livestock density to grassland productivity from the Eurostat data was calculated in each NUTS region.

## Grass-fed livestock numbers estimated from statistics and model simulations

In Europe, ruminant livestock are not only fed on grass, but also receive crop-feed and crop by-products. Thus, the number (LSU) of grass-fed livestock ($N_{obs}$) in each region can be calculated separately as:

$$N_{obs} = N_{beef} \times F_{beef} + N_{dairy} \times F_{dairy} + N_{sheep} \times F_{sheep} + N_{goats} \times F_{goats} \quad (13)$$

where $N_{beef}$, $N_{dairy}$, $N_{sheep}$ and $N_{goats}$ are the actual numbers (LSU) of beef cattle, dairy cattle, sheep, and goats given by regional census statistics; $F_{beef}$, $F_{dairy}$, $F_{sheep}$ and $F_{goats}$ are the fractions of grass-feed in the diet of each type of animal, available from the study by Wirsenius [3] (differentiating Western and Eastern Europe; S1 Table).

ORCHIDEE-GM simulates the potential livestock density $D_{opt}$ on a 25 km × 25 km grid over Europe. This model output was re-aggregated to the NUTS level ($D_{reg}$ for NUTS-0 to NUTS-2) and weighted by the actual grassland area of each grid cell from the CORINE Land Cover map (CLC2000 [56]). The simulated number of grass-fed livestock ($N_{sim}$) is given by:

$$N_{sim} = f_{temp} \times D_{reg} \times A_{temp} + f_{perm} \times D_{reg} \times A_{perm} + f_{rough} \times D_{reg} \times A_{rough} \quad (14)$$

where $A_{temp}$, $A_{perm}$, and $A_{rough}$ are the area of temporary grassland, permanent grassland, and rough grazing respectively in each NUTS region (data from Eurostat; names in italics from
Eurostat terminology, see above); $f_{temp}$, $f_{perm}$, and $f_{rough}$ are the fraction of $D_{reg}$ that grassland of different types are assumed to support, and are set to 100%, 80%, and 10% for temporary grasslands, permanent grasslands, and rough grazing lands respectively according to their productivity. Newly sown temporary grassland has high productivity, and thus is assumed to support the simulated potential livestock density, whereas less productive rough grazing grasslands are assumed to receive little management and to support only 10% of $D_{reg}$ (these grasslands are usually located in mountainous areas with steep slopes and limited accessibility). Combining modeled livestock density with statistical grassland-area data, we produced a simulated number of animals for each region, $N_{sim}$ that can be compared with $N_{obs}$ available from the statistics.

Long-term regional time series of grassland area and livestock numbers ($N_{obs}$) are not available for the Eurostat NUTS regions. Therefore, country average $N_{obs}$ (also called NUTS-0) from 30 countries from the FAO statistical database [5] were used to evaluate modeled trends of $N_{sim}$ over the period from 1961 to 2009.

All the data used in this study are listed in Table 1, as well as the corresponding model results for comparison.

## Model-data agreement metrics

To assess model-data agreement for the spatial distribution of grassland productivity, we use the Pearson's product-moment correlation coefficient ($r$) with P-value. $r$ quantifies the proportion of the total variance in the observed data that can be explained by the model, given by:

$$r = \frac{\sum_{i=1}^{n} (P_i - \bar{P})(O_i - \bar{O})}{\sqrt{\sum_{i=1}^{n} (P_i - \bar{P})^2 \sum_{i=1}^{n} (O_i - \bar{O})^2}} \quad (15)$$

where $P_i$ is modeled data, $O_i$ is observed data, $\bar{P}$ is modeled mean, $\bar{O}$ is observed mean, and $n$ is the number of data-points. The P-value is the probability that the correlation coefficient is zero (null hypothesis). When the P-value is lower than the conventional 5% ($P < 0.05$) the correlation coefficient is statistically significant. The correlation coefficient ($r$) is also used with the P-value to assess the model-data agreement on the interannual variability of productivity time series (including observed and modeled productivity and summer NDVI). For all time series, model output, statistics and summer NDVI), long-term trends were removed (detrended).

To evaluate grass-fed livestock numbers in Europe, ordinary least squares linear regression with intercept forced to zero is performed between $N_{sim}$ and $N_{obs}$. A regression coefficient

<table><caption>Table 1. Model results and data used for comparison.</caption><thead><tr><th>Model results</th><th>Data</th><th>Source</th></tr></thead><tbody><tr><td rowspan="2">Grassland productivity</td><td>Statistical productivity</td><td>[41], assembled by Smit et al. [20]</td></tr><tr><td>GIMMS NDVI</td><td>[61]</td></tr><tr><td>Net primary productivity (NPP)</td><td>Site NPP observations</td><td>NPP database (Campioli, unpublished)</td></tr><tr><td>Livestock density</td><td>Global livestock distribution maps</td><td>FAO [62]</td></tr><tr><td rowspan="3">Regional grass-fed livestock numbers</td><td>Grassland area</td><td>[41]</td></tr><tr><td>Livestock numbers</td><td>[41]</td></tr><tr><td>Livestock diet</td><td>[3]</td></tr></tbody></table>

doi:10.1371/journal.pone.0127554.t001
(slope) close to unity with a high coefficient of determination (R²) and with P < 0.05 indicates a good model-data agreement.

All the statistics in this study were calculated using Interactive Data Language (IDL) software (Version 8.3 (linux x86_64 m64). (c) 2013, Exelis Visual Information Solutions, Inc. http://www.exelisvis.com/ProductsServices/IDL.aspx).

## Uncertainties of the potential grassland productivity estimation

The uncertainties of predictions from process models could be rather large resulted from model forcing data [63, 64], parameter values uncertainty [65], as well as model structure related uncertainty [66–68]. At a large geographical scale (Europe), a comprehensive assessment of uncertainty, such as using factorial designs [69] and Monte Carlo-type stratified sampling approach [70], requires many model runs and are therefore prohibitive for complex models (such as ORCHIDEE-GM) with a large number of parameters, a half-hourly time step and a high computational demand [71]. In this study, we identified four model parameters that likely substantially contribute to uncertainties of productivity simulations [63, 65, 69, 72–74]. These four sources of uncertainties, define 16 combinations given minimum and maximum values that define a range (± 20% approximately) around the standard values used in the control simulation. The uncertain settings that are tested by sensitivity simulations are (1) the constant *a* (Eq 3) determining the shape of saturating curve, (2) the response of photosynthetic capacity to nitrogen addition (parameter *Nadd*<sub>max</sub>, Eq 3), (3) the maximum rate of Rubisco carboxylase activity (*Vcmax*<sub>opt</sub>) and the maximum rate of photosynthetic electron transport (*Jmax*<sub>opt</sub>), and (4) the prescribed maximum specific leaf area (*SLA*<sub>max</sub> [29]). Simulations with the 16 factors combinations at the full geographical scale of this study (9237 grids) would require enormous computational time beyond our capability. Thus, evenly spaced 368 grid cells over our study area were selected for the uncertainty analysis, which were able to well represent the spatial distribution, magnitude and interannual variability of grasslands’ productivity (see S1 File for detail). Complete simulations (as described in section “Simulation set-up”) were conducted at these grids with these factors combinations (with minimum and / or maximum values for each factors; Table 2). The standard deviation (SD) of the simulated results is used to characterize and assess the uncertainties of potential productivity, grass-fed livestock density, and their trends.

## Results

### Distribution of potential grassland productivity

ORCHIDEE-GM (Fig 3A) simulates the highest potential grass production (annually harvested forage biomass per unit of cut grassland area) in the Atlantic maritime-climate zone (The

<table><caption>Table 2. Key model parameters for grassland productivity simulations and their ranges.</caption><thead><tr><th>Model Parameters</th><th>Unit</th><th>Standard value</th><th>Minimum</th><th>Maximum</th><th>Description</th></tr></thead><tbody><tr><td>a</td><td>unitless</td><td>0.75</td><td>0.6</td><td>0.9</td><td>proportions of intensively managed grasslands</td></tr><tr><td>N<sub>addmax</sub></td><td>percent</td><td>60%</td><td>40%</td><td>80%</td><td>the saturate status of N addition effect on photosynthetic capacity</td></tr><tr><td>Vcmax<sub>opt</sub> / Jmax<sub>opt</sub></td><td>µmol m<sup>-2</sup> s<sup>-1</sup></td><td>55 / 110</td><td>44 / 88</td><td>66 / 132</td><td>Vcmax<sub>opt</sub>: the maximum rate of Rubisco carboxylase activity Jmax<sub>opt</sub>: the maximum rate of photosynthetic electron transport</td></tr><tr><td>SLA<sub>max</sub></td><td>m<sup>2</sup> gC<sup>-1</sup></td><td>0.048</td><td>0.0384</td><td>0.0576</td><td>the prescribed maximum SLA</td></tr></tbody></table>

Note: Factors are modified by ±20% of standard value (except for *N*<sub>addmax</sub>, which was modified by ± 20% of absolute value). For each combination, minimum or maximum value of each factor is used, which forms 2⁴ = 16 factor combinations.

doi:10.1371/journal.pone.0127554.t002
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Biological potential productivities simulated by ORCHIDEE-GM of (A) potential net primary productivity (harvested dry matter in forage) in cut grasslands and of (B) potential net primary grazed dry matter used for animal intake in pasture. Both fields are average values for the period 1995–2004.</figcaption></figure>

doi:10.1371/journal.pone.0127554.g003

Netherlands, Belgium, Luxembourg, and west of Ireland) with values reaching up to 10 tonnes DM ha⁻¹ yr⁻¹, and in the lower elevation parts of the Alps (e.g., southern Germany, west of France, northern part of Austria and Switzerland, northeast Italy, and Slovenia). The model simulates a lower productivity in Mediterranean countries (e.g., Greece, southern Italy, Spain, and Portugal) and in Scandinavia (e.g., Norway, northern Sweden and northern part of Finland,) and in the high altitude grid cells of the Alps. The modeled distribution of biomass eaten by grazing animals follows the same pattern as the production of harvested forage grass (see Fig 3B). However, biomass grazing amounts to only about 78% of harvested forage biomass. The difference comes from the lower herbage use efficiency of grazing (28% of NPP) compared with cutting (40% of NPP), given the similar NPP values simulated by ORCHIDEE-GM for cut (508 g C m⁻² yr⁻¹) and grazed grasslands (577 g C m⁻² yr⁻¹).

When the simulated productivity at the pixel level is aggregated over the Eurostat administrative regions (NUTS-2 to country level, depending on data availability, see Fig 4D), a significant positive spatial correlation ($r = 0.6$; $P < 0.01$) is obtained between simulated and observed productivity across the 167 NUTS regions. However, ORCHIDEE-GM tends to simulate higher potential productivity than the actual productivity (yield) reported in the Smit et al. [20] in most regions (Fig 4). This result is logical because the model simulates the potential (maximum) productivity of permanent cut grassland, whereas Eurostat productivities are based on actual harvest data. Exceptions are northern Spain, Norway, and northern Sweden where ORCHIDEE-GM simulates lower productivity than that from statistics. The positive difference between simulated potential and actual productivities [20] (Fig 4C) is the largest for Mediterranean (e.g. southern France and southern Italy) and east European regions.

Observed and modeled NPP are comparable ($R^2 = 0.76, P < 0.001$, Table 3), showing that ORCHIDEE-GM is well able to reproduce the NPP of productive grasslands without nitrogen limitation (e.g. sites CZ-mok and IT-bea; see list in Table 3), whereas it tends to overestimate NPP (up to 1.6-fold) at sites with relatively low productivity. However, given the fact that some local conditions (such as soil properties, topographic features, etc.) are not considered in the model and that few observations are available for this exercise, any conclusions on the bias of the model in simulating NPP can only be drawn with low confidence.
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Spatial distribution of (A): potential grassland productivity simulated by ORCHIDEE-GM from cut grasslands (M), (B): actual grassland forage productivity from Smit et al. [20] (S), (C): relative discrepancy between them expressed as log (1+(M-S)/S), and (D): the NUTS administrative units at which statistical data are available. Simulated and "observed" data from statistics represent both a 10-year average from 1995 to 2004.</figcaption></figure>

doi:10.1371/journal.pone.0127554.g004

Distribution of livestock density

The modeled potential grass-fed livestock density is compared with the FAO gridded data product [62] in Fig 5. The highest observed densities (≥ 0.6 LSU ha⁻¹) are found in the regions with a maritime climate bordering the North Atlantic and low elevation areas in the Alps; Fig 5B). These two regions are also simulated to have the highest potential livestock density. It should be kept in mind that we model potential livestock densities based on grassland carrying capacity in Europe, whereas the FAO data provide livestock density over land suitable for live-livestock production (including mosaics of cropland and pasture).

Distribution of potential mowing frequencies and the length of grazing period

Fig 6 shows the spatial distribution of potential mowing frequencies (Fig 6A) and the length of grazing period (Fig 6B) simulated by ORCHIDEE-GM. They follow the same pattern as the
<table><caption>Table 3. Comparison of observed and modeled annual NPP at seven temperate grassland sites.</caption><thead><tr><th>Site ID</th><th>Latitude</th><th>Longitude</th><th>Altitude (m)</th><th>Year</th><th>Observed NPP (g C m<sup>-2</sup>yr<sup>-1</sup>)</th><th>Modeled NPP (g C m<sup>-2</sup>yr<sup>-1</sup>)</th><th>References</th><th>Vegetation Source</th></tr></thead><tbody><tr><td rowspan="3">CH-jur</td><td rowspan="3">47° 33' N</td><td rowspan="3">07° 34' E</td><td rowspan="3">520</td><td>1994</td><td>453</td><td>694</td><td>[95]</td><td>Semi-natural calcareous grassland</td></tr><tr><td>1995</td><td>554</td><td>679</td><td>[95]</td><td>Semi-natural calcareous grassland</td></tr><tr><td>1996</td><td>495</td><td>691</td><td>[95]</td><td>Semi-natural calcareous grassland</td></tr><tr><td rowspan="4">CZ-kam</td><td rowspan="4">49° 43' N</td><td rowspan="4">15°58' E</td><td rowspan="4">1424</td><td>1976</td><td>457</td><td>539</td><td>[96]</td><td>Moist meadow</td></tr><tr><td>1977</td><td>628</td><td>624</td><td>[96]</td><td>Moist meadow</td></tr><tr><td>1978</td><td>428</td><td>552</td><td>[96]</td><td>Moist meadow</td></tr><tr><td>1979</td><td>409</td><td>596</td><td>[96]</td><td>Moist meadow</td></tr><tr><td rowspan="2">CZ-mok</td><td rowspan="2">49° 01' N</td><td rowspan="2">14° 46' E</td><td rowspan="2">430</td><td>2007</td><td>777</td><td>670</td><td>[97]</td><td>Nutrient-poor wet grassland</td></tr><tr><td>2007</td><td>924</td><td>938</td><td>[97]</td><td>Nutrient-rich wet grassland<sup>a</sup></td></tr><tr><td>DE-gri*</td><td>50° 57' N</td><td>13° 30' E</td><td>375</td><td>2004</td><td>403</td><td>646</td><td>[98]</td><td>Temperate grassland</td></tr><tr><td rowspan="2">IT-bea*</td><td rowspan="2">46° 00' N</td><td rowspan="2">13° 01' E</td><td rowspan="2">43</td><td>2007</td><td>1159</td><td>1025</td><td>[99]</td><td>N<sub>2</sub>-fixing grassland used as fodder</td></tr><tr><td>2008</td><td>1109</td><td>976</td><td>[99]</td><td>N<sub>2</sub>-fixing grassland used as fodder</td></tr><tr><td>SE-kje*</td><td>60° 10' N</td><td>17° 38' E</td><td>30</td><td>1982–1983</td><td>765</td><td>586</td><td>[100]</td><td>Meadow fescue ley, lucerne ley<sup>b</sup></td></tr><tr><td>UK-pge</td><td>51° 49' N</td><td>0° 21' W</td><td>128</td><td>1970–1985</td><td>400</td><td>513</td><td>[101]</td><td>Mixed sward</td></tr></tbody></table>

doi:10.1371/journal.pone.0127554.t003

production of harvested forage grass and the biomass eaten by grazing animals (Fig 3). In the Atlantic maritime-climate zone (The Netherlands, Belgium, Luxembourg, and west of Ireland) and in the lower elevation parts of the Alps (e.g., southern Germany, west and central part of France, and Slovenia), cut grasslands could be mown as frequent as four times per year, while pasture could be grazed by more than 8 month of a year. Low mowing frequencies (less than 2 times per year) and short grazing season were simulated in Mediterranean countries (e.g., Greece, southern Italy, Spain, and Portugal) and in the mountain area and the high latitude of

<figure><img src="image_8.png" /><figcaption>Fig 5. Comparison between (A): potential livestock density from ORCHIDEE-GM expressed in livestock unit (LSU) per ha of grassland, under the modeled optimal cut-grazing management scheme, and (B): observed livestock density from FAO [62].</figcaption></figure>

doi:10.1371/journal.pone.0127554.g005

<footnote><sup>a</sup> For the nutrient-rich wet grassland at site CZ-mok and N<sub>2</sub>-fixing grassland at site IT-bea and SE-kje, site simulations were carried out with an addition of 300 kg N ha<sup>-1</sup> assuming there is no N limitation at the sites.</footnote>

<footnote><sup>b</sup> For the grassland at site SE-kje, 2 NPP observations with meadow fescue ley (740 g C m<sup>-2</sup>yr<sup>-1</sup>) and lucerne ley (790 g C m<sup>-2</sup>yr<sup>-1</sup>) respectively were averaged to compare with result from ORCHIDEE-GM.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig 6.** Spatial distribution of potential mowing frequencies (A) and the length of grazing period (B) simulated by ORCHIDEE-GM. Both fields are average values for the period 1995–2004.</figcaption></figure>

doi:10.1371/journal.pone.0127554.g006

Scandinavia (e.g., the mountain area of Norway, northern Sweden and northern part of Finland,) and in the high altitude grid cells of the Alps.

## Temporal evolution of grassland production and livestock density

ORCHIDEE-GM simulates an annual increase (0.36% ± 0.06% per year, P < 0.001) of potential grassland productivity over the last five decades (Fig 7). Because the model mostly neglects other factors, this trend is driven by CO₂ concentration, climate change, nitrogen fertilizer applications and nitrogen deposition. The relative contribution of these drivers of the trend in potential productivity are 97%, -3% and 15%, respectively (Table 4). The total contribution of the three factors is larger than 100%, because these variables are not independent. The productivity trend is surprisingly close to the mean annual trend in the production of sown grasses estimated from measurements (0.35% per year; reviewed by Smit et al. [20]), although it is unclear

<figure><img src="image_10.png" /><figcaption>**Fig 7.** Temporal variations of modeled grassland productivity and of livestock density from model simulations and grassland agricultural statistics (see text) respectively. LSU: livestock unit. All variables are spatial averages over European grasslands masked by CORINE Land Cover map (CLC2000 [56]).</figcaption></figure>

doi:10.1371/journal.pone.0127554.g007
<table><caption>Table 4. The contribution of rise in atmospheric CO₂, climate change and nitrogen fertilization (including nitrogen deposition) trends on trends in potential productivity.</caption><thead><tr><th rowspan="3">Productivity trend<br>T<sub>E4</sub><br>(t DM ha<sup>-1</sup> yr<sup>-2</sup>)</th><th colspan="2">Rise in CO<sub>2</sub></th><th colspan="4">Contribution</th><th colspan="2">Total contribution</th></tr><tr><th rowspan="2">T<sub>E4</sub>-T<sub>E5</sub></th><th rowspan="2">Percentage</th><th colspan="2">Climate change</th><th colspan="2">N fertilization and<br>deposition trends</th><th rowspan="2">Value</th><th rowspan="2">Percentage</th></tr><tr><th>T<sub>E4</sub>-T<sub>E6</sub></th><th>Percentage</th><th>T<sub>E4</sub>-T<sub>E7</sub></th><th>Percentage</th></tr></thead><tbody><tr><td>1.62×10<sup>-2</sup></td><td>1.56×10<sup>-2</sup></td><td>97%</td><td>-0.05×10<sup>-2</sup></td><td>-3%</td><td>0.24×10<sup>-2</sup></td><td>15%</td><td>1.76×10<sup>-2</sup></td><td>109%</td></tr></tbody></table>

Note: T<sub>E4</sub>, T<sub>E5</sub>, T<sub>E6</sub> and T<sub>E7</sub> are productivity trends from simulation E4, E5, E6 and E7 respectively.

doi:10.1371/journal.pone.0127554.t004

whether the observed production gain can be attributed to genetic or to environmental (e.g., atmospheric CO₂ and climate) factors. A smaller trend was found in the simulated potential livestock density, which increased at a rate of 0.32% ± 0.06% per year (P < 0.001). The livestock density from FAOstat shows a similar but larger trend of 0.33% per year (P < 0.001) between 1961 and 1990, and then a strong decline (2.6% per year, P < 0.001) from 1991. However this decrease in livestock density comes with an increase in milk and meat productivity per head, as a result of the increasing use of crop-feed products including imported feed such as soybean [75]. In 1984 the European Community introduced milk production quotas that contributed to a reduction in the dairy cow population in Europe. This was followed in 1991 by the Nitrates Directive (91/676/EEC) [76] that restricted the application of animal manure in “nitrate vulnerable zones” to a maximum of 170 kg N ha⁻¹ yr⁻¹, effectively capping livestock density in pastures at some 1.7 LSU ha⁻¹. These two policies have caused a decrease in livestock numbers in Europe. The uncertainty of the simulated trends above comes from 1-sigma standard deviation of the trends from 16 sensitivity tests.

Changes in grassland productivity between 1973 and 2005 are shown for France, Germany, Italy and the Czech Republic (Fig 8). In Germany, the data show a steep decline of productivity after 1985, which may reflect a change in the definition of pastures [20]. Thus for Germany, only the normalized series after 1985 were shown when statistics are consistent. For the above-listed countries, ORCHIDEE-GM reproduces the drops in productivity observed during the 1976 and 2003 droughts. Further, ORCHIDEE-GM is also able to reproduce the interannual variability of production, although Italy is an exception. High correlations between time series of modeled and observed [20, 41] productivity are found in France (r = 0.69, P < 0.01), Czech Republic (r = 0.52, P < 0.01), and the best results in Germany (r = 0.85, P < 0.01). The modeled annual productivity is also significantly correlated with summer grassland NDVI in France (r = 0.39, P < 0.05), but not in Germany (r = 0.31, P = 0.11; Table 5). We use the coefficient of variation (CV) to express the magnitude of the interannual variability in modeled and Eurostat productivities, and the summer NDVI (Table 5). In general, the CV of modeled productivity is larger than the one of Eurostat productivity. The CV of modeled productivity (range 17–30%) is much larger than that of summer NDVI (2–3%).

## Ratio of livestock density to grassland productivity in Europe

Table 6 compares ORCHIDEE-GM estimates and FAO/Eurostat statistics for the ratio of livestock density to grassland production for different Köppen-Geiger climate zones (S2 Fig; derived from Peel et al. [77]). The modeled ratios in different climate zones are very similar (ranging between 1.58 and 1.75 LSU per 10 tonnes of DM, and do not show significant changes over the simulation periods (1961–2010). However, in the statistical data, this ratio varies across climate zones. Regions with cold and temperate climate have low stocking
<figure><img src="image_11.png" /><figcaption>**Fig 8.** Normalized temporal evolution of modeled and observed productivity of grasslands in four European countries between 1973 and 2005. Observed productivity is derived from Smit et al. [20] based on Eurostat agricultural statistics data (see text).</figcaption></figure>

doi:10.1371/journal.pone.0127554.g008

density to grassland productivity ratios (0.75 ± 0.49 and 0.75 ± 0.48 LSU per 10 tonnes of DM, respectively), much less than that modeled by ORCHIDEE-GM. The highest ratio (1.49 ± 1.24 LSU per 10 tonnes of DM), close to the simulated one, is found in regions with a Mediterranean climate (Table 6).

## Evaluation of grass-fed livestock numbers in Europe

Regressing the simulated numbers of grass-fed animals ($N_{sim}$ in Eq 14) against the Eurostat data ($N_{obs}$ in Eq 13), gives a (spatial) coefficient of determination $R^2 = 0.76$ ($P < 0.01$) for the NUTS-2 regions (slope = 1.00; NUTS-1 for Germany; Fig 9). Country scale (NUTS-0) comparison (Fig 10) shows the mean value and standard deviation of the observed and modeled grass-fed livestock numbers during the period 1990–2010. When the comparison is made for regions grouped into Köppen-Geiger climate zones (S2 Fig; [77]), $N_{sim}$ is comparable to $N_{obs}$ in regions with temperate climate (slope close to 1, $R^2 > 0.8$, $P < 0.01$; Fig 9). However, ORCHIDEE-GM tends to underestimate the number of grass-fed livestock in regions with cold high latitude climate (e.g., Norway) and in Mediterranean countries (e.g., Greece, Italy, Spain, Portugal).
<table><caption>Table 5. Mean value, standard deviation (SD), and coefficient of variation (CV) of modeled and statistical estimates of grassland productivity and NDVI time series for four European countries, as well as correlation coefficient (r) between each pair.</caption><thead><tr><th></th><th>Czech Republic</th><th>France</th><th>Germany</th><th>Italy</th></tr></thead><tbody><tr><td colspan="5">Modeled Production (M)</td></tr><tr><td>Mean (t DM ha<sup>-1</sup> yr<sup>-1</sup>)</td><td>5.82</td><td>6.41</td><td>7.44</td><td>2.29</td></tr><tr><td>SD (t DM ha<sup>-1</sup> yr<sup>-1</sup>)</td><td>1.50</td><td>1.48</td><td>1.26</td><td>0.70</td></tr><tr><td>CV (%)</td><td>26</td><td>23</td><td>17</td><td>30</td></tr><tr><td colspan="5">Statistical Production (S)</td></tr><tr><td>Mean (t DM ha<sup>-1</sup> yr<sup>-1</sup>)</td><td>3.26</td><td>4.14</td><td>6.3</td><td>1.54</td></tr><tr><td>SD (t DM ha<sup>-1</sup> yr<sup>-1</sup>)</td><td>0.45</td><td>0.51</td><td>0.47</td><td>0.19</td></tr><tr><td>CV (%)</td><td>14</td><td>12</td><td>7</td><td>12</td></tr><tr><td colspan="5">Summer NDVI (June–August) (N)</td></tr><tr><td>Number of point</td><td>3</td><td>688</td><td>309</td><td>18</td></tr><tr><td>Mean</td><td>0.84</td><td>0.82</td><td>0.83</td><td>0.74</td></tr><tr><td>SD</td><td>0.02</td><td>0.02</td><td>0.02</td><td>0.02</td></tr><tr><td>CV (%)</td><td>3</td><td>2</td><td>2</td><td>3</td></tr><tr><td colspan="5">r</td></tr><tr><td>M vs. S (1973–2005<sup>a</sup>)</td><td><strong>0.52</strong></td><td><strong>0.69</strong></td><td><strong>0.85</strong></td><td>0.12</td></tr><tr><td>M vs. N (1982–2010<sup>a</sup>)</td><td>0.22</td><td><strong>0.39</strong></td><td>0.31</td><td>0.20</td></tr></tbody></table>

Note: To represent interannual variability, long-term trends were removed in all series (detrended) before the calculation of *r*. M vs. S, *r* between modeled and statistical grassland productivity; M vs. N, *r* between modeled productivity and NDVI. *r* values in bold indicate a statistically significant correlation (with *P* < 0.05). Modeled productivity time series are from 1961 to 2010; statistical productivity time series are from 1973 to 2005 (for Germany from 1985 to 2005); NDVI time series are from 1982 to 2010.

doi:10.1371/journal.pone.0127554.t005

## Discussion

### Spatial distribution of grassland productivity

In this study, the grassland productivity simulated by ORCHIDEE-GM represents a biological potential estimate (in equilibrium with climate and CO₂ concentration), but neglects soil nutrient limitations on herbage production (not alleviated by nitrogen fertilization and deposition, a factor being accounted for in ORCHIDEE-GM simulations) and also ruminant nutrition limitations caused by low herbage quality. On the contrary, the statistical data compile actual production (often estimated from expert judgments) in each region under different types of management, which is often restricted by limiting soil nutrients and by low herbage quality. This difference is a major cause of productivity overestimation by the model. Another cause could be the possible overestimation of NPP by ORCHIDEE-GM (Table 3, see also [78]). NPP

<table><caption>Table 6. Ratio of livestock density to grassland productivity in Europe.</caption><thead><tr><th rowspan="2">Köppen-Geiger climate zones</th><th colspan="2">Livestock density / productivity LSU per 10 tonnes of DM</th></tr><tr><th>Model</th><th>Statistics</th></tr></thead><tbody><tr><td>Clod climate</td><td>1.58±0.11</td><td>0.75±0.49</td></tr><tr><td>Temperate climate</td><td>1.72±0.15</td><td>0.75±0.48</td></tr><tr><td>Mediterranean climate</td><td>1.75±0.13</td><td>1.49±1.24</td></tr></tbody></table>

Note: See S2 Fig for map of Köppen-Geiger climate zones.

doi:10.1371/journal.pone.0127554.t006

<footnote><sup>a</sup> The periods of time series that use for the calculation of *r* and *P*-value.</footnote>
<figure><img src="image_12.png" /><figcaption>Fig 9. Comparison of modeled (Nmodel) and ‘observed’ (Nobs) grass-fed livestock numbers (‘observed’ numbers are inferred by combining Eurostat agricultural statistics and additional data in Eq 13) at the scale of each NUTS-2 region for the period 1990–2010 using ordinary least squares linear regression forced through the origin. One point is the average over a NUTS-2 region. A regression coefficient (slope) close to unity with a high coefficient of determination (R²) with P < 0.05 indicates a good model-data agreement. Climate zones follow the Köppen-Geiger classification (see text).</figcaption></figure>

doi:10.1371/journal.pone.0127554.g009

established from ecological measurements, separating aboveground from belowground productivity, is a more accurate indication of annual grass growth than grassland production estimates based on statistics. However, given the limited number of NPP measurements used to validate the model, its accuracy for grassland NPP projections cannot be fully ascertained.

Nevertheless, the spatial distribution of modeled potential productivity agrees well with statistics (Fig 4). ORCHIDEE-GM realistically reproduces (i) high productivity in humid regions,
<figure><img src="image_13.png" /><figcaption>Fig 10. Comparison of modeled (Nmodel) and ‘observed’ (Nobs) (inferred from Eurostat and other data in Eq 13) grass-fed livestock numbers at the scale of NUTS-0 regions in Europe (countries) for the period 1990–2010. LSU: livestock unit.</figcaption></figure>

doi:10.1371/journal.pone.0127554.g010

oceanic regions, and in wet mountain regions at low to mid-elevation, and (ii) low productivity in dry regions with a Mediterranean climate heavily affected by drought. This pattern confirms the major role played by water limitation in the productivity of temperate grassland ecosystems [79–84]. Then again, site-specific properties such as soil texture, terrain slope, nutrient availability, and grass species are not taken into account in the current model version, which reduce its ability to capture the actual details of the spatial distribution of grassland productivity. For example, we have a large overestimate of productivity in southern Europe (Fig 4) possibly due to the model’s lack of consideration of soil nutrients (S3A Fig) and to the low aboveground productivity of some grass species. On the other hand, overestimation of grassland productivity in southeast European countries (e.g., Romania and Bulgaria) may be the result of less grazing management reducing herbage productivity [7] and steppic conditions in some areas. Detailed spatial information about soil fertility and landscape management intensity would be helpful for improving the model capacity to simulate grassland productivity. On the other hand, some data limitations in terms of quality, transparency and detail [20] could be another cause of model-data discrepancies.

## Decadal trends and interannual variability of grassland productivity

For European grasslands, ORCHIDEE-GM estimates an average annual increase of potential productivity over the last five decades (0.36% ± 0.06% per year, Fig 7) which in the model is explained by the effects of elevated CO₂, climate change, nitrogen fertilization, and nitrogen deposition. Changes in other management drivers (except nitrogen fertilization) in this period are not included in the model because the data (e.g., plant breeding) are not available. Meta-analyses revealed that nitrogen fertilization itself could stimulate grass biomass by around 50% [36, 37], and the positive biomass responses could be significantly enhanced when CO₂
enrichment and/or the addition of other nutrients (such as phosphorus) were taken into ac-
count [36, 85]. Such a strong effect of fertilization, however, was not observed in regional-scale
productivity statistics, nor in modeled potential productivity. We believe there are two reasons
for this lack of response: i) Intensive fertilization that would strongly improve productivity has
not been fully applied to all grasslands in Europe (a large part of the total grassland area was
fertilized with only 0–40 kg N ha⁻¹ yr⁻¹, S1 Fig). ii) Fertilization and the improvement of grass
species by plant breeding were introduced only a few decades ago in sown (i.e. temporary)
grasslands, and with apparently limited effects on productivity. On the other hand, it has been
shown that water stress is a major limitation on grassland productivity [79, 82, 83] and a small
number of grasslands are irrigated in Europe [86]. Thus, in contrast to arable agriculture (an-
nual crops), grasslands, which are still largely semi-natural in Europe, display changes in pro-
ductivity over recent decades that are primarily controlled by climatic and atmospheric factors
and not by management. Nevertheless, despite being increasingly decoupled from local grass
productivity (e.g., increased use of crops and compound feeds in ruminant diet [2]), the live-
stock production (i.e., the potential grass-fed livestock density determined by potential produc-
tivity in this study) is still significantly influenced by trends in atmospheric CO₂ concentration
and climate (Table 4). The effect of these trends have to be taken into account for future scenar-
ios, because a change of grassland productivity will imply significantly less of more feed com-
plements that strongly related to agricultural policy.

Land-use change could be another significant and non-modeled factor that impacts long-
term trends of grassland productivity in some EU countries. Transition to arable crops and
abandonment can both mark a land-use change affecting grasslands. For example, in the Czech
Republic, besides the impact of unfavorable climate factors captured by our model (Fig 8), grass-
lands were also abandoned during the transition to a market economy early 1990s, causing a
drop in productivity between 1989 and 1992. In Italy, the evolution of grassland productivity has
also been impacted by major land-use changes. The number of grazing animals decreased owing
to the abandonment of low-productive dry rangelands [87]. At the same time, grassland areas
with favorable conditions have been partly converted to fodder maize and to cash crops [87].

In certain regions, changes in the management of grasslands, such as fertilization, may have
reduced the variability of productivity between years. However, significant interannual vari-
ability of grassland production was observed (Fig 8) in some countries because of climatic vari-
ability. In this study, we have shown that ORCHIDEE-GM is able of capturing the interannual
variability of grassland productivity (Table 5).

Grassland productivity uncertainties from model parameters

The key model inputs and parameters considered for uncertainty assessment in this study cause an uncertainty of potential grassland productivity (on average 1-sigma error of 0.8 tonnes DM ha⁻¹ yr⁻¹), and of its decadal trends (on average 1-sigma error of 0.06% per year), which is 19% of the average potential grassland productivity (4.4 tonnes DM ha⁻¹ yr⁻¹), and 17% of its average trend (0.36% per year). Similar magnitude of uncertainties, expressed in percentage of average value, was found on the grass-fed livestock density (19%) and the decadal trends of it (13%). Within the total uncertainty on the trends of potential productivity, the uncertainties caused by the parameters related to the response of grass photosynthesis to nitrogen addition effect (constant *a* and *Nadd*<sub>max</sub>) are small (contributing an uncertainty of 0.03% per year for each of them). It implies that the maximum response of photosynthetic capacity to nitrogen addition and the shape of saturating curve (Eq 3) have very limited effects on the trends of potential productivity estimate. The larger uncertainties come from parameters representing photosynthetic and morphological plant traits (*Vcmax*<sub>opt</sub> / *Jmax*<sub>opt</sub>, and *SLA*<sub>max</sub>), which contribute
a productivity uncertainty of 0.05% per year for each of them. The uncertain values of these pa-
rameters could be one of the sources for model-data disagreement when simulating C fluxes at
sites [29]. However, these PFT-specific average plant functional traits in ORCHIDEE, in reality,
are highly site-specific, but not with large variation on average. To narrow the uncertainty of
the traits-related parameters, improved observation datasets on both mean value at community
level (contrary to species level) and on spatial distribution are required. Meanwhile, these traits
are tightly correlated with leaf nitrogen concentrations [88] suggesting a possible way to reduce
the uncertainty by fully coupling nitrogen and carbon cycles in terrestrial ecosystems [39].

Livestock density and its ratio with grassland productivity

Spatial distribution of grass-fed livestock density estimated by ORCHIDEE-GM is in agree-
ment with the general pattern of actual ruminant livestock distribution given by FAO [62] (Fig
5). However, because FAO data are densities expressed as livestock per unit of land area suit-
able for livestock production, data from regions with reduced grassland cover, or with low in-
tensity of grassland utilization by livestock, should not be compared with ORCHIDEE-GM
output. The actual livestock density is also restricted by soil quality, terrain slope, and by socio-
economic factors. Regions with low livestock density (< 0.2 LSU ha⁻¹) are distributed across
areas with high levels of soil workability constraints (i.e., soil conditions to ploughing, level 2–4
in S3A Fig; moderate constraints or (very) severe constraints) and high terrain slope (> 16% in
S3B Fig), including Norway, Sweden, Scotland, Mediterranean regions, and mountainous
areas. Those soil and terrain conditions prevent the use of mechanized mowing and, hence, the
harvesting of large stocks of conserved forage, which is a prerequisite if animals are to be
stocked at a high density during the grassland growing season. Moreover, the low livestock
density in eastern Europe may be restricted by socioeconomic conditions as well.

The model simulations of the ratio of livestock density to grass production are derived
under the assumption of maximum densities supported by the productivity of each grid cell
and the maximum utilization of the grassland productivity by grazing and mowing. This as-
sumption implies a fully utilized biological potential. When derived from statistics this ratio is
an indicator of the grassland utilization intensity. The ratio from statistics in Mediterranean re-
gions is close to that derived from model simulations, indicating an intensive use, or even an
over-utilization (i.e., overgrazing), of grasslands. The low value of the statistical ratio in tem-
perate Europe (especially in eastern Europe; data not shown) is related to an extensive use of
grasslands, resulting from unfavorable conditions in mountainous areas (Slovakia, Romania
and Bulgaria; S3B Fig), to the widespread abandonment of grassland areas in former commu-
nist countries [89], and to non-intensive management.

Regional numbers of grass-fed livestock

Grass-fed livestock numbers simulated by ORCHIDEE-GM are comparable with the data in-
ferred from Eurostat and FAO statistics, yet with some discrepancies (Figs 9 and 10). Model-
data misfits in the numbers of grass-fed livestock units could be caused by: i) the assumption
used in Eqs 4 and 5 that grasslands are either cut or grazed, whereas a number of grassland
areas in Europe have a mixed use (e.g., combining a spring or summer cut and several grazing
periods per year [7, 75], ii) the fact that cut grasslands sustain more head of livestock per hect-
are than grazed grasslands (since the fraction of the total aboveground dry matter harvested by
mowing is greater than with grazing; see Fig 3 and [32]). The proportion of cut grasslands
(1-F_opt, in Eq 4) in the real world could be higher than our optimal value in some regions (e.g.,
Denmark, northern Germany and Austria [90]), given the prevalence in these countries of
zero-grazing for dairy cattle. If so, grasslands in these regions are able to feed more livestock
than the numbers estimated by the model. iii) Animal diet composition in each region is different because of differences in the ratio of cropland to grassland area, in feed imports and exports, and in socioeconomic conditions. Thus, the coarse estimate of the fractional component of grass in animal diet used in our model-data comparison (only aggregated data for western Europe and eastern Europe were available in [3], the only data source about the grass-fed fraction of animals to the best of our knowledge) could be a source of inaccuracy in estimating grass-fed livestock numbers. iv) The factors limiting the utilization intensity for each kind of grassland ($f_{temp}$, $f_{perm}$ and $f_{rough}$ in Eq 14) are assumed fixed over all of Europe. However, in reality these parameters may be different for different regions [87]. For example, in the continental, boreal or Mediterranean zones, permanent grasslands tend not to be plowed and are less intensively managed, thus with an $f_{perm}$ lower than the default value (80%), which may induce a model overestimation of livestock density in these regions. v) In Mediterranean zones, where livestock is mostly sheep and goats, animals sometimes graze in shrubland or forest [91–93], which is not considered in our model. For the Mediterranean zones, the tendency of ORCHIDEE-GM to underestimate grass-fed livestock numbers (Fig 9) while overestimating productivity (Fig 4), could be partly explained by the particular farm structure (with sheep and goats as the major ruminant type [62]) and by overgrazing in some regions (e.g. [94]). Meanwhile, this result suggests that integrated productivity data (i.e., data in [20]) in Mediterranean zones might underestimate grassland productivity.

Accuracy of data is another important factor inducing discrepancies with model results: i) For some countries, grassland areas from FAOstat [5] and Eurostat [41] are significantly different (S2 Table). For example, grassland area from FAOstat [5] is larger than that from Eurostat [41] for permanent grassland in Bulgaria (by 39%), Greece (by 499%), Hungary (by 39%), Italy (by 29%), Spain (by 25%), Switzerland (by 77%) and United Kingdom (by 16%), and for temporary grassland in Germany (by 440%) and Spain (by 94%). This may partly explain the underestimation of grass-fed livestock numbers in these countries, because in this study grassland area was taken from Eurostat [41] (Eq 14). ii) For some countries, the area of grassland is very different from year to year. In Bulgaria (S3 Table) for instance, the areas of grasslands change between each survey, which may be due to changes in the criteria used for the surveys.

## Conclusion and Outlook

This study provides an estimation of potential productivity and potential livestock density over European grasslands by using the process-based ecosystem global vegetation model ORCHIDEE-GM with a representation of grassland management. When compared to statistical data from Eurostat or FAOstat, ORCHIDEE-GM can distinguish the spatial distribution of regions with low or high productivity and livestock density. However, the model tends to overestimate absolute values of productivity and livestock density because it was built and run to reproduce a biological potential. The model-data discrepancies reflect site specific properties such as soil quality, terrain slope, nutrient availability, grass species, management intensity, land-use change and socioeconomic factors that constrain productivity to be below the biological potential. Including these limitations in the model would probably allow future studies to estimate productivity more accurately.

The effects of elevated CO₂, climate change, nitrogen fertilization and nitrogen deposition are well reflected in the temporal evolution of productivity simulated by the model. But at regional scale the evolution of managed grassland productivity is also determined by land-use change, and by land management intensity changes, both reflecting agricultural policy and socio-economical drivers. ORCHIDEE-GM is able to reproduce the interannual variability of productivity in three of the four countries where long-term data are available (the exception is
Italy), making it likely that the model will be able to give useful predictions of the dynamics of potential productivity of European grasslands under future climate change.

In addition, grass-fed livestock numbers can be realistically reproduced by ORCHIDEE-GM using simple equations and parameterization. Regional model-data misfits exist due to uncertainties caused by different farming strategies (e.g., mowing or grazing, type of animal species), diet composition, management intensity, as well as the uncertainties in the measured data. ORCHIDEE-GM, with its capability to account for grassland management and livestock fed by grass, could contribute to our understanding of future changes in the European livestock sector and of climate impacts on grassland, as well as atmospheric feedbacks through GHG emissions and soil carbon balance.

## Supporting Information

**S1 File. Grids selection for productivity uncertainty analysis.**

(DOCX)

**S1 Fig. Spatial distribution of the amount of nitrogen addition (including both mineral and organic nitrogen fertilizer) over European grasslands for the year 2010.** The grey areas in EU27 states indicate that nitrogen addition was not applied. The blank grids for Croatia, Norway, and Switzerland indicate that the data were not available.

(TIFF)

**S2 Fig. Köppen-Geiger climate classification map of Europe adapted from Peel et al. [77].**

Cold climate, also as tundra climate (ET); Temperate climate, including moist mid-latitude climate with cold winters (Df), and with mild winters (Cf); Mediterranean climate (Cs). ET, Df, Cf, and Cs are original climate types used in Peel et al. [77].

(TIFF)

**S3 Fig. (A): The soil workability constraints, and (B): the median terrain slope class from the Global Agro-Ecological Zones assessment (GAEZ, http://www.fao.org/nr/gaez/en/).**

The soil workability constraints are classified into four levels: 1, no or slight constraints; 2, moderate constraints; 3, severe constraints; and 4 very severe constraints.

(TIFF)

**S1 Table. Fraction of grass fed in ruminant livestock diet (%, dry matter basis) [3].**

(DOCX)

**S2 Table. Grassland area of European countries from FAO [5] and Eurostat [41] (Unit: 1000 hectare).**

(DOCX)

**S3 Table. Area of grassland of Bulgaria from Eurostat [39] (unit: hectare).** The year indicates the year of farm surveys.

(DOCX)

## Acknowledgments

We would like to thank Dr. Bond-Lamberty and the anonymous reviewer for their valuable comments on the manuscript. We greatly thank Dr. Harm Smit and Dr. Stefan Wirsenius, who provided data on grassland productivity and on livestock diet, respectively. We thank Dr. Benjamin Poulter, Dr. Shushi Peng and Mr. Ulrich Eza for their help on processing downscaled climate data, NDVI data and fertilization data, respectively. We would like to thank Dr. John Gash for his effort on English language editing.
## Author Contributions

Conceived and designed the experiments: JC N. Viovy N. Vuichard PC JFS. Performed the experiments: JC N. Viovy N. Vuichard. Analyzed the data: JC. Contributed reagents/materials/ analysis tools: N. Vuichard MC KK RM AL. Wrote the paper: JC PC MC KK AL.

## References

1. Alexandratos N, Bruinjsma J (2012) World Agriculture Towards 2030/2050. The 2012 Revision FAO, Rome.

2. Hasha G (2002) Livestock Feeding and Feed Imports in the European Union—A Decade of Change. Electronic Outlook Report from the Economic Research Service. United States Department of Agriculture.

3. Wirsenius S (2000) Human use of land and organic materials: modeling the turnover of biomass in the global food system. Department of Physical Resource Theory, Chalmers University of Technology, Göteborg University: Göteborg, Sweden. 2000.

4. Leip A, Achermann B, Billen G, Bleeker A, Bouwman AF, de Vires W, et al. (2011) Integrating nitrogen fluxes at the european scale, in: Sutton M, Howard C, Erisman JW, Billen G, Bleeker A, van Grinsven H, Grennfelt P, Grizzetti B, editors. European Nitrogen Assessment. Cambridge University Press, Cambridge, UK, pp. 345–376. doi: 10.129/2003GB002060

5. FAOstat. http://faostat3.fao.org

6. Hadjigeorgiou I, Osoro K, Fragoso de Almeida JP, Molle G (2005) Southern European grazing lands: production, environmental and landscape management aspects. Livest Prod Sci 96: 51–59.

7. Marriott CA, Fothergill M, Jeangros B, Scotton M, Louault F (2004) Long-term impacts of extensification of grassland management on biodiversity and productivity in upland areas: A review. Agronomie 24: 447–461

8. Soussana J-F, Luescher A (2007) Temperate grasslands and global atmospheric change: a review. Grass Forage Sci 62: 127–134

9. Rötter R, van de Geijn SC (1999) Climate change effects on plant growth, crop yield and livestock. Climatic Change 43: 651–681.

10. Wand SJE, Midgley GF, Jones MH, PS Curtis (1999) Responses of wild C4 and C3 grass (Poaceae) species to elevated atmospheric CO₂ concentration: a meta-analytic test of current theories and perceptions. Global Change Biol 5: 723–741. doi: 10.1046/j.1365-2486.1999.00265.x

11. Ainsworth EA, Long SP (2005) What have we learned from 15 years of free-air CO₂ enrichment (FACE)? A meta-analytic review of the responses of photosynthesis, canopy properties and plant production to rising CO₂. New Phytol 165: 351–372. doi: 10.1111/j.1469-8137.2004.01224.x PMID: 15720649

12. Casella E, Soussana J-F (1997) Long-term effect of CO₂ enrichment and temperature increase on the carbon balance of a temperate grass sward. J Exp Bot 48: 1309–1321.

13. Aeschlimann U, Nösberger J, Edward PJ, Schneider MK, Richter M, Blum H (2005) Responses of net ecosystem CO₂ exchange in managed grassland to longterm CO₂ enrichment, N fertilization and plant species. Plant Cell Environ 28: 823–833.

14. Izaurralde RC, Thomson AM, Morgan JA, Fay PA, Polley HW, Hatfield JL (2011) Climate impacts on agriculture: implications for forage and rangeland production. Agron J 103: 371–381.

15. Morgan JA, Pataki DE, Korner C, Clark H, Del Grosso SJ, Grunzweig JM, et al. (2004) Water relations in grassland and desert ecosystems exposed to elevated atmospheric CO₂. Oecologia 140: 11–25. PMID: 15156395

16. Jones MB, Donnelly A (2004) Carbon sequestration in temperate grassland ecosystems and the influence of management, climate and elevated CO₂. New Phytol 164: 423–439. doi: 10.1111/j.1469-8137.2004.01201.x

17. Lee J (1983) The spatial distribution of grassland production in Europe. In: Proceedings of the 9th General Meeting of the European Grassland Federation, 11–20.

18. Hume CJ, Corrall AJ (1986) A quantitative study of the effects of climate variability and climatic change on herbage production from intensively managed grassland in Western European countries. Final Report, EEC Contract No. CLI-058-8 I-UK(H) (1986).

19. Topp CFE, Doyle CJ (2004) Modelling the comparative productivity and profitability of grass and legue systems of silage production in northern Europe. Grass Forage Sci 59: 274–292.
20. Smit HJ, Metzger MJ, Ewert F (2008) Spatial distribution of grassland productivity and land use in Europe. Agr Syst 98: 208–219.

21. Aranjuelo I, Irigoyen JJ, Pérez P, Martínez-Carrasco R, Sánchez-Díaz M (2005) The use of temperature gradient tunnels for studying the combined effect of CO₂, temperature and water availability in N₂ fixing alfalfa plants. Ann Appl Biol 146: 51–60.

22. Henry HAL, Cleland EE, Field CB, Vitousek PM (2005) Interactive effects of elevated CO₂, N deposition and climate change on plant litter quality in a California annual grassland. Oecologia 142: 465–473. PMID: 15558326

23. Tubiello FN, Soussana J-F, Howden SM (2007) Crop and pasture response to climate change. P Natl Acad Sci USA 104: 19686–19690. PMID: 18077401

24. Soussana J-F, Graux Al, Tubiello FN (2010) Improving the use of modelling for projections of climate change impacts on crops and pastures. J Exp Bot 61: 2217–2228. doi: 10.1093/jxb/erp100 PMID: 20410317

25. Fischer G, van Velthuizen H, Shah M, Nachtergaele FO (2002) Global agroecological assessment for agriculture in the 21st century: methodology and results. IIASA, Laxenburg, Austria.

26. Bouwman AF, Kram T, Klein_Goldewijk K (2006) Integrated Modelling of Global Environmental Change. An overview of IMAGE 2.4. Netherlands Environmental Assessment Agency, Bilthoven, The Netherlands.

27. Van Ittersum MK, Ewert F, Alkan Olsson J, Andersen E, Bezlepkina I, Brouwer F, et al. (2008) Integrated assessment of agricultural and environmental policies—a modular framework for the EU (SEAMLESS). Agri Syst 96: 150–165.

28. Vuichard N, Ciais P, Viovy N, Calanca P, Soussana J-F, (2007) Estimating the greenhouse gas fluxes of European grasslands with a process-based model: 2. Simulations at the continental level. Global Biogeochem Cy 21: GB1005, doi: 10.1029/2005GB002612

29. Chang JF, Viovy N, Vuichard N, Ciais P, Wang T, Cozic A, et al. (2013) Incorporating grassland management in ORCHIDEE: model description and evaluation at 11 eddy-covariance sites in Europe. Geosci Model Dev 6: 2165–2181 doi: 10.5194/gmd-6-2165-2013

30. Krinner G, Viovy N, de Noblet-Ducoudre N, Ogee J, Polcher J, Friedlingstein P, et al. (2005) A dynamic global vegetation model for studies of the coupled tmosphere-biosphere system. Global Biogeochem Cy 19: GB1015, doi: 10.1029/2003GB002199

31. Riedo M, Grub A, Rosset M, Fuhrer J (1998) A pasture simulation model for dry matter production, and fluxes of carbon, nitrogen, water and energy. Ecol Model 105: 141–183.

32. Vuichard N, Soussana J-F, Ciais P, Viovy N, Ammann C, Calanca P, et al. (2007) Estimating the greenhouse gas fluxes of European grasslands with a process-based model: 1. Model evaluation from in situ measurements. Global Biogeochem Cy 21: GB1004, doi: 10.1029/2005GB002611

33. Graux Al, Gaurut M, Agabriel J, Baumont R, Delagarde R, Delaby L, et al. (2011) Development of the Pasture Simulation Model for assessing livestock production under climate change. Agr Ecosyst Environment 144: 69–91.

34. Evans JR (1989) Photosynthesis and nitrogen relationships in leaves of C plants. Oecologia (Berlin) 78: 9–19.

35. Garnier E, Salager J-L, Laurent G, Sonié L (1999) Relationships between photosynthesis, nitrogen and leaf structure in 14 grass species and their dependence on the basis of expression. New Phytol 143: 119–129. doi: 10.1046/j.1469-8137.1999.00426.x

36. Xia J, Wan S (2008) Global response patterns of terrestrial plant species to nitrogen addition. New Phytol 179: 428–439. doi: 10.1111/j.1469-8137.2008.02488.x PMID: 19086179

37. Le Bauer DS, Treseder KK (2008) Nitrogen limitation of net primary productivity in terrestrial ecosystems is globally distributed. Ecology 89: 371–379. PMID: 18409427

38. Bardgett RD, Wardle DA (2003) Herbivore-mediated linkages between aboveground and below-ground communities. Ecology 84:2258–2268.

39. Zaehle S, Friend A (2010) Carbon and nitrogen cycle dynamics in the O-CN land surface model: 1. Model description, site-scale evaluation, and sensitivity to parameter estimates. Global Biogeochem Cy 24, GB1005, doi: 10.1029/2009GB003521

40. Bai Y, Wu J, Clark CM, Naeem S, Pan Q, Huang J, et al. (2010) Tradeoffs and thresholds in the effects of nitrogen addition on biodiversity and ecosystem functioning: evidence from inner Mongolia Grasslands. Global Change Biol 16: 358–372. doi: 10.1111/j.1365-2486.2009.01950.x

41. Eurostat-AFF. Eurostat Agriculture, forestry and fisheries. http://epp.eurostat.ec.europa.eu

42. Eurostat Glossary: Livestock unit (LSU), http://ec.europa.eu/eurostat/statistics-explained/index.php/Glossary:Livestock_unit_(LSU)
43. IPCC (Intergovernmental Panel on Climate Change) (2006) 2006 IPCC Guidelines for National Greenhouse Gas Inventories, Prepared by the National Greenhouse Gas Inventories Programme, Eggleston HS, Buendia L, Miwa K, Ngara T, Tanabe K, editors. Published: IGES, Japan.

44. IUCN/UNEP/WWF (1991) Caring for the Earth: A Strategy for Sustainable Living, Gland, Switzerland.

45. Beer C, Weber U, Tomelleri E, Carvalhais N, Mahecha MD, Reichstein M (2014) Harmonized European Long-Term Climate Data for Assessing the Effect of Changing Temporal Variability on Land—Atmosphere CO2 Fluxes. J Climate 27: 4815–4834. doi: 10.1175/JCLI-D-13-00543.1

46. Dee DP, Uppala SM, Simmons AJ, Berrisford P, Poli P, Kobayashi S, et al. (2011) The ERA-Interim reanalysis: configuration and performance of the data assimilation system. QJR Meteorol Soc 137: 553–597. doi: 10.1002/qj.828

47. Weedon GP, Gomes S, Viterbo P, Österle H, Adam JC, Bellouin N, et al. (2010) The WATCH Forcing Data 1958–2001: a meteorological forcing dataset for land surface- and hydrological models. WATCH Tech. Rep. 22, 41p (available at www.eu-watch.org/publications).

48. Weedon GP, Gomes S, Viterbo P, Shuttleworth WJ, Blyth E, Österle H, et al. (2011) Creation of the WATCH Forcing data and its use to assess global and regional reference crop evaporation over land during the twentieth century. J Hydrometeorol 12: 823–848, doi: 10.1175/2011JHM1369.1

49. New M, Lister D, Hulme M, Makin I (2002) A high-resolution data set of surface climate over global land areas. Climate Research 21:1–25.

50. Leip A, Britz W, Weiss F, de Vries W (2011) Farm, land, and soil nitrogen budgets for agriculture in Europe calculated with CAPRI. Environ Pollut 159: 3243–3253, doi: 10.1016/j.envpol.2011.01.040 PMID: 21420769

51. Leip A, Weiss F, Lesschen JP, Westhoek H (2014) The nitrogen footprint of food products in the European union. J Agric Sci 152: S22–33

52. Leip A, Marchi G, Koeble R, Kempen M, Britz W, Li C (2008) Linking an economic model for European agriculture with a mechanistic model to estimate nitrogen and carbon losses from arable soils in Europe, Biogeosciences 5: 73–94, doi: 10.5194/bg-5-73-2008

53. AGRESTE statistics. http://agreste.agriculture.gouv.fr

54. Tenkorang F, Lowenberg-DeBoer J (2008) Forecasting long-term global fertilizer demand. Nutr Cycl Agroecosys 83: 233–247.

55. Eurostat (2007) Regions in the European Union—Nomenclature of territorial units for statistics—NUTS 2006/EU-27.

56. Büttner G, Feranec J, Jaffrain G, Mari L, Maucha G, Soukup T (2004) The CORINE Land Cover 2000 Project. EARSeL eProceedings 3: 331–346.

57. Scurlock JMO, Cramer W, Olson RJ, Parton WJ, Prince SD (1999) Terrestrial NPP: Toward a consistent data set for global model evaluation. Ecol Appl 9: 913–919.

58. Luyssaert S, Inglima I, Jung M, Richardson AD, Reichstein M, Papale D, et al. (2007) CO₂ balance of boreal, temperate, and tropical forests derived from a global database. Global Change Biol 13: 2509–2537. doi: 10.1111/j.1365-2486.2007.01439.x

59. Eurostat (2011) Handbook for annual crop statistics (Following provisions of Reg. 543/2009).

60. Paruelo JM, Epstein HE, Lauenroth WK, Burke IC (1997) ANPP estimates from NDVI for the Central Grassland Region of the United States. Ecology 78: 953–958.

61. Tucker CJ, Pinzon JE, Brown ME, Slayback DA, Pak EW, Mahoney R, et al. (2005) An extended AVHRR 8-km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. Int J Remote Sens 26: 4485–98

62. Wint GRW, Robinson TP (2007) Gridded livestock of the world 2007. Rome, Italy: Food and Agriculture Organization of the United Nations (FAO). 131 pp.

63. Jung M, Vetter M, Herold M, Churkina G, Reichstein M, Zaehle S, et al. (2007) Uncertainties of modeling gross primary productivity over Europe: a systematic study on the effects of using different drivers and terrestrial biosphere models. Global Biogeochem Cy 21, GB4021, doi: 10.1029/2006GB002915

64. Zhao Y, Ciais P, Peylin P, Viovy N, Longdoz B, Bonnefond JM, et al. (2012) How errors on meteorological variables impact simulated ecosystem fluxes: a case study for six French sites. Biogeosciences 9: 2537–2564, doi: 10.5194/bg-9-2537-2012

65. Zaehle S, Sitch S, Smith B, Hatterman F (2005) Effects of parameter uncertainties on the modeling of terrestrial biosphere dynamics. Global Biogeochem Cy 19: 3020, doi: 10.1029/2004GB002395

66. Kramer K, Leinonen I, Bartelink HH, Berbigier P, Borghetti M, Bernhofer C, et al. (2002) Evaluation of six process-based forest growth models using eddy-covariance measurements of CO₂ and H₂O fluxes at six forest sites in Europe. Global Change Biol 8: 213–230. doi: 10.1046/j.1365-2486.2002.00471.x
67. Moral R, Moreno-Caselles J, Perez-Murcia MD, Perez-Espinosa A, Rufete B, Paredes C (2005) Char- acterization of the organic matter pool in manures. Bioresour Technol 96: 153–158. PMID: 15381211

68. Moorcroft PR (2006) How close are we to a predictive science of the biosphere? Trends Ecol Evol 21: 400–407. PMID: 16815439

69. White MA, Thornton PE, Running SW, Nemani RR (2000) Parameterization and sensitivity analysis of the BIOME-BGC terrestrial ecosystem model: net primary production controls. Earth Interactions 4: 1–85, 2000.

70. McKay M, Beckman R, Conover W (1979) A comparison of three methods of selecting values of input variables in the analysis of output from a computer code. Technometrics 21: 239–245.

71. Campolongo F, Kleijnen J, Andres T (2000) Screening Methods, in Sensitivity Analysis, edited by Saltelli A., Chan K., and Scott E. M., pp. 65–80, John Wiley, Hoboken, N. J.

72. Knorr W, Heimann M (2001) Uncertainties in global terrestrial bio- sphere modeling: 1. A comprehen- sive sensitivity analysis with a new photosynthesis and energy balance scheme. Global Biogeochem Cy 15: 207–225.

73. Knorr W, Kattge J (2005) Inversion of terrestrial ecosystem model parameter values against eddy co- variance measurements by Monte Carlo sampling. Global Change Biol 11: 1333–1351. doi: 10.1111/ j.1365-2486.2005.00977.x

74. Kattge J, Knorr W, Raddatz T, Wirth C (2009) Quantifying photosynthetic capacity and its relationship to leaf nitrogen content for global-scale terrestrial biosphere models. Global Change Biol 15: 976– 991. doi: 10.1111/j.1365-2486.2008.01744.x

75. Peeters A (2012) Past and future of European grasslands. The challenge of the CAP towards 2020. In: Goliński P, Warda M, Stypiński M: Grassland—a European Resource? Grassland Sci Eur 17: 17– 32.

76. Nitrates Directive (91/676/EEC) http://eur-lex.europa.eu/

77. Peel MC, Finlayson BL, McMahon TA (2007) Updated world map of the Köppen-Geiger climate classi- fication, Hydrol Earth Syst Sc 11: 1633–1644. doi: 10.5194/hess-11-1633-2007

78. Guenet B, Moyano FE, Vuichard N, Kirk GJD, Bellamy PH, Zaehle S, et al. (2013) Can we model ob- served soil carbon changes from a dense inventory? A case study over England and Wales using three versions of the ORCHIDEE ecosystem model (AR5, AR5-PRIM and O-CN). Geosci. Model Dev 6: 2153–2163, doi: 10.5194/gmd-6-2153-2013

79. Le Houerou HN, Bingham RL, Skerbek W (1988) Relationship between the variability of primary pro- duction and the variability of annual precipitation in world arid lands. J Arid Environ 15: 1–8

80. Silvertown J, Dodd ME, McConway K, Potts J, Crawley M (1994) Rainfall, biomass variation, and community composition in the Park Grass Experiment. Ecology 75: 2430–2437.

81. Briggs JM, Knapp AK (1995) Interannual variability in primary production in tallgrass prairie: climate, soil moisture, topographic position and fire as determinants of aboveground biomass. Am J Bot 82: 1024–1030.

82. Knapp AK, Briggs JM, Koelliker JK, (2001) Frequency and extent of water limitation to primary produc- tion in a mesic temperate grassland. Ecosystems 4: 19–28.

83. Nippert JB, Knapp AK, Briggs JM (2006) Intra-annual rainfall variability and grassland productivity: can the past predict the future? Plant Ecol 184: 65–74.

84. Harpole WS, Potts DL, Suding KN (2007) Ecosystem responses to water and nitrogen amendment in a California grassland. Global Change Biol 13: 2341–2348. doi: 10.1111/j.1365-2486.2007.01447.x

85. Elser JJ, Bracken MES, Cleland EE, Gruner DS, Harpole WS, Hillebrand H, et al. (2007) Global analy- sis of nitrogen and phosphorus limitation of primary producers in freshwater, marine and terrestrial ecosystems. Ecol Lett 10: 1135–1142. PMID: 17922835

86. Wriedt G, van der Velde M, Aloe A, Bouraoui F (2009) A European irrigation map for spatially distribut- ed agricultural modeling. Agr Water Manag 96: 771–789.

87. Peeters A (2009) Importance, evolution, environmental impact and future challenges of grasslands and grassland-based systems in Europe. Grassland Sci 55: 113–125. doi: 10.1111/j.1744-697X. 2009.00154.x

88. Ordoñez JC, van Bodegom PM, Witte J-PM, Wright IJ, Reich PB, Aerts R (2009) A global study of re- lationships between leaf traits, climate and soil measures of nutrient fertility. Global Ecol Biogeogr 18: 138–149.

89. Joyce CB (2013) Ecological consequences and restoration potential of abandoned wet grasslands. Ecol Eng In Press, Available online 24 June 2013.
90. Van den Pol-van Dasselaar A, de Haan MHA, Philipsen AP (2012) Combining grazing and automatic milking. In: Goliński P, Warda M, Stypiński M: Grassland—a European Resource? Grassland Sci Eur 17: 279–281

91. Henkin Z, Gutman M, Aharon H, Perevolotsky A, Ungar ED, Seligman NG (2005) Suitability of Mediterranean oak woodland for beef herd husbandry. Agr Ecosyst Environ 109: 155–161

92. Casasùs I, Bernués A, Sanz A, Villalba D, Riedel JL, Revilla R (2007) Vegetation dynamics in Mediterranean forest pastures as affected by beef cattle grazing. Agr Ecosyst Environ 121: 365–370.

93. Papanastasis VP, Yiakoulaki MD, Decandia M, Dini-Papanastasi O (2008) Integrating woody species into livestock feeding in the Mediterranean areas of Europe. Anim Feed Sci Tech 140: 1–17.

94. Yassoglou NJ, Kosmas C (2000) Desertification in the Mediterranean Europe. A case in Greece. RALA Report No. 200:27–33.

95. Leadley PW, Niklaus PA, Stocker R, Korner C (1999) A field study of the effects of elevated CO₂ on plant biomass and community structure in a calcareous grassland. Oecologia 118: 39–49. doi: 10.1007/s004420050701 PMID: 20135159

96. Fiala K (2010) Belowground plant biomass of grassland ecosystems and its variation according to ecological factors. Ekológia (Bratislava) 29: 182–206.

97. Káplová M, Edwards KR, Květ J (2011) The effect of nutrient level on plant structure and production in a wet grassland: a field study. Plant Ecol 212: 809–819

98. Hussain MZ, Grünwaldb T, Tenhunen JD, Li YL, Mirzae H, Bernhofer C, et al. (2011) Summer drought influence on CO₂ and water fluxes of extensively managed grassland in Germany. Agr Ecosyst Environment 141: 67–76.

99. Alberti G, Delle Vedove G, Zuliani M, Peressotti A, Castaldi S, Zerbi G (2010) Changes in CO₂ emissions after crop conversion from continuous maize to alfalfa. Agr Ecosyst Environ 136: 139–147.

100. Paustian K, Andrén O, Clarholm M, Hansson AC, Johansson G, Lagerlof J, et al. (1990) Carbon and nitrogen budgets of four agro-ecosystems with annual and perennial crops, with and without N fertilization. J Appl Ecol 27: 60–84.

101. Jenkinson DS, Harkness DD, Vance ED, Adams DE, Harrison AF (1992) Calculating net primary production and annual input of organic matter to soil from the amount and radiocarbon content of soil organic matter. Soil Biol Biochem 24: 295–308.