Impact of droughts on the carbon cycle in European vegetation: a probabilistic risk analysis using six vegetation models

. We analyse how climate change may alter risks posed by droughts to carbon ﬂuxes in European ecosystems. The approach follows a recently proposed framework for risk analysis based on probability theory. In this approach, risk is quantiﬁed as the product of hazard probability and ecosystem vulnerability. The probability of a drought hazard is calculated here from the Standardized Precipitation–Evapotranspiration Index (SPEI). Vulnerability is calculated from the response to drought simulated by process-based vegetation models. We use six different models: three for generic vegetation (JSBACH, LPJmL, ORCHIDEE) and three for spe-ciﬁc ecosystems (Scots pine forests: BASFOR; winter wheat ﬁelds: EPIC; grasslands: PASIM). The periods 1971–2000 and 2071–2100 are compared. Climate data are based on gridded observations and on output from the regional climate model REMO using the SRES A1B scenario. The risk analysis is carried out for ∼ 18 000 grid cells of 0.25 × 0.25 ◦ across Europe. For each grid cell, drought vulnerability and risk are quantiﬁed for ﬁve seasonal variables: net primary and ecosystem productivity (NPP, NEP), heterotrophic respiration (Rh), soil water content and evapotranspiration. In this analysis, climate change leads to increased drought risks for net primary productivity in the Mediterranean area: ﬁve of the models estimate that risk will exceed 15 %. The risks increase mainly because of greater drought probability; ecosystem vulnerability will increase to a lesser extent.


Introduction
In the coming decades, climate change is expected to lead to higher temperatures and to greater temporal variability in precipitation and temperature in many parts of the world (Seneviratne et al., 2006). Such changes will alter the frequency, intensity and duration of droughts, with consequent effects on the productivity and carbon balance of ecosystems (Reichstein et al., 2013). Meanwhile, increasing CO 2 has been suggested to increase plant water use efficiency, and to save soil moisture by reducing stomatal conductance (e.g. Drake et al., 1997;Morgan et al., 2004;Keenan et al., 2013), a mechanism which could partly alleviate the sensitivity of plants and ecosystems to drought. How these changes will affect the vulnerability of vegetation to future drought and the associated risks has not yet been rigorously quantified (Van Oijen et al., 2013b).
Risk analysis is applied in many fields, and terminology is inconsistent (Brooks, 2003;Schneiderbauer and Ehrlich, 2004). Generally, "hazard", "vulnerability" and "risk" are distinguished (DHA, 1992). "Hazard" refers to damaging conditions which may on occasion prevail. "Vulnerability" refers to the sensitivity of the impacted system to those damaging conditions. "Risk" is commonly defined as an expectation value for loss induced by the hazardous conditions. Although most work on risk analysis distinguishes three such terms, the definitions -in particular of vulnerability -are generally too imprecise to place the terms in a mathematical relationship to each other (Ionescu et al., 2009). Risk analyses thus tend not to show how exactly risk should be decomposed into its two constituent terms, hazard and vulnerability. This hampers the use of risk analysis to ecosystems when we are interested not only in quantifying risk, or changes in risk, but also in identifying the main causal factors. In the present study, we aim to disentangle two effects of climate change in a risk analysis framework. First, to what extent does climate change lead to increased probability of extreme, hazardous weather? Second, how does climate change affect the vulnerability of ecosystems to such extreme conditions? For this, we require a quantitative method that allows risk to be decomposed.
Recently, we proposed a method for probabilistic risk analysis (PRA) which is able to decompose the risk into two constituent terms (Van Oijen et al., 2013b). The PRA is designed to analyse the effect of any environmental variable on any system variable of interest. An example could be the impact of drought on ecosystem carbon flux. The method allows for either variable to be derived from measurement or from modelling. "Hazardous conditions" are defined as those where the environmental variable is more extreme than a given threshold, and their probability of occurrence is denoted as P (H ). We define vulnerability (V ) as the difference between the expectation values for the system variable under non-hazardous and hazardous conditions. Risk (R) is defined, as is commonly done, as the expectation of loss: the difference between the actual average of the system variable and its value under continuously non-hazardous conditions. With these precise definitions, we achieve the desired decomposition: R = P (H ) × V . All three terms are expected to change over time. Climate change, for example, directly affects the probability of hazardous conditions, and changes in ecosystems alter their vulnerability.
We note that risk decomposition as defined in our PRA framework is not equivalent to "fault tree analysis" (FTA) as practised in engineering projects (Rausand, 2011). FTA quantifies the failure probability of the different components in a human-made system and this is modelled using discrete probability distributions. This approach is not suitable for our purposes in ecology where response variables such as carbon fluxes are not binary: fluxes do not necessarily "fail" when there is a drought but can change to any given degree. Therefore our framework for risk analysis uses continuous probability distributions and we define vulnerability as a function of expectation values and not discrete probabilities. We have not been able to find the PRA equations (given below in Sect. 2.5.1 and more fully in Van Oijen et al., 2013b) in the engineering literature, although of course there is conceptual similarity between the fields.
We apply the PRA framework here to the impact of droughts on vegetation in Europe, both in recent  and future (2071-2100) years. Our focus is on the responses of the major carbon fluxes between ecosystems and the atmosphere: net primary productivity (NPP), heterotrophic respiration (Rh) and net ecosystem productivity (NEP). The different carbon fluxes in ecosystems are tightly linked to each other (e.g. Raich and Schlesinger, 1992;van Oijen et al., 2010), and we can therefore analyse how drought risks to carbon sequestration in the form of NEP result from risks to the other two carbon fluxes, given that NEP = NPP − Rh. In addition, we examine the primary responses of soil water content (SWC) and evapotranspiration (ET), giving five system variables in total. The water variables SWC and ET are included to determine whether risks to NPP are mainly due to reductions in water availability or in water use. We quantify how risks to these five variables under future conditions differ from the present across Europe, and proceed by inspecting the two constituent terms of risk to pinpoint the main cause of change.
We define drought using the Standardized Precipitation-Evapotranspiration Index (SPEI, Vicente-Serrano et al., 2010). The SPEI is a localized measure of drought: it is based for every site on the local long-term frequency distribution of precipitation minus potential evapotranspiration. The SPEI quantifies how rare a given difference of precipitation and potential evapotranspiration is, with respect to this frequency distribution. Because one aim of our analysis is to determine in which parts of Europe climate change is expected to lead to a change in the probability of drought (P (H )), we use the same reference period for the SPEI calculations for both the present and future risk analysis, choosing 1901-2010 for this M. Van Oijen et al.: Impact of droughts on the carbon cycle in European vegetation 6359 purpose. We define as droughts those SPEI values, calculated for the half-year period from April to September, which are less than −1, i.e. more than one normalized standard deviation below average (Sect. 2.2).
Apart from the SPEI, various other drought indices exist (Joetzjer et al., 2013). For example, the Standardized Precipitation Index (SPI) is a simpler index, requiring only information about the time series of precipitation at a site. In contrast to this, the often used Palmer Drought Severity Index (PDSI) requires extra information about soil water holding capacity to allow estimation of a local water balance. The SPEI does not require such detailed information, yet accounts for the impact of other meteorological variables besides precipitation on the severity of a drought. There is longstanding and active debate in the literature about the strengths and weaknesses of different drought indices (e.g. Vicente-Serrano et al., 2010;Sheffield et al., 2012;Joetzjer et al., 2013). Much of this debate centres on whether the indices can detect whether ongoing climate change has already altered drought prevalence (e.g. Trenberth et al., 2014). In contrast, our main concern here is the comparison of the impacts of recent droughts with those from expected droughts a century into the future. This period is expected to see greater changes in weather than the past century, making the choice of drought index less critical. This paper quantifies current and future vulnerability and risk of change in major European carbon fluxes in response to drought. We present results from the PRA using six different vegetation models and discuss how the results depend on the way drought impacts are simulated. Three of the models simulate a combination of vegetation types in each grid cell, and three are ecosystem-specific (Scots pine, winter wheat, grassland). Comparison of the results from the different models allows a modest assessment of the uncertainty of the risk analysis with respect to model assumptions. We evaluate the PRA using remotely sensed Normalized Difference Vegetation Index (NDVI) data, and carry out a sensitivity analysis to determine the robustness of the results with regard to timing and severity of drought as well as climatic variability.

