A novel probabilistic risk analysis to determine the vulnerability of ecosystems to extreme climatic events

We present a simple method of probabilistic risk analysis for ecosystems. The only requirements are time series—modelled or measured—of environment and ecosystem variables. Risk is defined as the product of hazard probability and ecosystem vulnerability. Vulnerability is the expected difference in ecosystem performance between years with and without hazardous conditions. We show an application to drought risk for net primary productivity of coniferous forests across Europe, for both recent and future climatic conditions.


Introduction
Climate change is expected to affect both average weather patterns and the occurrence of extreme events (IPCC 2012). Process-based simulation models are key tools for quantifying the risks posed to ecosystems by such events. However, there is no agreed method for analysing the results produced with such models in order to quantify vulnerabilities and risks. We propose a simple method here.
Terminology for risk analysis was developed under the auspices of the United Nations (DHA 1992). Risk was then defined as 'expected loss', to be calculated as the 'product of hazard and vulnerability'. Hazard was defined as 'the probability of occurrence of a potentially damaging phenomenon within a given time period and Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. area', and vulnerability as the degree of loss that results if the phenomenon does occur (DHA 1992). To facilitate quantitative analysis, we here make the further distinction between hazard as the potentially damaging phenomenon itself, and the probability of the hazard actually occurring. So risk is zero if the probability of hazard or the vulnerability is zero, and risk is only large when both components are large. A similar definition, 'risk = probability × consequence', was more recently used by the Intergovernmental Panel on Climate Change (IPCC 2012). These definitions were set up primarily to facilitate risk analysis for hazards threatening human life, but they can be made operational for ecosystem modelling too.
The risks we quantify here are those related to the carbon cycle of ecosystems, e.g. the expected loss of net primary productivity (NPP). The hazards we are interested in relate to meteorological variables, so we quantify the probability of hazardous weather conditions such as droughts and heat waves. Vulnerability, in this context, is the impact that hazardous weather has on the carbon fluxes. In the following, we shall refer to the quantification of risk and its two components, the probability of hazardous conditions and the vulnerability, as 'risk analysis'.
Our method of risk analysis can be applied to time series of empirical data and to ecosystem modelling results, but we focus on the latter here. In case the time series come from modelling, the hazards, vulnerabilities and risks to be quantified are restricted to the environmental input variables and ecosystem output variables of the model in question.
Environmental and ecosystem variables are characterized by uncertainty. Environmental variables may be measured or modelled, e.g. using GCMs, so they reflect uncertainties about measurement or modelling error. Simulated ecosystem variables reflect uncertainties in the inputs, parameters and structure of the ecosystem model. To account for these uncertainties, we base our risk analysis method on probability distributions for all variables. Our approach thus is as an example of probabilistic risk analysis (PRA; Bedford and Cooke 2001;IPCC 2012). The proposed methodology is inspired by ideas from the literature on PRA (going back to Kaplan and Garrick 1981), but we have tried to make it more general and more rigorously probabilistic. PRA tends to focus on risks associated with discrete hazardous events (tsunamis, storms, epidemics, etc) whereas we analyse the responses of ecosystems to any value that environmental variables such as rainfall or temperature can take. Our basic tools will therefore be continuous probability distributions rather than discrete ones. And although PRA commonly expresses the occurrence of hazardous events probabilistically, based on observed frequencies, it generally ignores the fact that the response to any such event is not precisely known. In contrast, we also express vulnerability to environmental stress probabilistically.
In short, this letter aims to make PRA operational for ecosystems analysis where both the environment and the ecosystem are imperfectly known and variable, so risks, hazards and vulnerabilities must be defined in probabilistic terms. We give the definitions in the form of mathematical equations and show that they allow writing risk as the product of hazard probability and ecosystem vulnerability. We conclude with an example illustrating how the method can be used to analyse how environmental change may alter risks posed by environmental change to carbon fluxes in ecosystems. The specific goal of the analysis is to quantify how drought risks to carbon sinks in coniferous forest may change across Europe, and whether the underlying causes will mainly be changes in drought frequency or changes in forest vulnerability.
In section 2.1, we present the underlying theory in a formal way; concrete suggestions for practical implementation are provided in section 2.2, and examples of application follow in sections 3 and 4.

