A few extreme events dominate global interannual variability in gross primary production

Understanding the impacts of climate extremes on the carbon cycle is important for quantifying the carbon-cycle climate feedback and highly relevant to climate change assessments. Climate extremes and fires can have severe regional effects, but a spatially explicit global impact assessment is still lacking. Here, we directly quantify spatiotemporal contiguous extreme anomalies in four global data sets of gross primary production (GPP) over the last 30 years. We find that positive and negative GPP extremes occurring on 7% of the spatiotemporal domain explain 78% of the global interannual variation in GPP and a significant fraction of variation in the net carbon flux. The largest thousand negative GPP extremes during 1982–2011 (4.3% of the data) account for a decrease in photosynthetic carbon uptake of about 3.5 Pg C yr−1, with most events being attributable to water scarcity. The results imply that it is essential to understand the nature and causes of extremes to understand current and future GPP variability.


Introduction
Climate extremes such as droughts, heat waves or intense precipitation events are increasingly perceived as key players in the Earth system and are expected to increase in the wake of climate change (IPCC 2012, Barriopedro et al 2011, Sillmann et al 2013. The associated impacts on terrestrial ecosystems 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. are diverse and difficult to determine (Smith 2011, Ma et al 2012, Reichstein et al 2013. Some of these aspects are directly observable, for instance via crop damage (Piao et al 2010), and thus related to food scarcity (Rosenzweig et al 2001), or forest destruction by windthrow or fire (Chen et al 2011). Several studies have shown, however, that climate extremes can influence the terrestrial biosphere in less easily detectable ways, for instance by altering the carbon budget via a reduction of primary productivity (Piao et al 2010, Barriopedro et al 2011, Lewis et al 2011, Ma et al 2012. This aspect is crucial for an assessment and attribution of extreme anomalies in the carbon cycle. Hence, in this study, we aim to provide a global quantification of extreme anomalies in the photosynthetic carbon uptake, understand its implications for the carbon balance and attribute large-scale negative extremes to climatic drivers or fires. In contrast to earlier studies which analyzed climate extremes and their effects on the carbon cycle (Potter et al 2003, Reichstein et al 2007, Zhao and Running 2010, Barriopedro et al 2011, Lewis et al 2011, Ma et al 2012, our work is based on two major innovations. First, we adopt an 'impact perspective' (Smith 2011, Reichstein et al 2013 by directly evaluating extreme changes within the biosphere. In particular, we analyze four different data products of gross primary production (GPP) covering the last 30 years, two data-driven GPP estimates and two global carbon-cycle models. On the data-driven side we use the upscaled eddy-covariance-based GPP estimates (MTE, Jung et al 2011)  Our analysis then focuses on the largest extreme events where the impact of an event is determined by the GPP anomalies integrated over the spatiotemporal extent of the event. To estimate the impact of GPP extremes on the terrestrial net carbon balance we rely on net ecosystem exchange (NEE) estimates from Carbon-Tracker (Peters et al 2007) and estimates for the residual land sink (RLS, defined as the atmospheric increase in CO 2 minus emissions from fossil fuels minus net emissions from changes in land use minus oceanic uptake (Le Quéré et al 2013)).

Materials and methods
The core ingredients of this paper are two recently developed tools to first identify large-scale extreme events in variables characterizing the state of the biosphere and second attribute them to climatic drivers (Lloyd-Hughes 2012, Zscheischler et al 2013. The new element in the identification step is to combine individual extremes to three-dimensional extreme events to provide a robust assessment of large-scale extremes even when data sets contain uncertainties. In the attribution step we essentially test whether environmental conditions during a particular extreme event in GPP were extremely different compared to all other times in the available time period. After introducing the different data sets, in this section we explain in detail how we use these two approaches on our data sets.

Data
To identify extreme events that are relevant to the terrestrial biosphere, we rely on four different data sets describing gross primary production (GPP; all used datasets with their  abbreviations are listed in table A.1). MTE (Jung et al 2011) involves training a model tree ensemble at site level using FLUXNET (a global network of eddy-covariance observations in tandem with site level meteorology, Baldocchi et al 2001) to extrapolate to large spatiotemporal domains. We use a fully data-driven upscaling product that relies mainly on a composite of different remote sensing fAPAR products but also uses climate data from ERA interim. (fAPAR is the fraction of absorbed photosynthetically active radiation, a satellite remote sensing proxy for photosynthetic activity.) MOD17+ (Running et al 2004) is derived using the same model structure as the MODIS GPP data stream (Running et al 2000) linking shortwave incoming radiation, minimum temperature and vapor pressure deficit. The parameterization of the MOD17+ model follows Beer et al (2010); i.e., it is based on Bayesian inversion against GPP time series from FLUXNET. The model parameters are calibrated against GPP time series from the FLUXNET measurement network through a Bayesian data model synthesis. The terms in the MODIS-MOD17 biome-specific look-up table are used as priors. For regionalizing the model parameters we stratify the results of the in situ calibration per vegetation type and bioclimatic class. As climatic drivers we use the ERA-interim dataset and the same composite of fAPAR products as in MTE (Jung et al 2011).
LPJmL (Sitch et al 2003) is a dynamic global vegetation model (DGVM) with a fully coupled carbon and water cycle. Vegetation productivity, i.e. GPP, is derived by a process-based photosynthesis scheme that adjusts carboxylation capacity and leaf nitrogen seasonally and within the canopy profile. For the present study, LPJmL is run in its natural vegetation mode driven by ERA-interim temperature, radiation and precipitation.
OCN (Zaehle and Friend 2010) is a land-surface model derived from the ORCHIDEE DGVM (Krinner et al 2005), which prognostically simulates foliar area and N content and employs a two-stream radiation scheme coupled to the process-based calculation of photosynthesis in light-limited and light-saturated chloroplasts within each canopy layer.
To attribute negative extreme events in GPP to (extreme) climatic drivers and fire events we use temperature (T ) and precipitation (P) from ERA-Interim (Dee et al 2011), burned area (BA) and CO 2 emissions from fires (FE, Giglio et al 2010), and the water availability index (WAI). WAI is a surrogate for soil moisture and was computed using daily precipitation and potential evapotranspiration data from bias-corrected ERA interim and a map of plant available water holding capacity from the Global Harmonized World Soil Database (FAO/IIASA/ISRIC/ISSCAS/JRC 2012). Using WAI or variants of it as a proxy for soil moisture has a long tradition in ecosystem modeling (Federer 1982, Aber and Federer 1992, Prentice et al 1993, Kleidon and Heimann 1998. The parameter used in the calculation of the WAI was taken as the median value of 15 sites from the synthesis study of in situ data by Teuling et al (2006). We use low values of water availability as an index for drought.
The spatial resolution for all the above data sets is 0.5 • . T , P, WAI and all GPP data sets are available monthly from 1982 to 2011, BA and FE from 1997 to 2010. Table A.1 gives a summary of all used data sets of this gridded format.
To relate extremes in GPP to net ecosystem exchange (NEE) we downloaded monthly data from CarbonTracker for the years 2000-2011 (CT2011 oi, http://carbontracker.noaa.g ov, Peters et al 2007) and aggregated them to global anomalies. To obtain the residual land sink (RLS) we rely on the Global Carbon Project (GCP, www.globalcarbonproject.org) and use yearly data from 1982-2011 (Le Quéré et al 2013).

Preprocessing
We follow the suggestions from Zscheischler et al (2013). In particular, for T , WAI and all GPP data sets we first subtract the linear trend and mean annual cycle per grid cell. We divide the less smooth variables P, BA and FE at each pixel by the sum of the respective time series at that grid cell. We call the resulting values after preprocessing anomalies. The goal of the preprocessing is to obtain time series which depict deviations from the mean behavior.

Extreme event identification
In accordance with the IPCC climate extreme classifications , we define extremes as the occurrence of certain values in the tails of the probability distribution of the anomalies. We adapt the methodology described in Zscheischler et al (2013) and define extremes to be outside a certain threshold q, which is defined by a percentile on the absolute values of the anomalies (figure A.2). We choose the thresholds for each of the four data sets such that extremes (positive and negative together) comprise 1%, . . . , 10% of the anomalies, respectively. We then define an extreme event by spatiotemporally contiguous points whose values are larger than q (positive extremes) and smaller than −q (negative extremes), respectively. The size of an extreme event is defined by the integral of the corresponding anomalies of GPP over time and space, i.e. its unit is grams of carbon. Hence, 'large' events are determined by their integrated impact of GPP anomalies over their spatiotemporal domain.
To compute correlations with monthly time series of other data sets or map hot spots of extreme events we aggregate extreme events of each GPP dataset individually either in space or in time. For an aggregation in space, at each time step all anomalies which happen to be part of some extreme event (according to the definition under consideration) are summed together (see, e.g., figure 1), yielding one time series of the overall impact of the respective extreme events. This time series of GPP extreme events can then be correlated with e.g. the global GPP or NEE anomaly in order to obtain the fraction of explained variance. Similarly, for an aggregation in time, at a specific location (pixel), all anomalies which happen to be part of some extreme event are summed together (see, e.g., figure 3), yielding a global map of the cumulative impact of the extreme events under consideration. We further define the asymmetry between negative and positive GPP extremes as the quotient of the size of n negative events over n positive events. Throughout the text, we usually report the fraction exceeding unity in percent.

Attributing drivers to GPP extremes
We intend to identify drivers that possibly caused negative extreme events in GPP. To this end, for a certain negative GPP extreme event we compute the median of a driver variable over the spatiotemporal domain of the event. By shifting the event in time and computing such medians for each possible time step, we obtain a test statistic for each driver variable and each GPP extreme (see also Zscheischler et al 2013). We then compute p-values for the 100 largest GPP extreme events using the first-percentile definition by considering the drivers T , P, WAI, BA, and FE. For T and P we take both right-and left-sided p-values into account (heat waves and cold spells, excessive and exceptionally low precipitation). For WAI we consider only the left-sided p-value (water scarcity), for BA and FE only the right-sided p-value (fire events). We count a driver as a cause if its p-value is smaller than 0.1 (see table A.5 for a summary of significant drivers). Other thresholds are possible but do not affect the conclusions of this study. If a driver's p-value is below 0.1 one or several months before the GPP extreme event, we assume that the GPP extreme is a lagged response to the extreme driver. In this study we take time lags of a maximum of three months into account.

A few extremes explain most of the global interannual variability
The size distribution of extreme events can often be well approximated by a power law relationship f (x) ∼ x −α , where x is the size of an extreme event and α is the so-called scaling exponent (Ghil et al 2011). It has been recently shown that extremes in fAPAR can be likewise approximated well by a power law (Zscheischler et al 2013, Reichstein et al 2013). We also find power laws for the size distribution of extremes in GPP (here, the size of an extreme is its impact on GPP, figure A.3). The scaling exponent α is below 2 for most data sets and percentiles, indicating that a few exceptionally large extremes dominate the whole distribution of extremes (Fisher et al 2008).
Assuming that local extreme events in GPP can be of global relevance, the discovered power law behavior suggests that only a few events cause globally relevant impacts. To investigate this hypothesis, we correlate time series of integrated GPP extremes with the global GPP anomaly. We find that 78% (±5) of the global GPP anomaly can be explained by the largest 200 tenth-percentile GPP extreme events (positive and negative, figure 1(A) for the case of MTE; figure A.4(A) for all data sets). These extreme events occur on only 7% (±0.6) of the spatiotemporal domain (figure A.4(B)), revealing the strong spatial heterogeneity inherited from the power law distribution. By choosing a subset of the largest events we at the same time focus on the more robust features of the data sets.
The question is, however, how the identified extremes in GPP translate into anomalies of the land carbon balance. We find that same extreme events (the 200 largest tenth-percentile extreme events in GPP) explain 8.4% (±0.03, p < 0.001) of the variability in monthly anomalies of global NEE over the years 2000-2011. In addition, in the case of MTE, the same extreme events explain 22% of the annual variability in RLS for the entire time period ( p < 0.01, figure 1(B)). MTE is the data set closest to actual observations (it uses site level CO 2 flux observations compiled in FLUXNET in tandem with satellite derived data of fAPAR), which might explain its prominent role in yielding significant correlations with RLS. Very strong correlations between extremes in GPP and variability in NEE, or even the atmospheric CO 2 growth

Negative extremes are larger than positive extremes
The decrease of GPP associated with negative extremes is globally only partly compensated by positive GPP extremes. Negative GPP extreme events are on average up to 25% larger compared to positive extreme events if we consider the 1000 largest events on both tails of the distribution (figure 2(A)). While the data spread is large (ranging from 0% in LPJmL to 80% in MOD17+ for the first-percentile definition, figure 2(A), tables A.2 and A.3), all models show an asymmetry towards more decrease in gross land carbon uptake (negative extremes) for nearly all used percentiles. Moreover, this asymmetry is stronger if one considers fewer of the largest events and more extreme percentiles (figures 2(A) and A.5). According to the 'slow in, rapid out' principle for net ecosystem exchange (Körner 2003), the observed asymmetry might originate in e.g. a few large disturbance events such as droughts, fires or insect outbreaks that lead to often an instantaneous decrease in carbon uptake (negative GPP extremes). Excess productivity (positive GPP extremes, e.g. due to regrowth), in contrast takes place on much longer timescales (Bengtsson et al 2003, Dore et al 2008. The observed asymmetry could also be originated in asymmetric driving mechanisms (e.g. climate extremes). Skewed distributions in drivers due to amplifying feedbacks such as those between droughts and heat waves (Mueller and Seneviratne 2012) might cause disproportionally large decreases in carbon uptake.
We will focus on negative GPP extremes for most of the remainder. The cumulative effect of negative extremes depends strongly on the percentile used. Averaged over the four GPP data sets, the 1000 largest negative extremes based on the first-percentile definition yield a decrease in gross land carbon uptake of about 0.7 Pg C yr −1 , in contrast to 2.2 Pg C yr −1 on the fifth-and 3.5 Pg C yr −1 on the tenth-percentile definition (figure 2(B); for comparison, yearly photosynthetic carbon uptake ranges between 110 and 160 Pg in our data sets). Both the estimates for the two data-driven approaches as well as the estimates from the terrestrial biosphere models agree well among each other. However, the estimates of the models exceed the data-based estimates by a factor of more than two ( figure 2(B), table A.2). The difference in the size of the extremes between these two sets is mainly due to different magnitudes in the anomalies (table A.4) but the timing of extremes is highly correlated (figure A.6). It has been shown before that some models tend to overestimate interannual variability (Keenan et al 2012), while data-driven approaches rather underestimate interannual variability (Jung et al 2011). We hence assume that data-driven and process-driven estimates of GPP describe a reasonable envelope for the real size of the extremes.
The geographical hotspots of negative GPP extremes are found in Northeastern Brazil (matching with earlier observations, Potter et al 2003), Central South America, Southeastern Australia, South Africa, Kenya, Tanzania and South Central United States (up to more than 180 g C m −2 yr −1 decrease in GPP, figure 3(A)). In regions of Eurasia and North America the absolute decrease in GPP during extreme events is smaller (up to 100 g C m −2 yr −1 ) because the average GPP is lower. The individual patterns for the four data sets broadly agree with each other, although there are significant differences at a regional scale ( figure A.7). MOD17+ is the only data set showing extremes in the Amazon. MTE shows a hot spot in Eastern China, which is not seen by the other data sets ( figure A.7). Overall, the hot spots agree well with negative extreme events on remotely sensed variables describing the state of the biosphere such as fAPAR (figure 6 in Zscheischler et al 2013) and the enhanced vegetation index (EVI, figure A.8, obtained with the same methodology). To address the relative changes in gross land carbon uptake due to negative GPP extremes, we divide the globe into the 26 IPCC regions (IPCC 2012) (figure A.9). Northern Australia (region 25) on average experiences the largest relative decrease in gross land carbon uptake (nearly 7%) followed by the Central United States (region 8) with around 5.5% decrease (figure A.10). While we find most positive and negative GPP extremes in similar areas, particularly in Europe and in tropical areas negative GPP extremes are dominating ( figure 3(B)).

Most negative extremes are driven by water scarcity
Various environmental drivers can lead to extreme reductions in GPP (Reichstein et al 2013). We investigate extreme temperatures, extreme precipitation, droughts or fires as possible drivers using recently developed tools (Zscheischler et al 2013) and find that negative GPP extremes are most often associated with anomalous low values of water availability (between 58 (MTE) and 93 (MOD17+) events out of the 100 largest, figure 4). Associations with extreme temperature and precipitation are not found that often. GPP in MOD17+ is most susceptible to fire: 18 out of the 63 negative GPP extremes that fall in the range where fire data is available (29%) can be associated with fire (figure 4, table A.5). Differences in the number of events associated with a certain driver reflect varying climate sensitivity between the four data sets. Overall, we associate in each data set more than 70 of the 100 largest negative GPP extremes to extreme temperatures or precipitation, low water availability, or fires (table A.5, last column). In line with regional studies (Zhao and Running 2010, Ma et al 2012, we conclude that on a global scale, water deficit is the main driver for negative GPP extremes (figure 4). High temperature extremes and heavy precipitation events play a subordinate role at the global scale although they may be of regional relevance. The data products and models used, however, may also underestimate the effects of e.g. extreme high temperatures on vegetation activity and legacy effects. Other causes for extreme anomalies in GPP might be large disturbances such as pest outbreaks, extreme winds, and human deforestation.

Conclusions
Understanding the behavior and characteristics of extreme events in climate and carbon fluxes are important for the projections of the current land carbon sink (Reichstein et al 2013) and future vulnerability assessments (IPCC 2012). Directly studying extreme responses in the biosphere opens innovative possibilities for research within the disturbance and carboncycle science communities. Our study reveals that a few extreme events in GPP events explain most of its interannual variability and contribute a small but significant amount to the variability in the net carbon balance, although they occur on only a small fraction of the land surface. These results imply that to understand current and future GPP variability it is essential to understand the nature and causes of extremes. Our results also highlight the importance of hydrometeorological extremes for carbon-cycle variability at the global scale. Thus, expected alterations of the water cycle in future (Sillmann et al 2013, Huntington 2006) will likely have a large effect on carbon-cycle extremes and global interannual variability.     A symmetric threshold q is set such that (e.g.) 90% of the data anomalies fall in between −q and q. Those values which exceed the threshold are defined to be extreme. In this example the extremes are defined using the tenth-percentile definition (not to scale).   Shown is the quotient of accumulated carbon of the n largest negative events over the n largest positive events for different percentiles (1%, . . . , 10%, n = 10 (squares), 100 (diamonds), and 1000 (stars)) averaged over the four data sets MTE, MOD17+, LPJmL, and OCN.   . Shown are the 1000 largest negative tenth-percentile extremes in EVI. Since EVI is an index between 0 and 1, the units do not have a direct interpretation in gross land carbon uptake. Extremes were computed using the same approach as for GPP.