Materials and methods
In this section, we describe the atmospheric drivers (Sect. 2.1), the SPEI drought index (Sect. 2.2), the six models (Sect. 2.3), the Normalized Difference Vegetation Index (NDVI) vegetation greenness data (Sect. 2.4) and our new method for probabilistic risk analysis (PRA, Sect. 2.5).

Climate and CO 2
The vegetation models used in this study are driven with daily climate data spanning the period 1901-2100 at a 0.25 × 0.25 • longitude/latitude grid. The climate data are described in detail by Beer et al. (2014). Here, a short sum-mary follows. For the period 1901-1978 the data set is based on WATCH forcing data (Weedon et al., 2011). For the period 1979-2010, ECMWF ERA-Interim reanalysis data are used, bias-corrected against WATCH observation-derived climate forcing data using the method described by Piani et al. (2010). From 2011 until 2100, climate data are used on the basis of outputs from the regional climate model REMO using the SRES A1B scenario, with boundary conditions from the ECHAM5 general circulation model, as provided by the EU project ENSEMBLES (Hewitt and Griggs, 2004). This data set has also been bias-corrected against the 1901-2010 WATCH data set following Piani et al. (2010).
We refer to the above mentioned climate data for the period 1901-2100 as the "control climate" and most of our study will focus on the implications of climate change inherent in that data set. Besides these control climate runs, we conducted model experiments using an additional artificial climate data set (Beer et al., 2014), referred to as the "reduced variability climate", which was constructed such that day-to-day and year-to-year variability of the meteorological variables are smaller than in the control while the 30-year long-term averages of annual and seasonal values are being conserved (Beer et al., 2014).
The time series for the atmospheric concentration of CO 2 is based on data from ice-core records and NOAA atmospheric observations for the period 1901-2010 and on the SRES (Nakicenovic and Swart, 2000) A1B scenario for 2011-2100. The concentration increases from 296 ppm in 1901 to 710 ppm in 2100.

Standardized Precipitation-Evapotranspiration Index (SPEI)
To calculate time series of drought in each grid cell, we use the Standardized Precipitation-Evapotranspiration Index (SPEI, Vicente-Serrano et al., 2010). This index quantifies the degree of drought as a function of the difference D between precipitation and potential evapotranspiration. The SPEI is a normalization of this difference, defined as the normally distributed variable (mean = 0, standard deviation = 1) whose cumulative probability density function coincides with that of D. The probability distribution for D is estimated from monthly data in a reference period, which we chose to be the period for which we had observational weather data, 1901-2010. With these choices, for every grid cell the average SPEI value over the years 1901-2010 is zero and about 68 % of SPEI values are between −1 and +1. Using a long reference period, here 110 years, to calculate SPEI values is advisable when the focus of the study is on rare extreme events. The SPEI can be calculated for drought events of any duration, but here we only used values calculated for periods of half a year. In most cases, the summer half of each year was used, from April to September, which in Europe captures the season of highest vegetation greenness across all latitudes (see Sect. 2.4). As a  Table 1. Average values of model outputs for April-September, averaged over 1971-2000. NPP is net primary productivity (g C m −2 d −1 ), Rh is heterotrophic respiration (g C m −2 d −1 ), NEP is net ecosystem productivity (g C m −2 d −1 ), SWC is soil water content (mm), ET is evapotranspiration (mm d −1 ), and n is the number of grid cells simulated.

Model
Vegetation sensitivityanalysis, we assessed the consequences of shifting the start of the half-year to March, February or January. The Penman-Monteith equation (Monteith, 1995) with invariant canopy conductance was used to calculate potential evapotranspiration. Figure 1 shows maps for the average value of the April-September SPEI in the years 1971-2000 and 2071-2100 plus an example for the European heat wave year 2003. The two 30-year time periods are the periods for which we carry out risk analysis in this paper. The SPEI in the year 2003 is shown to illustrate how a single extreme year can differ from the average. 2003 had a very wet and cold spring in north Eastern Europe, but a very hot and dry summer in Western and Central Europe Vetter et al., 2008). The years 1971-2000 are all within the reference period for SPEI calculation , so the mean values of SPEI shown in Fig. 1a are close to zero throughout Europe. A century later, northern Europe is expected to have become even wetter and the south even drier (Fig. 1c). The future persistent spring and summer drought in the south, with average SPEI values −2 to −3, is comparable to the lowest values observed over temperate western regions in the currently still exceptional year 2003 (Fig. 1b).