2.
A method for probabilistic risk analysis using ecosystem models

Theory
We refer to the environmental variables which may or may not attain hazardous conditions as 'env', and to the ecosystem variables at risk as 'sys'. Our risk analysis is based on the probability distributions for these variables, denoted as P(env) and P(sys). The degree to which the environment accounts for the state of the system is represented by the conditional probability distribution P(sys|env). The three distributions are linked through the law of total probability: P(sys) = P(env)P(sys|env) d env. (1) Average environmental conditions and ecosystem behaviour are given by the expectations E(env) and E(sys). Our aim is to use the three probability distributions to define 'hazards', 'vulnerabilities' and 'risks' such that risk equals the product of hazard probability and vulnerability. We define hazard as a factor that can cause damage. Because any environmental variable can, on occasion, be too low or too high for optimal ecosystem performance, each constitutes a hazard. We shall say that an environmental variable is hazardous if it is in the range of values whose negative impacts we want to study. How to set that range will be discussed in section 2.2.
Vulnerability is defined here as a property of the conditional response distribution P(sys|env). It is the difference in expected system performance between hazardous and non-hazardous environmental conditions: where the conditional expectation values are calculated as: Equation (2) implies that vulnerability is defined as the average difference in system performance between 'good' and 'bad' conditions. Finally, risk is defined as the product of the probability of hazardous conditions and the ecosystem's vulnerability to such conditions: If we expand the vulnerability term in the last equation, using the definition given in equation (2), we can derive an alternative but mathematically equivalent formula: This alternative formula emphasizes that risk can be seen as the average loss that a system experiences due to occasional hazardous conditions. In conclusion, our definition of terms both obeys the common meaning of risk as 'expected loss due to a hazard' (equation (5)) and it allows decomposition of risk in the two terms given in equation (4).

Implementation of the method in ecosystem modelling studies
In practice, when we perform computer modelling of ecosystems, we do not calculate the theoretical distributions and integrals described in section 2.1, but we approximate them using results from model simulations. We propose the following six-step procedure for implementation of the risk analysis: (1) Select a relevant ecosystem model.
(2) Choose from the model's input variables (env) those that constitute hazards. (3) Set the threshold criteria for hazardousness.
(4) Choose from the model's output variables (sys) those that are considered at risk. (5) Set the boundary conditions for time and space and perform the model simulations. (6) Use the resulting input-output relationships and the equations of section 2.1 to calculate E(sys|env nonhazardous), E (sys|env hazardous), P (env hazardous), vulnerability and risk.
In step 1, we may decide to use more than one ecosystem model, which would make the analysis more comprehensive by taking into account uncertainty about model structure.
In step 2, the easiest choice would be to focus on individual input variables such as precipitation or temperature. Alternatively, we may decide to use composite environmental variables, such as weather indices which combine the values of multiple input variables, e.g. the Standardized Precipitation-Evapotranspiration Index (SPEI; Vicente-Serrano et al 2010). The analysis can also be performed using combinations of environmental variables without formally combining them in one index. The env-variable would then be multidimensional and P(env) would be a multivariate probability distribution. For example, the combined hazard probability of low rain and high temperature can be defined as P(env hazardous) = P(rain < c1 AND temperature > c2), where c1 and c2 are constants chosen by the analyst.
In step 3, we may select criteria for hazardousness from the literature, e.g. by choosing a common definition of meteorological drought. We can also base our criterion on the frequency distribution of the env-variable. An example would be to consider as hazardous any annual rainfall less than the current 25% quantile. In that approach, P(env hazardous) would be 0.25 under current climatic conditions, but it would change if the climate becomes drier or wetter. Note that if we would let the threshold vary over time, by always defining hazardousness with respect to the contemporary weather frequency distribution, P(env hazardous) would remain constant, and the analysis of risk in terms of its two components would become meaningless. We may however vary the reference climatic conditions across space. Amounts of rain that constitute a drought in wet regions may be above average elsewhere, but we can define as hazardous those weather conditions that are locally extreme.
In step 5, we define the region, spatial resolution and time period for which we carry out model simulations underpinning the risk analysis. We may be interested in the spatial distribution of hazards, vulnerabilities and risks, how they change over time, and how they differ between ecosystem types. In such cases, the ecosystem-specific environment and system distributions P(env) and P(sys) can be determined separately for individual spatial grid cells and by inspection of model outputs over limited periods of time, e.g. thirty years. Compared to P(env) and P(sys), the response distributions P(sys|env) will differ less between grid cells-because they are primarily determined by underlying universal mechanisms as represented in the models. If so, and if the selected env-variables cover the key sensitivities of the system, P(sys|env) can be derived from modelling results across larger regions or strata consisting of multiple similar grid cells.
In step 6, we use the modelling results to calculate risk and its components. This is done by interpreting the frequencies of input and output values as approximations of their probabilities. We suggest carrying out the calculations in the following order:

Example 1: analysis of virtual modelling results
The following is a simple example using virtual data. Say we have carried out simulations for a certain site for ten different time periods, and we have the results shown in table 1, sorted in order of increasing rainfall levels.
In this example, we have just one env-variable (rainfall) and one sys-variable (NPP). We assume here that P(env) and P(sys) can be approximated by the frequency distributions of the values in the table. For a serious analysis, we would need more than ten simulation results. We see that E(env) = 3 mm d −1 and E(sys) = 7.2 g m −2 d −1 . Say we define the hazardous range of env as values less than the average of 3 mm d −1 , so that the first three rows in the table represent hazardous conditions and the remaining seven non-hazardous ones. The expectation values for NPP then are 9 and 3 g m −2 d −1 for non-hazardous and hazardous conditions, respectively. Vulnerability is then calculated, using equation (2), as 9 − 3 = 6 g m −2 d −1 . The probability of hazardous conditions, P(env hazardous), is three out of ten = 0.3. Risk is the probability of hazardous conditions times the vulnerability (equation (4)), i.e. 0.3 × 6 = 1.8 g m −2 d −1 . We could also have calculated the risk using equation (5), as the expected NPP under non-hazardous conditions minus average NPP = 9 − 7.2 = 1.8 g m −2 d −1 . So the risk analysis has the following results: The vulnerability of this system is high, with drought reducing NPP by two-thirds on average. However, the risk is less because droughts only occur 30% of the time.