Models
Six different process-based vegetation models are used in this study. These include three ecosystem-specific models: BAS-FOR for Scots pine forests, EPIC as applied to winter wheat fields and PASIM for grassland. These models are referred to as "ecosystem models". The other three models simulate the different types of vegetation present in each grid cell using a functional-type approach: JSBACH, LPJmL and OR-CHIDEE. These are called "generic models".
All models are run using the same climatic input and the same 0.25 × 0.25 • spatial resolution. Only the generic models are applied to all grid cells (Table 1)

BASFOR (Scots pine)
BASFOR ) is a forest model that runs with a daily time step and requires as environmental inputs: weather, nitrogen deposition, CO 2 concentration, soil carbon and nitrogen content, and soil water retention characteristics. For the present risk analysis, we use the model as parameterized for Scots pine based on forest inventory data from four different European countries (Van Oijen et al., 2013a). Information on soils and on atmospheric nitrogen deposition is based on work by Cameron et al. (2013). The model is run with a fixed forest rotation length of 80 years with regularly spaced thinnings, for eight age-cohorts (Cameron et al., 2013). Trees respond to drought by reducing carbon uptake, accelerating senescence and changing allocation in favour of the roots, but mortality is not simulated. The model has been used for the analysis of carbon and nitrogen fluxes in Norway spruce forest (Van Oijen et al., 2011) and for assessing the impact of climate change on carbon sequestration in forests in the United Kingdom (Van Oijen and Thomson, 2010).

EPIC (winter wheat)
EPIC (Environmental Policy Integrated Climate; Williams et al., 1989) is an agro-ecosystem model running with a daily time step. It simulates crop development and yield, hydrological, nutrient and carbon cycles and a wide range of crop management activities. Required input data include minimum and maximum temperature ( • C), precipitation (mm), global radiation (MJ m −2 ), CO 2 concentration, soil properties (bulk density (g cm −3 ), base saturation (%), percentages of soil organic carbon, sand, silt and clay), crop heat requirement and crop management (planting dates, tillage operations, fertilizer application). In the present study, winter wheat is simulated in 1 year mono-crop rotations on all the available European cropland. It is simulated as a winter crop, with varieties based on growing season length and heat requirements distinguished for Atlantic, alpine, boreal, continental and Mediterranean climatic zones. The simulations are performed with current fertilizer input levels. A detailed description and evaluation of this implementation is provided by Balkovič et al. (2013).
Plant growth is simulated using the concept of radiationuse efficiency (Monteith, 1977). Potential biomass increase is calculated as a function of photosynthetically active radiation and leaf area index. CO 2 fertilization increases radiationuse efficiency and stomatal resistance, thereby reducing transpiration (Stockle et al., 1992). Crop yield is calculated via a harvest index, determining yield as a fraction of aboveground biomass. Daily heat unit accumulation determines leaf area growth and senescence, canopy height, nutrient uptake, harvest index and date of harvest. Temperature also determines potential ET, which is calculated with the Hargreaves method. Actual ET is governed by leaf area index and the water content and depth of the root zone. Both actual ET and biomass increase are reduced when the plantavailable soil water level is less than 25 % of the maximum. Heterotrophic respiration is also temperature and moisture dependent following the carbon and nitrogen cycle model of Izaurralde et al. (2006).
Previously we have evaluated and used the EPIC model for hindcasting assessments of the impact of extreme dry and wet weather on crop production (Van der Velde et al., 2010, 2012).

PASIM (grassland)
PASIM, the Pasture Simulation model (Riedo et al., 1998;Vuichard et al., 2007) is a multi-year biogeochemical model that simulates water, carbon and nitrogen cycles in grassland systems at the plot scale on a daily to sub-daily time step. Soil processes are based on the CENTURY model of Parton et al. (1988). Photosynthetically assimilated carbon is either respired or allocated dynamically to one root compartment and three shoot compartments. Accumulated aboveground biomass is removed by either cutting or grazing, or enters a litter pool. Soil organic carbon (SOC) is represented in three pools (active, slow and passive) with different potential decomposition rates, while above and below-ground plant residues and organic excreta are partitioned into structural and metabolic pools. The nitrogen cycle considers three types of nitrogen inputs to the soil via atmospheric nitrogen deposition, fertilizer nitrogen addition, and symbiotic nitrogen fixation by legumes. The inorganic soil nitrogen is available for root uptake and may be lost through leaching, ammonia volatilization and nitrification/denitrification, the latter processes leading to nitrous oxide (N 2 O) gas emissions to the atmosphere. Management includes nitrogen fertilization, mowing and grazing and can either be set by the user or optimized by the model (Vuichard et al., 2007;Graux et al., 2011). In this model version, nitrogen fixation is simulated by assuming a constant legume fraction, an algebraic method for SOC equilibrium search (Lardy et al., 2011) is used to reduce computation time, and grassland nitrogen fertilization corresponds to current farming practices (Weiss and Leip, 2012). Stomatal conductance depends on soil water content and atmospheric vapour pressure and CO 2 concentration. Drought-induced stomatal closure leads to a decreased photosynthetic rate and increases in vegetation temperature as well as autotrophic and heterotrophic respiration.

JSBACH (generic vegetation)
The land surface scheme JSBACH (Raddatz et al., 2007) simulates land-atmosphere exchanges of energy, water and carbon at 30 min temporal resolution. The representation of canopy processes is based on the BETHY model which couples energy and water balance with photosynthesis through stomatal conductance (Knorr, 2000) to which a module for phenology and a simple carbon cycle scheme including several pools for vegetation, litter and soil carbon have been added (Raddatz et al., 2007). Canopy conductance is additionally constrained by soil water availability (Raddatz et al., 2007). In the version used for this study, soil hydrology is represented by a five layer scheme (Hagemann, 2014). The fractional coverage of land by several plant functional types is prescribed based on GLC2000 (Bartholomé and Belward, 2005) and MODIS VCF land cover information (Hansen et al., 2003). Disturbances are not explicitly simulated but taken into account by the mean residence time of carbon pools (Raddatz et al., 2007).

LPJmL (generic vegetation)
LPJmL (Sitch et al., 2003;Bondeau et al., 2007;Gerten et al., 2004) is a dynamic global vegetation model that simulates process-based terrestrial vegetation dynamics and physiology and land-atmosphere carbon and water fluxes. Its capability to adequately represent ecosystem responses to changing precipitation has been demonstrated by Gerten et al. (2008). Vegetation is represented as 10 plant functional types (PFTs) which differ in their bioclimatic limits and morphological, phenological and physiological parameters. Bioclimatic limits determine the ability of a species to establish under a given climate, ecophysiological parameters determine the ability to survive and to become dominant species. Each plant functional type is represented by an average individual covering a fraction of the grid cell. For each average individual of the 10 PFTs, LPJmL simulates photosynthesis and respiration, the allocation of accumulated carbon to the plant's compartments (leaves, stem, roots, reproductive organs). Vegetation structure and dynamics are updated on an annual time step, while photosynthesis, evapotranspiration and soil water exchange are calculated daily. The soil consists of two layers with water content driven by precipitation, snowmelt, percolation, runoff and evapotranspiration. PFTs differ in their fraction of roots in the upper and lower soil layer representing shallow-and deep-rooting plants.
Evapotranspiration is calculated for each PFT as a function of soil moisture supply and atmospheric demand. If atmospheric demand is higher than water supply, canopy conductance is reduced until transpiration equals the supply. This directly influences photosynthesis rates by reduced diffusion of CO 2 into the leaf intercellular space. Canopy photosynthesis and thus gross primary productivity are calculated for each PFT based on a modified Farquhar photosynthesis model (Farquhar et al., 1980;Collatz et al., 1991) driven by absorbed photosynthetically active radiation, temperature and leaf intercellular CO 2 concentration under the assumption of optimal nitrogen availability. Net primary productivity (NPP) is the difference between gross primary productivity GPP and carbon lost from growth and maintenance respiration. Net ecosystem productivity is derived from the difference of NPP and heterotrophic respiration (Rh). Heterotrophic respiration results from the decomposition of litter and soil organic matter and is temperature and moisture dependent (for details see Sitch et al., 2003).