Example 2: analysis of drought risks across Europe using a forest model
We now give a more complex example. We use the dynamic model BASFOR to analyse risks posed by drought to coniferous forests across Europe. Because we focus on presenting the method for probabilistic risk assessment, rather than on the properties of any single model, we keep the following model description brief. BASFOR is a forest model that requires as input meteorological variables (radiation, temperature, precipitation, wind speed, atmospheric humidity), atmospheric CO 2 concentration and N-deposition rate. The model runs on a daily time step and simulates the dynamics of pools of carbon, nitrogen and water in tree organs and soil. For details about the model, including comparison with measurements in various applications, see (Van Oijen et al 2005, 2011. Here we consider the model's predictions of NPP (g C m −2 d −1 ). In a recent comparison of models for growth of Scots pine (Pinus sylvestris L) across Europe, the NPP-predictions from BASFOR were found to be the closest to independent observations from forest inventories and ecological research sites (Cameron et al 2013).
The model is applied on a latitude-longitude grid of 0.25 × 0.25 • across Europe, using two time periods: 1971-2000 and 2071-2100. For the first period, CO 2 concentrations and daily weather data were based on observations as assembled by the WATCH project and downscaled from 0.5 • using CRU CL 2.0. For the second period, the forcing variables were generated by the MPI-Remo model assuming the SRES A1B scenario for greenhouse gas emissions. Data on N-deposition across Europe for the year 2000 were downloaded from EMEP.
In this example, we define hazardous conditions as annual average precipitation being less than the 25% quantile (Q25). The absolute threshold for 'hazardousness' will thus be lowest in the driest regions. The results of the risk analysis over Europe are displayed in figure 1. Figures 1(A) and (B) show that NPP under both 'good' and 'bad' conditions (annual average precipitation larger or smaller than Q25) in the period 1971-2000 is highest in Central Europe. Vulnerability, the difference in NPP between good and bad years, is highest in eastern Spain and around the Black Sea ( figure 1(C)). The probability of bad conditions is everywhere 0.25 because of the Q25 criterion. Therefore the risk map (figure 1(E)) is just a rescaling of the vulnerability map, with risks being proportional to vulnerability.
We repeated this risk analysis for the years 2071-2100 (figure 2), but kept the threshold for hazardousness as the Q25 for 1971-2000. Across much of Europe, NPP is expected to increase under both non-hazardous and hazardous conditions (figures 2(A) and (B)), mainly because of elevated [CO 2 ]. However, the increases are not of the same magnitude in all areas and drought vulnerability is expected to increase in Southern Europe but decrease further north ( figure 2(C)). Drought probability changes in a similar way, with increases predominantly in the Mediterranean area and decreases in Northern Europe ( figure 2(D)). Overall, future drought risk will be highest in Spain and Anatolia ( figure 2(E)). These regions are already at relatively high drought risk in 1971-2000 ( figure 1(E)), but this is increased further due to the regional increases in both hazard probability and vulnerability.

Terminology of risk analysis
Many theoretical expositions and applications of risk analysis can be found in the literature (e.g. Adams 1995;Bedford and Cooke 2001). However, terminology may be confusing. An example is the concept which we refer to as 'vulnerability', defined by means of equation (2) as the difference between expected system performance under non-hazardous and hazardous conditions. The same concept is referred to by other authors as 'susceptibility', 'sensitivity', 'consequence' or 'loss'. Also, in the IPCC-context, the term vulnerability is defined in a different way, explicitly including people's ability to cope with a given situation (IPCC 2012). For ecosystems, the expression 'biophysical vulnerability' is sometimes used (Brooks 2003). Overviews of the many definitions in use for 'risk' are given by Brooks (2003) and Schneiderbauer and Ehrlich (2004). The definitions are generally not precise enough to allow mathematical analysis, e.g. where risk is 'a function of probability and magnitude of different impacts' (IPCC 2001). We have tried to circumvent terminological problems by the precise definitions and equations given in section 2.1.

Evaluation of the proposed method for risk analysis
The examples given in sections 3 and 4 show that risk analysis can be carried out using the inputs and outputs from a dynamic ecosystem model. Such models provide vast quantities of numerical results, especially when they are applied across a large region such as Europe and over a long time period. The risk analysis may help in extracting important information from the analysis. The method is internally consistent, grounded in probability theory (equations (2)-(5)), and allows for straightforward quantification and decomposition of ecosystem risks. Admittedly, the examples given in this paper were limited in scope. The European forestry study in section 4 used only one scenario for climate change, one ecosystem model, one env-variable and one sys-variable. Nevertheless, even this simple setup allowed us to address questions such as 'will risks increase in the future for the given climate change predictions?', 'If yes, will that be because of changes in the probability of hazardous conditions or because of changes in ecosystem vulnerability?', 'Which part of Europe will be most affected?' With minor additional effort, additional envand sys-variables could be included to analyse the risks more comprehensively. However, for some questions we would need to extend the exercise to other models, e.g. if we want to know which ecosystem types are at greatest risk.
The probability distributions underpinning our analysis represent both temporal variability of env and sys and uncertainty because of incomplete knowledge about inputs, parameters and model structure. In every risk analysis, we therefore need to ask whether all sources of variability and uncertainty have been accounted for. The example of section 4 used only one climate change scenario, so uncertainty about future climate, including the frequency and intensity of extreme events, was not represented. The environment distribution P(env) thus only accounted for the inter-annual variability present in this single scenario for Europe. The response distribution P(sys|env) also expresses both variability and uncertainty. Responses to given values of an environmental variable will vary because of interactions with the other environmental variables and because of differences in the state of the ecosystem itself. P(sys|env) should also account for uncertainty about model structure and parameterization. These sources of uncertainty can to some extent be quantified by the use of multiple models and by sampling from each model's probability distribution for parameters. The example of section 4 did not address these uncertainties.
The quality of the analysis depends on how well the models simulate ecosystem response to environmental conditions. This can be evaluated in a standard way, by comparing simulated versus observed sys-variables for as many environmental conditions as possible and combining this with model intercomparison. However, because the probability distributions in our vulnerability assessment are partly reflections of uncertainties about data and models, rather than of some true stochasticity in nature, comparison of simulated distributions and observed frequency distributions may be misleading.
The response of ecosystems to environmental conditions at any moment in time depends on the system state at that moment, and not on its past history. However, over time, the ecosystem state will keep changing as a result of internal processes directly or indirectly affected by the environment. Phenotypic plasticity allows individual organisms to tolerate climate variability. Genetic variability and natural selection induce long-term population-scale adaptation to climatic variability, e.g. in forest stands. Plant community structure, soil biota and food webs can be modified by climate change. Management options, especially in agricultural systems, may accelerate adaptation to climatic variability. The quality of the risk analysis will depend on the extent to which the model is able to calculate how much these various adaptation processes decrease ecosystem vulnerability.

Scope for application of the risk analysis
The risk analyst can choose from many envand sys-variables, only constrained by the available data or model input-output relationships. The variables studied in the risk analysis may be the original model variables or derivatives. For example, if an ecosystem model incorporates no explicit dependence on the SPEI drought index, we can still calculate the SPEI from the input data and use the model to quantify the associated risk.
There is also wide choice concerning the timescale of the envand sys-variables. We can calculate the risks posed by summer heat waves to ecosystem productivity in the same year, the following year, or the following decade, (see also Rolinski et al 2013). Such choices do not alter the applicability of the equations of section 2.1.
Our framework further allows wide choice for the definition of hazardous conditions. Conclusions from risk analysis may strongly depend on this choice (Arbez et al 2002), so it may be advisable to explore a range of criteria of hazardousness in any analysis. In the example of section 4, we used only one criterion for drought (the annual Q25 of rainfall), but adding other quantiles, timescales and variables would have provided a more robust analysis. It may also be informative to define hazardousness using disjoint ranges of environmental variables. For example, if we want to analyse the risk posed to NPP by both extremely low and high values of any env-variable, we could set P(env hazardous) = P(env < Q10 OR env > Q90). Compared to the one-sided risk analysis, this would result in both higher hazard probability and higher vulnerability in those regions where both extremes can occur.
The interaction between environment and ecosystems is complex. Risks associated with drought in any given year may be small if the ecosystems have had the opportunity to adapt to similar droughts in preceding years. Alternatively, risks may increase if the drought follows directly after another drought from which the ecosystem has not yet recovered. Such complexities do not limit the applicability of the risk analysis. We can carry out the analysis for stand-alone droughts and for droughts occurring in sequences. This only requires that such sequences actually do occur in the time series of our environmental data.

Conclusions and outlook
We have presented a simple probabilistic method for calculating hazards, vulnerabilities and risks using dynamic ecosystem models. The method only requires consideration of the model's time series of inputs and outputs. The method is easy to implement, yet allows formal decomposition of risk into hazard probability and ecosystem vulnerability. We aim to apply the approach presented here in project Carbo-Extreme, to investigate the impact of extreme weather events on carbon fluxes of diverse ecosystem types across Europe, for different scenarios of environmental change, and using a range of ecosystem-specific and generic models. To help understand the dynamics of carbon sources and sinks, the analysis will be carried out for net ecosystem exchange (NEE) and its two components: NPP and heterotrophic respiration.