ORCHIDEE (generic vegetation)
ORCHIDEE  is a land surface model that simulates carbon, water and energy exchanges within ecosystems, and between the land and the atmosphere. The land vegetation is represented by 12 plant functional types (PFTs). The definition of a PFT accounts for physiognomy (tree or grass), leaf types (needle-leaf or broad-leaf), phenology (evergreen, summer-green or rain-green), the photosynthetic pathway for crops and grasses (C 3 and C 4 ) and the climatic regions (boreal, temperate and tropical). Exchanges of energy and water are calculated following Ducoudré et al. (1993) and de Rosnay and Polcher (1998) at 30 min time steps.
Potential transpiration is limited by aerodynamic resistance above and within the canopy and stomatal resistance (Ducoudré et al., 1993), and potential evaporation is proportional to the humidity difference between air and air at saturation for soil temperature. Bare soil evaporation is the maximum upward hydrological flux permitted by diffusion if this flux is smaller than potential evaporation (d'Orgeval et al., 2008). When root zone soil moisture falls below 40 %, stomatal conductance, photosynthesis and evapotranspiration decrease linearly, reaching zero for entirely dry soils (McMurtrie et al., 1990). Atmospheric dryness also affects these processes, following Ball et al. (1987). The 11 layer scheme for the upper 2 m of soil is based on the CWRR model (Bruen, 1997). Layers are thinner near the surface (De Rosnay et al., 2002). Vertical water flow is based on the one-dimensional Darcy equation (De Rosnay et al., 2002). The bottom boundary condition is set to be free drainage.
GPP is calculated every 30 min using the Farquhar et al. (1980) and Collatz et al. (1992) leaf scale equations for C 3 and C 4 PFTs respectively. Net primary productivity (NPP), plant growth, maintenance respiration, mortality, litter and soil organic matter decomposition are calculated on a daily time step following GPP. Water stress alters carbon allocation and accelerates leaf senescence . The model accounts for fires but not for vegetation dynamics. Drought has lagged effects on leaf area index (LAI) through reduced production of carbon reserves for the following year, but does not directly increase tree mortality. Heterotrophic respiration (Rh) is coupled positively to NPP on time scales of weeks to months via the short-lived litter and labile soil carbon pools. Decomposition increases with warmer temperatures, according to a Q 10 law, and with relative soil moisture .

Normalized Difference Vegetation Index (NDVI)
The Normalized Difference Vegetation Index (NDVI, Tucker, 1979) is a radiometric measure of vegetation photosynthetic capacity. It compares the reflectance difference between photosynthetically active radiation (PAR, ∼ 400-700 nm) and near-infrared radiation (NIR, ∼ 700-1100 nm). The index is indicative of the abundance of chlorophyll content in the vegetation canopy (Myneni et al., 1992), and is highly correlated with gross primary productivity (Veroustraete et al., 1996). The index is derived from satellite observations, made globally over long time periods. It is used in the present study for two purposes: to determine the timing of the growing season at different latitudes in Europe and to provide independent data for testing the performance of the six models including the PRA.
We use the latest version of the Global Inventory Modeling and Mapping Studies (GIMMS) NDVI data set (NDVI3g) generated from the Advanced Very High Resolution Radiometer (AVHRR). The AVHRR instruments aboard a series of NOAA polar-orbiting meteorological satellites (NOAA 7,9,11,14,(16)(17)(18)(19) have observations in multiple spectral bands from 1981 to the present. As the first set of instruments that have non-overlapping spectral bands from visible (VIS) and near-infrared (NIR) ranges, NDVI can be calculated as (NIR−VIS) / (NIR+VIS). The GIMMS data set has been recalibrated recently to improve data quality at northern high latitudes, and an overall uncertainty of ±0.005 NDVI units is achieved independent of time frame (Pinzon and Tucker, 2014). It has also been corrected for calibration, viewing geometry, volcanic aerosols, and other effects not related to vegetation change (Pinzon and Tucker, 2014;Pinzon et al., 2005;Tucker et al., 2005).
The GIMMS NDVI3g data set contains NDVI values from 1981 to 2012 with 15 day temporal and 8 km spatial resolution, of which the years 1982-2010 are used in the present study. The data set has been used to explain diminishing vegetation seasonality in the northern high latitudes over the past 30 years (Xu et al., 2013), vegetation greening and browning trends (De Jong et al., 2013), latitudinally asymmetric vegetation growth trends with small increases in the Southern Hemisphere and larger increases at high latitudes in the Northern Hemisphere (Mao et al., 2013), global trends of seasonality change in natural areas (Eastman et al., 2013) and an earlier start of the growing season in Fennoscandia (Høgda et al., 2013). Figure 3a shows a map of the average value of the NDVI across Europe in the years 1982-2010, months April-September. The lowest NDVI values are in the cold north and the dry Mediterranean. The same map further shows the division in three latitudinal bands ("North", "Mid", "South") that we use for summarizing results later in the study. The unbroken lines in Fig. 3b show the intraannual progression of NDVI averaged over the years per latitudinal band. Seasonal peak NDVI is reached by May in the south, a month later at mid-latitude and yet another month later in the north. Based on these observations, we used a common time interval of April to September for the calculation of the SPEI at all latitudes. This period is close to the period of April-August in which different kinds of extreme events, but predominantly water scarcity, were shown to reduce gross primary productivity in Europe the most in the years 1982-2011(Zscheischler et al., 2014 Fig. 8). To indicate inter-annual variation, the five dashed lines in Fig. 3b represent the individual years 2001, 2002, 2004 and 2005, with 2003 in bold. The variation between years is much lower than between months, even when, as here, an extreme year like 2003 is included.

Theory for PRA
Following the procedure outlined in greater detail by Van Oijen et al. (2013b), we distinguish environmental and system variables, denoted as env and sys. Sys variables are expected to be sub-optimal when env variables reach levels that constitute hazardous conditions. Numerically, risk (R) is defined as the expectation of system loss, i.e. the amount by which average system performance is less than it would be under continuously non-hazardous conditions: where E(sys) is the overall expectation value of the sys variable, and E(sys|env non-hazardous) is the expectation value of the sys variable when conditions are not hazardous. The probability of env variables assuming hazardous values is denoted as P (hazardous), or P (H ) in short. Vulnerability (V ) has, in contrast to most previous applications of risk analysis, a precise quantitative definition: with the obvious interpretation of terms. This definition of vulnerability as the difference in sys performance between non-hazardous ("good") and hazardous ("bad") conditions allows us to write R alternatively as: Equations (1) and (3) are mathematically equivalent, giving the same estimates of risk. The last equation shows that the risk to sys variables can be decomposed into two terms: the probability of hazardous environmental conditions and the vulnerability of the system.

Application of PRA in this study
The env variable in the present study is the SPEI, and hazardous conditions are mostly defined as SPEI being less than −1 (Sect. 2.2). We carry out the PRA for five different sys variables: net primary productivity (NPP), heterotrophic respiration (Rh), net ecosystem productivity (NEP), soil water content (SWC) and evapotranspiration (ET). Values for the sys variables are generated by the six different vegetation models. PRAs are carried out for each grid cell separately, and for two time periods: 1971-2000 and 2071-2100. P (H ) is calculated as the fraction of the thirty years in each period with SPEI lower than the threshold. Likewise, the various expectation values E(sys|env) are calculated from the frequency distribution of sys values over the thirty years. The main outcomes of the PRA are vulnerabilities and risks, typically expressed in absolute values (e.g. g C m −2 d −1 for the carbon fluxes). In the last two figures of the paper, we express them as fractions of system performance in non-hazardous years.
An additional PRA is carried out using NDVI as the sys variable. Drought vulnerability and risk for NDVI are not of interest themselves, but comparison of V and R for NDVI with those for model outputs provides a coarse empirical check of the accuracy of the model-based analysis.

Sensitivity analysis of the PRA
After the main PRAs outlined above, we subject NPP output from one of the generic models, LPJmL, to additional PRAs to assess the robustness of our results. First, we assess how vulnerability and risk estimates change when a different SPEI threshold is chosen to demarcate hazardous and non-hazardous conditions. We vary the SPEI threshold over the sequence −2, −1.5, −1, −0.5, 0. We then test the idea that a lag time in the response of sys to env should be accounted for. The half-yearly periods for calculating NPP are kept at April-September, but the half-yearly SPEI periods are shifted to be one, two or three months earlier. In a third and final sensitivity analysis, we vary the climate scenario by reducing inter-annual weather variability, while conserving 30 year long-term averages of annual and seasonal values (see Sect. 2.1).

The basic PRAs
Each of the six models has been used to calculate PRAs for each grid cell (the number of which varied between models, Table 1). Per grid cell, ten PRAs have been carried out by each model to assess five variables (NPP, Rh, NEP, SWC, ET) in two time periods . These PRAs represent the core results of this study.
To clarify the approach, we show in Fig. 4 all the details of one example: the PRA for NPP in 1971-2000 according to model LPJmL. The different steps in the PRA are represented by the different panels of the figure. Panel A shows the average NPP in the non-hazardous years, i.e. the years with SPEI for months 4-9 above −1. Panel B shows the same for the hazardous years. Vulnerability, shown in the next panel, is calculated as the difference between the two: C = A-B. The frequency of years with SPEI below −1, i.e. the probability of hazardous conditions, P (H ), is depicted in panel D. Finally, risk, shown in the last panel, is the product of vulnerability and P (H ), so E = C × D. In this example, vulnerability is largest in the south, but P (H ) varies little across Europe. According to this particular vegetation model, the risk to NPP posed by drought (SPEI < −1) is therefore highest in the Mediterranean because the vegetation is more vulnerable to exceptional droughts there than elsewhere, and not because of a higher drought frequency. But note that SPEI < −1 represents far drier conditions in the Mediterranean than further north. Figure 5 summarizes the risk analyses for NPP per latitudinal band for all six models, and for both time periods. The bars show relative vulnerability, i.e. vulnerability divided by the NPP at the same location under contemporaneous nonhazardous conditions. The black part of each bar is the risk, also in relative terms. The ratio of the black part and the whole bar, i.e. risk divided by vulnerability, is equal to the hazard probability, P (H ) (Eq. 3). Most models indicate that relative vulnerabilities and risks increase from north to south, in both time periods. Note that in some cases, mainly found in the north, the values are negative. These are conditions where dry years (SPEI < −1) are in fact the more favourable for NPP, likely because of climatic variables that co-vary with SPEI. During dry years in the period 1971-2000, radiation in the north was on average reduced by only 0.1 %, which is unlikely to be of any significance, but temperature was 0.2 • C higher, possibly extending growing seasons and stimulating NPP.
Overall the relative risks to NPP are higher in 2071-2100 than in the earlier period, especially in the south according to all six models. The main reason is an increase in hazard probability rather than an increase in vulnerability.
The three ecosystem models show very similar results, especially in the most sensitive region, the south. This is despite the models simulating different vegetation types and partly different sets of grid cells. In contrast, the three generic models differ strongly from each other, with vulnerabilities and risks to NPP decreasing in the order LPJmL > ORCHIDEE > JSBACH.
The results of the PRAs for all five model output variables and both time periods are summarized in Table 2. Note that the values in the table are in absolute units, in contrast to Fig. 5. The highest and lowest values of vulnerability and risk, across all models, are indicated for each variable in bold italic typeface.
The vulnerabilities and risks for Rh differ strongly between models. In over half the cases, the ecosystem models show negative vulnerability and risk for Rh (i.e. drought increases soil respiration). The vulnerabilities found by the generic models are always positive (i.e. drought reduces respiration), likely because in those models, Rh anomalies are tightly coupled with NPP trough fast carbon pools. NEP can be calculated as NPP minus Rh, and likewise NEP vulnerability and risk can be calculated from those for NPP and Rh. Therefore the results for NEP vary strongly between the models as well, with vulnerabilities that are always less than for NPP in the case of the generic models, but often the opposite for the ecosystem models.
The vulnerability and risk for SWC also vary considerably among the models, with lowest values estimated by the two generic models that have low average SWC (Table 1): LPJmL and ORCHIDEE. But all models agree that SWC vulnerabilities do not differ much between latitudinal bands or time periods. Climate change does induce a spatial gradient to the risk, with highest values in the south. Vulnerability and risk for ET show greater spread among the ecosystem models (with PASIM simulating relatively low values in both time periods) than among the generic models. Climate change increases risks mainly in the south.

Sensitivity analysis of the PRA
Results of the assessment of the robustness of the PRA to various assumptions are shown in Fig. 6. The figure shows relative values of vulnerability and risk for NPP as simulated by model LPJmL for both time periods. The left column shows how results depend on the choice of the environmental threshold, which is varied from −2 to 0 by half-unit steps. The bars in the middle are for the default SPEI threshold of −1, thus showing the same results as in Fig. 5. When the threshold is increased, a greater fraction of years is classified as hazardous, so risk increases relative to vulnerability. Vulnerability itself decreases, mainly in the first time period. In the second time period, vulnerability is nearly independent of the choice of the threshold.
In the second sensitivity analysis, shown in the middle column of Fig. 6, the timing of the half-year long SPEI period has been shifted by different amounts: 0, −1, −2, −3 months. So in this case the rightmost bars, with shift equal to zero months, repeat the results of Fig. 5. Overall, vulnerability and risk are highest for the zero shift, i.e. drought periods that cover the same six months (April-September) as are chosen for NPP.
In the final sensitivity analysis, shown in the right column, the models have been run with two climate scenarios: the control scenario used in all other assessments, and the reduced variability scenario (Sect. 2.1). The results are very robust against this change, for both time periods.

Comparisons with NDVI
To test the performance of the models, we examine how the spatial patterns over Europe simulated by the models correlate with the spatial pattern of NDVI observations. Table 3 shows the correlations across all grid cells of model outputs with NDVI. With the single exception of evapotranspiration as simulated by PASIM, the correlations are positive, and strongly so for the fluxes simulated by the generic models (r > 0.6).
A second test is aimed at evaluating not the original model outputs but the vulnerabilities and risks derived from them. To this end, we have applied the risk analysis to NDVI and report in Table 4 the correlations between V and R for NDVI Table 2. Drought vulnerability and risk of the carbon and water balance across Europe. Averages for two time periods and three latitudinal bands, calculated by six models. N is above 55 • N, M is between 45 and 55 • N, S is below 45 • N. NPP is net primary productivity (g C m −2 d −1 ), Rh is heterotrophic respiration (g C m −2 d −1 ), NEP is net ecosystem productivity (g C m −2 d −1 ), SW is soil water content (mm), ET is evapotranspiration (mm d −1 ). A bold font is used for the most extreme vulnerabilities and risks for each variable.

Model
Time   with those for NPP and ET. We have restricted this test to NPP and ET because changes in NDVI reflect changes in leaf area which has direct effects only on the fluxes of carbon and water through the plants. With the exception of ET risk calculated by PASIM, all correlations are again positive, but generally less strong than for the original sys variables (Table 4).

Theory and current application of probabilistic risk analysis
The method of PRA that we have applied here was proposed in a recent publication (Van Oijen et al., 2013b). We argued there that the strengths of the method included its quantitative nature, its simplicity, the possibility of choosing any env variables and sys variables of interest, as well as the Figure 5. Drought risk analysis for NPP in two time periods , for the six different models. Top row: above 55 • N; mid row: between 45 and 55 • N; bottom row: below 45 • N (see Fig. 3a). The whole bars depict vulnerability (V ), the black parts are risks (R). V and R are expressed in relative terms, as fractions of NPP under non-hazardous conditions. Drought probabilities, P (H ), are equal to R/V . Table 3. Correlations of NDVI with model outputs at the European scale (where n is model-specific, see Table 1). NDVI and outputs are averages for April-September over multiple years (NDVI: 1982(NDVI: -2010(NDVI: , outputs: 1971(NDVI: -2000. Variables as in Table 2 applicability to any environmental thresholds and any time lag between env and sys. These strengths seem to be borne out by the comprehensive application of the method in the present study. The main advantage of the method may be that it allows decomposition of risk into two multiplicative terms: the probability of hazardous conditions and the vulnerability of the system. Our analysis focused on the half-year period of April to September, which was estimated from NDVI observations Table 4. Correlations of vulnerabilities and risks calculated for NDVI with those calculated for model outputs of NPP and of ET at the European scale (where n is model-specific, see Table 1). All values are averages for April-September over multiple years (NDVI: 1982(NDVI: -2010(NDVI: , outputs: 1971(NDVI: -2000.

MODEL
Vuln ( Fig. 3b) as the main period of high carbon uptake and water use by vegetation in the different European climatic zones. Climate change toward drier summers may affect the annual pattern of vegetation activity, but shifts to early spring or late autumn will be constrained by the lower levels of solar radiation. The risk analysis could be refined by using slightly different time periods for different parts of Europe, more closely matching the growing season at each location. NDVI is a convenient proxy for identifying this period for large regions, given its global availability and high spatial resolution. Moreover, the high spatiotemporal coverage of NDVI allowed us to test our models and risk analyses in various ways as shown in Tables 3 and 4, and discussed below.
In contrast to NDVI, use of the comparatively sparse set of European eddy covariance towers, which also have footprints considerably smaller than the areas of grid cells, would not have afforded the same capability. Because we focused on half-year periods, temporary interruptions of ecosystem activity by droughts of shorter duration are not included in our analyses, especially if vegetation in models recovered quickly after the drought ended. But such short-term droughts would then not play a significant role in the long-term carbon balance of the vegetation. Similarly, use of longer periods than the growing season, e.g. whole years, would have caused dilution of the drought signal by precipitation in the off-season.
The quality of our ecosystem vulnerability assessment, in particular for future conditions, is limited by how well the six different models represent processes of adaptation to environmental change. The models simulate vegetation change (migration, fire disturbance in some models, acclimatization and physiological adaptation) only to a limited degree. Three of the six models are dynamic vegetation models which allow for replacement of plant functional types by others when the environment changes, but migration is not explicitly simulated. However, in Europe, land use change, forest management and landscape fragmentation are likely to limit future ecosystem migration, except in the very few regions where ecosystems are not managed. Physiological adaptation is simulated to some extent by the models (e.g. stomatal closure with increased atmospheric CO 2 concentration and drought, increased allocation to roots when soil resources become limiting, temperature optimum of photosynthesis). It is conceivable that the vegetation types will adapt more strongly to the new climatic conditions than can be foreseen at this stage. However, our main results seem to be robust because additional adaptation processes would be likely to reduce vulnerability, whereas our risk decomposition already identified increased hazard probability as the greater threat (Sect. 3.1). Also, changes in vulnerability would have to be extreme and much more favourable in the Mediterranean area than elsewhere to overturn our further prediction that the southern part of Europe is at the greatest risk. Regarding the drought hazard itself, a consistent drying trend in Southern Europe, with increased drought extremes, is also predicted in the recent work of Jacob et al. (2014).

Vulnerabilities and risks in the years 1971-2000
We found that relative drought vulnerability of NPP in the early period, 1971-2000, was highest in the south (Fig. 5). This result seems fairly robust as it was found by all six models. Numerically there was variation, with LPJmL estimating higher vulnerabilities than the other models, and JS-BACH lower ones. The exceptional vulnerabilities of these two generic models can partly be explained by their widely differing estimates of the average amount of water in the root zone (SWC) of the vegetation (Table 1), but SWC varied strongly for the other models, too. NPP is the difference between GPP and autotrophic respiration, and high drought sensitivity of GPP as simulated by LPJmL was already reported by Zscheischler et al. (2014).
Despite the relatively high NPP vulnerabilities in the south, the risks were commonly estimated as being low (< 10 %) due to the low probability of drought. The very similar responses of the three ecosystem models (vulnerability < 15 %, risk < 5 %) were not expected, given that the models simulate very different ecosystem types, and model structures are very different. This may be the result of compensating factors. Wheat and grass are more shallow-rooting than trees, making them more drought sensitive in the short term, but they can more easily be irrigated or re-sown and suffer less from carry-over effects in the case of consecutive drought years.
Unlike for NPP, vulnerability of Rh did not show a consistent latitudinal trend for most of the models (Table 2). For all latitudes, the ecosystem models estimated lower Rh vulnerabilities than the generic models did. In fact, many estimates by the ecosystem models were negative, indicating that the droughts increased soil respiration. It is likely that this is the consequence of differences between the models in stress response functions, with soil organic matter turnover being mainly regulated by temperature in the ecosystem models, and mainly by water availability in the generic models. The systematic disagreement between the two model types is surprising, and it points to significant uncertainty about the regulation of organic matter turnover. However, despite model differences, risks for Rh were always assessed as small (< 0.35 g C m −2 d −1 ) because of low P (H ).
The estimates for drought vulnerability and risk to NEP were generally highest in the south, but the values varied considerably among models. The drought response of NEP tended to follow that of NPP more closely than that of Rh for most models except JSBACH. JSBACH simulated much greater drought-induced reductions of Rh than of NPP and was therefore the only model to predict increases of NEP under drought even in the south and in both time periods. To quantify the degree to which the vulnerability of NEP was generally more closely related to NPP vulnerability than to Rh vulnerability, we carried out a correlation analysis for each model. The vulnerability correlation coefficients for NEP-NPP ranged from 0.70 to 0.96 across all models, whereas the NEP-Rh correlation coefficients were in the range −0.65 to 0.11. The results for risk were similar, confirming that drought effects on NPP were stronger determinants of changes in net carbon flux than drought effects on Rh.
The models seemed more consistent in their estimates of vulnerability and risk for SWC and ET, but as these variables were mainly quantified to help understand the underlying causes of differences in carbon flux risks, which were small in 1971-2000, we postpone discussion of them to the case of climate change (2071-2100), treated next.

Impact of climate change on vulnerabilities and risks
Climate change, studied here by comparing results for 2071-2100 with those for 1971-2000, affected vulnerabilities and risks across Europe. In Table 2, we used bold italic typeface to highlight the extreme values of vulnerability and risk for each of the five examined sys variables. Invariably, the highest values were found in the period 2071-2100. For carbon fluxes this result can partly be explained by CO 2 sensitivity: with a general increase in carbon fluxes because of elevated CO 2 , differences between wet and dry conditions may be expected to increase as well. Averaged across the models, the increases in NPP, Rh and NEP were 17, 22 and 10 %, respectively. However, vulnerabilities generally increased less than that, and even decreased slightly in mid-Europe and the north, possibly because of elevated CO 2 increasing drought tolerance by inducing stomatal closure.
In contrast, climate change is expected to increase risks strongly, particularly in the south, according to the scenario and models used in this study. The underlying cause is mainly the probability of drought, which increases far more than vulnerabilities do. Risks exceeding 0.25 g C m −2 d −1 are predicted for carbon sequestration (NEP) in the south by most models, amounting to reductions in carbon sequestration of 20 to 80 %. In contrast, the model JSBACH shows different results, even predicting negative risk for NEP at all latitudes. JSBACH differs from the other models in that it predicts a far greater drought vulnerability of Rh than of NPP. Because the model predicts that carbon uptake is reduced less than carbon loss under future droughts, it thus predicts increased carbon sequestration in drought years. This is a counter-intuitive result, which needs to be evaluated more closely in future work.
The southward increasing risks to carbon fluxes that are predicted by all models except JSBACH, are mirrored by increasing risks to SWC and ET. In relative terms (compare Table 2 with the European average values given in Table 1), risks to SWC and ET tend to vary to a similar degree, suggesting that the risk to water use by vegetation is controlled by soil water availability rather than by stomatal regulation or atmospheric feedbacks.

Comparison of model results with NDVI
All of the models have been used in different applications, and their predictions have been compared to observations (see Sect. 2). However, this study is the first where the models have been used in PRA, except for a preliminary application using forest model BASFOR (Van Oijen et al., 2013b). Testing a probabilistic risk analysis is difficult: to test prob-abilities, you need many observations from which to derive frequencies. NDVI is useful for this due to the available long time series of observations: monthly 1982-2010 at high spatial resolution.
Although NDVI measures a variable that is different from the various model outputs, it is correlated with LAI and vegetation photosynthetic activity (Myneni et al., 1995), so we may expect that our results, if at all realistic, show some degree of correlation with NDVI. We show the correlations, over all grid cells, between the multi-annual averages for the seasonal (April-September) means for NDVI and each of the five model outputs (Table 3). Because these calculations are at the European scale, we should in fact expect reasonably strong positive correlations, reflecting the contrast between productive vegetation in central Europe versus the low LAI and vegetation activity in the dry Mediterranean area and the cold far north. This is borne out by the data in Table 3 which show generally high correlations for the four fluxes, with as sole exception the negative correlation of NDVI with ET predicted by PASIM. The relatively low correlations for the grassland and crop ecosystem models PASIM and EPIC are explained by the role of management which strongly determines grassland and crop biomass, e.g. spatially variable nutrient inputs which lead to high grassland and wheat yields in north Western Europe (Smit et al., 2008;Balkovič et al., 2013).
Changes in NDVI are expected to correlate well with drought impacts on leaf area and photosynthetic activity, so NDVI data should pick up the signal from severe droughts Reichstein et al., 2007;Zaitchik et al., 2006;Bevan et al., 2014). We therefore also applied the PRA to NDVI and tested for correlations between vulnerability and risks for NDVI with those for the fluxes of carbon and water through the plants (NPP and ET) ( Table 4). Note that, across Europe, the risk to NPP is not proportional to the level of NPP itself. Some regions with low NPP have high risk (Mediterranean), whereas other low-NPP regions have low or even negative risk (Scandinavia). This is a complex spatial pattern to reproduce, but we do see in Table 4 that vulnerability and risk for NDVI are correlated positively with those for both NPP and ET, again with the partial exception of PASIM. The latter exception may be related to the fact that grassland LAI can be partly uncoupled from weather conditions by grazing and cutting (compare Smit et al., 2008). These positive correlations thus give some degree of confidence in the validity of the risk calculations by the different models.

Sensitivity analysis of the PRA
We carried out three sensitivity analyses (Fig. 6). These were restricted to examining how the PRA for NPP in 1971-2000 and 2071-2100 as simulated by LPJmL would have changed with certain different choices in the relationship between env and sys.  Fig. 3a. Left column: sensitivity of the PRA to the SPEI thresholds used to define drought. Middle column: sensitivity of the PRA to shifting the start of the six month period over which the SPEI is calculated. Right column: sensitivity of the PRA to reducing climatic variability. Whole bars represent vulnerability, the black shaded parts are the risks. Both are expressed in relative terms, as fractions of NPP under non-hazardous conditions. The first sensitivity analysis addressed the threshold level for SPEI below which conditions are considered hazardous. Vulnerability decreased with increasing threshold, but it was less sensitive to these changes than the probability of exceeding the threshold, P (H ), so risks tended to increase with threshold level. In fact, if sys increases or decreases monotonously with env, then the vulnerability as defined here cannot change much when we select a different threshold. Vulnerability is even perfectly constant if the sys-env relationship is linear and the distribution of env (and thus sys too) is uniform over its range. In that case, the value of the vulnerability equals half the range over which the sys variable varies: V = 0.5 × (max(sys)−min(sys)). The finding that the vulnerability of NPP to drought was much more independent of the choice of the SPEI threshold in the years 2071-2100 compared to 1971-2000 thus indicates that with climate change the relationship of NPP to SPEI becomes more linear.
In this behaviour, where V becomes a constant in the linear case, it resembles a simple sensitivity or elasticity coefficient. However, it is a more useful concept than that in that it can handle non-linear relationships of any form. It is even possible to define hazardous conditions as the union of both very small and very large env values if we believe that the relationship of sys to env has the shape of an optimum curve.
In the second sensitivity analysis we varied the lag time between sys and env, as shown in the middle column of Fig. 6. Overall, vulnerabilities and risks decrease if we shift the drought period to earlier months and thus analyse the sys state with a lag time of up to three months. This suggests that our choice of fully contemporaneous sys and env timings, i.e. April to September for both, was the best choice, revealing the most pronounced impact of droughts on the vegetation.
In the third and final sensitivity analysis we examined how using a different climate scenario, with reduced weather variability in the future, would have affected the PRA. The differences were found to be very small, with only minor reductions in risk under the reduced variability scenario. However, we would need to repeat our PRA with various other climate scenarios to establish robustness against climate prediction error more confidently. -Climate change is expected to lead to large drought risks to primary productivity in the Mediterranean area, according to the six different vegetation models studied, which confirm earlier analyses (e.g. Schröter et al., 2005).
-The risks will increase mainly because of greater drought probability; ecosystem vulnerability will increase to a lesser extent.
-Future carbon sequestration (NEP) will also be at risk predominantly in the south because NPP will be more affected than Rh. Risks to NEP will exceed 0.25 g C m −2 d −1 according to most models, amounting to reductions in carbon sequestration of 20 to 80 %.
The current study can be expanded in various ways. First, in order to study the impact of extreme climatic conditions, other env variables than the SPEI can be considered, e.g. variables or a constellation of variables (Seneviratne et al., 2012) that quantify the strength of heat waves, frost periods or storms. On the side of the vegetation, we focused here on the carbon fluxes, but mortality, in particular of trees, may be an important vegetation property at risk that can be studied with our methods. Analysis of the responses of the ecosystem and generic models during past extreme weather events (e.g. 2003 Western Europe, 2010 Eastern Europe), within the risk analysis framework used here, could enhance the assessments of future risks and vulnerabilities. It may further be useful to attempt to investigate how vulnerabilities and risks would be reduced if different ecosystem management methods were adopted, or what the scope is for reducing risks by plant breeding.
The study we carried out here involved a total of close to 10 8 simulated vegetation years. However, this large number did not include a comprehensive uncertainty analysis. We did compare results between models (which can be seen as a form of uncertainty analysis with respect to vegetation model structure and indeed revealed large uncertainty about NEP in particular) and we examined climate uncertainty to some extent in the sensitivity analysis. But these efforts do not include all sources of uncertainty. Our risk analysis method is probabilistic by design, so it can easily accommodate uncertainties about vegetation response in the conditional distribution, P (sys|env), and climate uncertainty in the P (env). Further sources of uncertainty to be assessed in future work include parameter uncertainty for the individual models and uncertainty about model inputs relating to soils, management and land cover.