Climate-mediated spatiotemporal variability in terrestrial productivity across Europe

Quantifying the interannual variability (IAV) of the terrestrial ecosystem productivity and its sensitivity to climate is crucial for improving carbon budget predictions. In this context it is necessary to disentangle the influence of climate from impacts of other mechanisms underlying the spatiotemporal patterns of IAV of the ecosystem productivity. In this study we investigated the spatiotemporal patterns of IAV of historical observations of European crop yields in tandem with a set of climate variables. We further evaluated if relevant remote-sensing retrievals of NDVI (normalized difference vegetation index) and FAPAR (fraction of absorbed photosynthetically active radiation) depict a similar behaviour. Our results reveal distinct spatial patterns in the IAV of the analysed proxies linked to terrestrial productivity. In particular, we find higher IAV in water-limited regions of Europe (Mediterranean and temperate continental Europe) compared to other regions in both crop yield and remotesensing observations. Our results further indicate that variations in the water balance during the active growing season exert a more pronounced and direct effect than variations of temperature on explaining the spatial patterns in IAV of productivity-related variables in temperate Europe. Overall, we observe a temporally increasing trend in the IAV of terrestrial productivity and an increasing sensitivity of productivity to water availability in dry regions of Europe during the 1975–2009 period. In the same regions, a simultaneous increase in the IAV of water availability was detected. These findings suggest intricate responses of carbon fluxes to climate variability in Europe and that the IAV of terrestrial productivity has become potentially more sensitive to changes in water availability in the dry regions in Europe. The changing sensitivity of terrestrial productivity accompanied by the changing IAV of climate is expected to impact carbon stocks and the net carbon balance of European ecosystems.


Introduction
One of the largest sources of uncertainty in modelling the future global climate and carbon cycle changes is the response (hereafter called sensitivity) of the terrestrial ecosystem productivity to climate change and variability.Comprehensive understanding of the sensitivity of the interannual variability (IAV) in different productivity terms to climate will provide crucial insights into future features of the terrestrial carbon balance and its climate feedbacks (Cox et al., 2013;Govindasamy et al., 2005).During the past decades, great efforts were devoted to investigating the IAV in biosphereatmosphere net carbon exchange and the underlying mechanisms (Kaplan et al., 2012;Moors et al., 2010;Schimel et Published by Copernicus Publications on behalf of the European Geosciences Union. X. Wu et al.: Climate-mediated variability in terrestrial productivity across Europe al., 2001).From these studies, it is evident that the terrestrial carbon uptake responds to climate variations and trends at a global scale (Heimann and Reichstein, 2008).The factors controlling terrestrial productivity, its magnitude and spatiotemporal variability are still poorly quantified (e.g.Keenan et al., 2012) and understanding them better is a prerequisite for constraining the variability of net carbon fluxes.
It has been demonstrated that several processes including land use changes, natural and/or anthropogenic disturbances, climate change and extreme events, as well as CO 2 variations contribute to the variability of terrestrial productivity (Houghton, 2000;Kaplan et al., 2012;Schimel et al., 2001;Zscheischler et al., 2014).Among those processes, land use changes and CO 2 / N fertilization contribute primarily to long-term changes in terrestrial productivity, i.e. on decadal to centennial timescales (Houghton et al., 1999;Jungclaus et al., 2010;Kaplan et al., 2012;Schimel et al., 2001;Zaehle et al., 2010).In contrast, the effects of climate variations are probably the most important drivers of IAV in terrestrial productivity on shorter (e.g.seasonal to decadal) timescales (Barford et al., 2001;Houghton, 2000).This categorization has been supported by data from diverse ecosystems, including forests (Richardson et al., 2007;Tian et al., 1998;Yuan et al., 2009), grassland (Craine et al., 2012;Flanagan et al., 2002;Jongen et al., 2011;Suyker et al., 2003), cropland (Lei and Yang, 2010), and tundra (Schuur et al., 2009).However, few studies have systematically investigated the sensitivity of terrestrial carbon uptake to interannual variations in climate during past decades, especially on regional to continental scales (Schwalm et al., 2010).
Despite some progress, the sensitivity of terrestrial productivity to climate IAV continues to be a significant source of uncertainty hindering accurate carbon cycle predictions across biomes.For instance, northern high-latitudinal vegetation not only showed non-linear responses to climate variability in recent decades, but also showed biome-specific characteristics of feedback response (Fang et al., 2005;Jeong et al., 2013).These studies indicate that the climate sensitivity of the IAV of terrestrial productivity may even depend upon climate state -a characteristic with important implications in the context of projected climate change (IPCC, 2012).
Insights into the past IAV of ecosystem productivity can be obtained from comparing multiple historical proxies.In this study, we investigate the IAV of crop yield records covering the past four decades and their sensitivity to interannual climate variability in Europe.In order to obtain a more robust picture of the underlying dynamics, we accompany the study with an analogue analysis of remote-sensing observations.Specifically, we aim to (1) identify the spatiotemporal patterns in IAV of terrestrial productivity across Europe from crop yields and remotely sensed observations, and (2) assess the possibly changing sensitivity of terrestrial productivity to climate IAV over the past decades.

EUROSTAT yield data
We obtained the annual EUROSTAT regional crop yield statistics (according to administrative boundaries at the NUTS 2 level) for the period 1975-2009 from the European Commission (http://epp.eurostat.ec.europa.eu).This data set documents all primary crop types, such as barley, wheat, grain maize, etc.However, data for most of these crops are scarce.In this study we thus focus on four major crops with rather continuous spatiotemporal coverage: barley, wheat, grain maize and potatoes.In the EUROSTAT database, yield is reported as the amount of dry matter suitable for consumption (Moors et al., 2010).
We first examined the length of the available yield records for each crop and region, thereby discarding time series below 10 yr in length.We subsequently filled gaps in the crop yield time series for the remaining regions using a simple linear interpolation method.Only short data gaps (length of individual data gap ≤ 3 yr) for regions with relatively few discontinuities (sum of gaps < 15 % of the total time series) were filled via linear interpolation, whereas regions with longer gaps were discarded.Our evaluation showed that this gap filling exerts only minor effects on crop yield IAV (data not shown).Finally, we rasterized the yield data for each kind of crop from the regions at NUTS2 level to obtain 0.5 • × 0.5 • gridded yield data matching the resolution of the climate, land cover, and other productivity proxy data (see below).We calculated the gridded average crop yield for each crop type, ȳ, using the following equation: , where A k is the area fraction of the kth region in one grid and C k is the crop yield of region k obtained from EURO-STAT.In cases where extreme values were present in C k (i.e.above or below three standard deviations of the crop yield series for region k during 1975-2009), we assigned the corresponding A k as zero to minimize possible bias introduced by unrealistic yield records prior to the rasterizing process.We only retained grid points with a cropland fraction ≥ 5 % based on the International Geosphere Biosphere Programme (IGBP) vegetation classification scheme (Loveland et al., 2000).

NDVI, up-scaled GPP, and FAPAR
Bimonthly normalized difference vegetation index (NDVI) and fraction of absorbed photosynthetically active radiation (FAPAR) data derived from the Advanced Very High Resolution Radiometer (AVHRR) of the National Oceanic and Atmospheric Administration (NOAA) at a spatial resolution of 8 km during the period 1982-2008 were obtained from the NASA GIMMS (Global Inventory Modeling and Mapping Studies) group.The GIMMS NDVI and FAPAR data sets have been thoroughly corrected for orbit drift, clouds and atmospheric aerosols (Tucker et al., 2005), and have already proved useful to identify long-term changes in vegetation greenness/activities (Myneni et al., 1997;Nemani et al., 2003;Zhou et al., 2001).In this study, we resampled the GIMMS NDVI and FAPAR data into monthly data with a spatial resolution of 0.5 • × 0.5 • .The cumulated NDVI and FAPAR during the active growing season is regarded as a good indicator for vegetation activity (e.g.Barichivich et al., 2013).Note that for comparison with other data streams, the NDVI and FAPAR data were averaged across all vegetation types present in each grid cell.
Gridded monthly gross primary production (GPP) data over Europe at 0.5 • × 0.5 • spatial resolution during the 1982-2008 period were obtained from Jung et al. (2011).This data set was constructed by integrating in-situ measurements of FLUXNET data from the La Thuile 2007 synthesis (www.fluxdata.org)with satellite remote-sensing data and meteorological reanalysis (Beer et al., 2010).This observation-based estimation has been demonstrated to add confidence and spatial detail to the global terrestrial GPP (Beer et al., 2010).The gridded GPP product is not independent of FAPAR, given that FAPAR is a key variable in the upscaling of point-wise FLUXNET measurements.

Climate data
Monthly temperature (TMP) data from 1975 to 2009 at 0.5 • × 0.5 • spatial resolution were derived from the CRU (Climatic Research Unit) TS (time series) 3.1 gridded global data set (http://www.cru.uea.ac.uk/).The 0.5 • × 0.5 • gridded water availability index (WAI) data from 1975 to 2009 were estimated based on a simple formula modified from Kleidon and Heimann (1998).Monthly reanalysis data of the number of consecutive frost days (CFD) and the maximum of daily maximum temperature within a month (MMT) for more than 8000 stations over the globe were acquired from the European Climate Assessment (ECA) project (http: //www.ecad.eu).We then constructed monthly gridded data for CFD and MMT at a spatial resolution of 0.5 • × 0.5 • for the period 1975-2009 by performing a cubic spline spatial interpolation (Tabor and Williams, 2010).

Data preprocessing
Long-term trends in specific crop yield series during the period 1975-2009 in different regions of Europe are likely caused by technical progress (e.g.breeding), land use policy changes and changes in management practices.The latter could alter the photosynthetic capacity among different crops, even at the same location (Moors et al., 2010).Therefore, we fitted cubic smoothing splines with a frequency cut-off of 50 % at 32 yr to each series of crop yield and produced a set of dimensionless indices for different kinds of crops with a mean of one.This was achieved by dividing each crop yield time series by the fitted values, which is analogous to standard dendrochronological procedures, where climate unrelated low-frequency trends occur as an effect of tree geometry (Babst et al., 2012;Cook and Peters, 1981).
Accordingly, we used the same normalization method to subtract long-term trends in NDVI, GPP, FAPAR and climate data on each grid.The resulting indices mainly preserve high-frequency (i.e.interannual) variability, but not the trends or low-frequency signals (Babst et al., 2012;Cook and Peters, 1997).
The growing season we refer to here is the season of active growth, which is vital for vegetation activity and crop production.Although phenology differs across Europe, even for the same kind of crop (Sacks et al., 2010), we attempted to define the growing season in a consistent way.Our approach is supported by observational data showing that the differences in both the planting and the harvest dates for crops are generally within 1 month between central and southern Europe (Sacks et al., 2010).The active growing season is estimated in this study at monthly rather than daily resolution, which could further mitigate the differences in the active growing-season intervals among regions in Europe.Consequently, we defined the active growing season as March-June for barley and wheat (hereafter, spring crops) and as March-September for grain maize and potato (hereafter, summer crops).

Spatiotemporal patterns of the IAV for crop yield, NDVI, GPP and FAPAR
We evaluated the IAV of detrended crop yield, and growingseason NDVI, GPP and FAPAR (called NDVI gs , GPP gs , and FAPAR gs hereafter) for each grid by calculating the coefficient of variation (CV), which has been demonstrated as an effective measure of year-to-year variability (Cao et al., 2003;Galloway, 1985).Regions with extensive irrigation (with irrigated area fraction ≥ 15 %) were masked out in our analysis to minimize the effects of human management on the IAV of productivity proxies.All data sets are compared on an annual basis.We subsequently quantified spatial patterns of the IAV for different productivity proxies.Spearman's rank correlation coefficients were used to assess the relationships between the IAV of crop yield, NDVI gs , GPP gs and FAPAR gs for each grid.Correlations of the IAV in these different observations allowed us to investigate possible intrinsic common signals and hence provide insight into the general patterns of IAV of terrestrial productivity across proxies.
The temporal changes in the IAV of crop yield, NDVI gs , GPP gs and FAPAR gs in each grid were investigated by

Climate relations of IAV of productivity proxies
Response function and generalized linear model analyses were performed to assess the principal relationships between the IAV of different productivity proxies and the corresponding IAV of detrended climate variables including TMP, WAI, CFD and MMT.The response function analysis is a multiple regression technique using principal components of climate data to estimate the values of dependent variables (Fritts and Wu, 1986;Guiot, 1991).CFD only during spring (March-May) was considered in this analysis because frost damages primarily affect vegetation activity and crop yield during this season.

Changing climate sensitivity of IAV of productivity proxies
For each grid, we estimated the interannual sensitivity of each productivity proxy to each climate factor, γ IAV , based on the least-squares regression formula: where CI j is the detrended time series of each productivity proxy for period j , T j is the detrended climate variable (either TMP or WAI) for period j , and ε j is taken as random noise.Strong linkage between CI j and T j yields higher absolute values of γ IAV j .In this analysis, only the growing season's mean temperature and total growing-season WAI were taken into account, since these two climate variables are the most important factors determining the IAV of terrestrial productivity (see result section).Temporal trends in γ IAV for each grid were evaluated using moving least-squares regressions with a 10 yr sliding window (1 yr lag).The Theil-Sen slope method was again applied to estimate the magnitude and significance of the trends in the changes of γ IAV (Yue et al., 2002).
The distributions of differences in γ IAV of carbon proxies to TMP and WAI between all 10 yr sliding windows in 1982-1995 and those in 1995-2008 (difference matrix of 15 elements × 15 elements for γ IAV ) were analysed for different climate zones to further illustrate regional differences in γ IAV trends.

Spatiotemporal patterns of IAV in crop yield, NDVI, GPP and FAPAR
The IAV of different productivity proxies shows a consistent spatial pattern, except for potato yield (Fig. 1).Evidence from historical yield observations reveal a much larger (p < 0.05) IAV in southern (e.g.south of 45 • N) than in central Europe, especially for spring crops (∼ 68 and ∼ 80 % for barley and wheat, respectively) (Fig. 1).Similar patterns are also observed in the IAV of both March-June-and March-September-summed NDVI, GPP and FAPAR (Supplement Fig. S1).Frequency distribution analyses confirm that the IAV of crop yield is 36.5 and 70.8 % and 30.6 and 97.8 % larger (p < 0.05) in drier climate zones (i.e.climate zone Cs) than that in the Cf and Df climate zones for barley and wheat, respectively (Supplement Fig. S2).However, we did not find such patterns in the yield IAV for the two summer crops: maize and potato (Supplement Fig. S2).
Spearman's rank correlation analysis shows a good relationship between the IAV of different productivity proxies during the 1982-2008 period, especially for regions/crops that have a larger IAV in yield records (Figs. 2, 3).Notably, there is a tendency towards much stronger relationships between IAV of crop yield and IAV of NDVI gs , GPP gs and FAPAR gs in the drier regions of Europe (Supplement Fig. S3).In these dry regions, a much higher fraction (∼ 45.5-60.7 %) of grid cells shows significant (p < 0.1) and positive correlations (∼ 0.4-0.9) between the IAV of productivity proxies than in other regions (∼ 15-30 %) (Fig. 4).
Despite large differences in the spatial patterns of the linear trends in temporal changes of crop yield IAV, the IAV of most crops (except potato) displays a markedly increasing trend in water-limited regions (e.g. the Cs climate zone and some regions in eastern Europe) (Fig. 5).In these dry regions, the linear trends in CV of accumulative growingseason NDVI, GPP and FAPAR also show a consistent increase during both March-June and March-September seasons (Supplement Figs.S4, S5).

Relationships between spatiotemporal IAV of productivity proxies and climate
A response function analysis reveals that the IAV of terrestrial productivity in most parts of temperate and Mediterranean Europe responds strongly to the IAV of water availability and negatively to the IAV of temperature, as illustrated  by a significant positive (negative) productivity IAV response coefficient to WAI (TMP) IAV (Fig. 6, Supplement Figs.S6-S11).This conclusion is confirmed by the results of a generalized linear model, which reveals that the spatial patterns in IAV of productivity proxies generally show significant and positive relationships to the spatial patterns in IAV of WAI in water-limited regions (such as Cs climate zone) and even in some of the temperate humid climate regions (Table 1).The response of different productivity proxies to IAV of CFD did not show unequivocal spatial patterns (Fig. 6, Supplement Figs.S6-S11).In contrast, all productivity proxies show a consistently significant negative response to the IAV of MMT in most parts of Europe (Fig. 6, Supplement Figs.S6-S11).

Changes in the climate sensitivity of productivity proxies
Despite the great spatial heterogeneity, a generally negative crop yield sensitivity γ IAV to mean growing-season temperature IAV was obtained during the 1975-2009 period, particularly in eastern and southern Europe (Fig. 7, Supplement Fig. S12).Consistent with this finding, a positive γ IAV of IAV of crop yield to total growing-season WAI was obtained (Fig. 7, Supplement Fig. S12).Uniform patterns in the general climate sensitivity of productivity proxies such as NDVI, GPP and FAPAR were also observed (Supplement Fig. S13).Notably, there is a pronounced change in γ IAV of IAV of crop yield in Europe towards present (Fig. 8).Specifically, there is a generally increasing crop yield γ IAV for both TMP and WAI in water-limited regions in Europe (Fig. 8, Supplement Figs.S14, S15).Increasing climate sensitivity of IAV of cumulated March-September NDVI, GPP and FAPAR to mean March-September temperature and total March-September WAI is also observed in these regions (Supplement Figs.S16-S18).Additionally, there seems to be a generally increasing (but spatially variable) climate sensitivity of IAV of terrestrial productivity to changes in WAI in southern and eastern Europe (Figs.7, 8, Supplement Figs.S12-S18).Further analysis illustrated that the IAV of productivity in dry-summer temperate and Mediterranean climate zones (i.e. the Cs climate zone) is more sensitive to changing temperature and water availability in the more recent 1995-2008 period than during the 1982-1995 period (Fig. 9, Supplement Fig. S19).

Spatiotemporal patterns of IAV of terrestrial productivity and its linkage to climate variations
Evidence from most crop yield observations and remotesensing retrievals clearly indicates that the IAV of terrestrial productivity in Europe is spatially heterogeneous, with much higher IAV in water-limited regions (e.g.southern and eastern Europe) compared to humid oceanic regions (Figs.1-4, Supplement Figs.S1, S2).Variations in the water balance (approximated here by the WAI simple index) exert much more important and direct controls (signified by stronger rank correlations) on the spatial patterns in IAV of terrestrial productivity in temperate Europe (i.e.climate zones Cs and Cf) (Table 1) than temperature variations.Water availability is a crucial factor controlling terrestrial productivity at regional scales, especially in arid and semiarid regions (Austin et al., 2004;Huxman et al., 2004), whereas temperature affects productivity in such regions mainly by modifying the water availability, evaporative demand and vapour pressure deficit (Williams et al., 2013;Zhao and Running, 2010).Despite these common patterns, the spatiotemporal patterns of terrestrial productivity IAV across Europe derived from various proxy data sets differ to some extent (Figs. 1, 3, Supplement Figs.S1, S2).In particular, the IAV of summer crop yield are not consistent with the IAV in the other productivity proxies.Summer crops are expected to be strongly controlled by intensive human management (e.g.irrigation), whose effects on IAV of crop yields cannot be completely eliminated even though we filtered our analysis in less intensively irrigated areas, and different phenological development stages compared to natural or less intensively managed vegetation.They have also a shorter growing season, which makes crop yield sensitive to climate anomalies only during short periods.
The IAV of terrestrial productivity proxies in temperate Europe appears to be more sensitive to variations of extreme temperature conditions (such as MMT), rather than to variations in mean temperature (Table 1).This conclusion is confirmed by increasing evidence from ground-truth observations and remote sensing (Babst et al., 2012;Ciais et al., 2005;Reichstein et al., 2007).Importantly, there is an increasing trend in the IAV of productivity proxies (except summer crop yields) especially in water-limited regions in   southern and eastern Europe (Fig. 8, Supplement Figs.S14,  S18).The underlying mechanisms for the increasing IAV of productivity in such regions remain unresolved but could at least partially be attributed to increasing variations in climate (Supplement Fig. S21) and associated changes in vegetation phenology (e.g.Maignan et al., 2008;Zhang et al., 2006).One possible reason for increased IAV of spring crops could be linked to the general average yield increase of these crops over the past decades, i.e. a larger average yield may exhibit a larger risk of being negatively affected by adverse climate conditions in a given year.The variations in vegetation carbon uptake and phenology induced by changing climate play an important role in regulating interannual variability of the growing season's averaged productivity (Wu et al., 2013).Remarkably, there seems to be increased asynchrony between water availability and growing-season length in water-limited regions (Austin et al., 2004).Hence, subtle  changes in water balance are expected to exert great effects on the productivity in such regions.This finding suggests that the productivity is becoming more vulnerable to the changing climate IAV, especially in water-limited regions.Inference of how a more vulnerable productivity could or could not translate into a lesser carbon sink remains speculative.However, if the productivity-sink function is concave, a more vulnerable productivity could end up decreasing the carbon sink.Additionally, for crops that are harvested each year, a more vulnerable productivity could have negative economic impacts but only minor effects on the soil carbon balance, which is controlled by other factors such as tillage, and nongrowing-season climate conditions.Furthermore, the regionspecific responses of productivity-related proxies to climate variations imply that changes in global climate zones and hence vegetation distributions could indirectly alter the relationships between climate variability and carbon balance.This could be expressed, e.g. as a northward extension of water-limited areas.
Notably, climate variations can only explain 20-40 % of the spatial variance in the IAV of the analysed proxies in most of temperate Europe (Table 1).Other factors, such as soil water holding capacity, lagged effects, human management (irrigation, fertilization), specific land cover patterns and climate regimes, likely contribute significantly to the spatial IAV of terrestrial productivity (Ciais et al., 2010;Gervois et al., 2008).However, the magnitude and individual contribution of these processes to the IAV of the carbon cycle remain largely unknown.

Changing sensitivity of productivity proxies to climate variability in Europe
The changing sensitivity of terrestrial productivity proxies to climate is expected to affect the carbon balance and associated climate feedbacks (Cox et al., 2013).We observed a consistently significant negative γ IAV of the spring's crop yield, NDVI, FAPAR and GPP to mean growing-season temperature and positive γ IAV to total growing-season WAI, in temperate continental and Mediterranean Europe (Fig. 7, Supplement Fig. S12), indicating that the terrestrial productivity in these regions is sensitive to interannual variations of the water balance.The change of γ IAV of terrestrial productivity is ambiguous in western Europe and shows different patterns compared to other parts of Europe (Fig. 7, Supple-ment Fig. S12).This is probably linked to the oceanic climate regime and the lack of a single dominant limitation on productivity.This finding is supported by previous studies reporting a divergent response of western European vegetation to climatic fluctuations during the past decade compared with the other parts of Europe (Zhao and Running, 2010).Overall, the IAV of productivity in the water-limited southern and eastern regions in Europe, which are mainly dominated by short-rooted grasslands and croplands, show a higher sensitivity to water availability than in other regions which are less water-limited (Reichstein et al., 2007;Schwalm et al., 2010).Interestingly, we observed an increasing sensitivity of productivity to water availability in temperate continental and Mediterranean Europe (Fig. 8,, which may be attributed to the recently increased IAV of water availability in these regions (Supplement Fig. S21) due to the non-linear threshold-like response of productivity to water availability and to strong land-atmosphere coupling (Seneviratne et al., 2006)   The trends in time series are estimated by the Theil-Sen slope method (Yue et al., 2002).Significant trends were marked by black points.Blank region indicates that there is intensive irrigation (with irrigated area fraction ≥ 15 %), no crop yield records or a very low (< 5 %) fraction of cropland in those pixels.terrestrial productivity in these regions is becoming more vulnerable to changes in water conditions.Notably, simulations predict a considerable enhancement of IAV of European climate, associated with higher risks of heat waves and droughts throughout the 21st century (Meehl and Tebaldi, 2004;Schär et al., 2004), suggesting that climate-driven productivity variability will continue to increase in Europe.Evapotranspiration tends to increase with climate warming and will likely lead to a more negative water balance  in dry regions or to an extension of the total water-limited area if precipitation does not increase concurrently (Supplement Fig. S20) (Christensen et al., 2007;Kjellström et al., 2011).The increasing sensitivity of productivity to climate IAV accompanied by an increasing water deficit (Supplement Fig. S20) and changes in plant phenology could lead to a weakening trend in productivity and a weaker carbon sink in this region (Hu et al., 2010;Rivier et al., 2010;Wu et al., 2013).Assimilation is more drought sensitive than respiration and could further weaken the terrestrial carbon sink in such regions (Schwalm et al., 2010), or even turn them into net carbon sources to the atmosphere.

Conclusions
Evidence from different proxies indicates that the IAV of terrestrial productivity in Europe is spatially heterogeneous and much higher in water-limited regions.Despite considerable regional differences, the IAV of terrestrial productivity and its sensitivity to climate variation showed pronounced temporal changes in Europe during past decades.In water-limited regions there has been an increase in the IAV of terrestrial productivity and in the sensitivity of terrestrial  productivity to climate variations, which is at least partly attributable to the increased IAV of water availability in these regions.These findings emphasize that the carbon cycle in water-limited regions in Europe is becoming more vulnerable to changes in water availability.The increasing climate sensitivity of carbon uptake accompanied by the predicted increase in water deficit could lead to weaker carbon sinks.Although the terrestrial productivity in Europe responds directly to the changing IAV of climate, it could also be co-regulated by other factors, such as climate regimes, land use patterns and human footprints.We anticipate that atmospheric inversion and land-surface modelling will be helpful to further attribute changes in the climate sensitivity of productivity IAV.Similarly, further development of proxy/historical data sets (e.g.tree-ring indices, remotesensing observations) should be increasingly regarded as a unique and valuable empirical baseline of carbon cycle processes over the past centuries.
The Supplement related to this article is available online at doi:10.5194/bg-11-3057-2014-supplement.

Figure 1 .
Figure 1.Spatial patterns of CV of detrended crop yield data during the 1975-2009 period for barley (a), wheat (b), grain maize (c) and potatoes (d).Blank region indicates that there is intensive irrigation (with irrigated area fraction ≥ 15 %), no crop yield records or a very low (< 5 %) fraction of croplands in those pixels.

Figure 2 .
Figure 2. Spatial distribution of Spearman correlation coefficients between the IAV of barley yield and IAV of the corresponding growing season's (March-June for barley) summed NDVI (a), GPP (b), and FAPAR (c).Significant relationships are marked by black points.Blank region indicates that there is intensive irrigation (with irrigated area fraction ≥ 15 %), no crop yield records or a very low (< 5 %) fraction of croplands in those pixels.

Figure 3 .
Figure 3. Relationships between barley yield IAV and IAV of other proxies.Scatter plot of the relationships between the correlation coefficients of IAV of barley yield and IAV of the growing season's (March-June for barley) summed NDVI, GPP, and FAPAR.Different colours indicate different Köppen-Geiger climate zones.

Figure 4 .
Figure 4. Frequency distributions of correlation coefficients of IAV of barley yield and IAV of the growing season's (March-June) summed NDVI, GPP, and FAPAR in Europe (green bars) and climate zones: Cf (yellow bars), Cs (light blue bars), and Df (brown bars).Percentage of pixels where there are significant (p < 0.1) relationships between IAV of barley yield and IAV of the growing season's summed NDVI, GPP, and FAPAR are shown in the inset figures.

Figure 5 .
Figure 5. Spatial patterns of the trends in crop yield interannual variability.The estimated trends in the moving CVs of the detrended crop yield calculated using a 10 yr sliding window with 1 yr lag during the 1975-2009 period for barley (a), wheat (b), grain maize (c) and potatoes (d).The trends in time series are estimated by the Theil-Sen slope method(Yue et al., 2002).Significant trends were marked by black points.Blank region indicates that there is intensive irrigation (with irrigated area fraction ≥ 15 %), no crop yield records or a very low (< 5 %) fraction of cropland in those pixels.

Figure 9 .
Figure 9. Probability density function for the differences of moving climate sensitivities of barley yield to mean growing-season temperature (blue line) and total growing-season water availability index (red line) between 1995-2008 and 1982-1995 for Europe (a), and climate zones Cf (b), Cs (c), and Df (d).The inset figure shows the percentages of pixels suffering the four different trend types (significantly increasing, IS; increasing but not significant, INS; decreasing but not significant, DNS; significantly decreasing, DS) for climate sensitivities of barley yield to mean growing-season temperature (blue bars) and total growing-season water availability index (red bars).

Wu et al.: Climate-mediated variability in terrestrial productivity across Europe calculating
(Kottek et al., 2006)respective time series using a 10 yr window with 1 yr lag during the 1975-2009 period (if applicable).The choice of the moving window length is a compromise between a large sample size within each window and a sufficient number of windows to evaluate changing IAV.The sign and magnitude of the monotonic trends in moving CVs for the productivity proxies in each grid are estimated by the Theil-Sen slope method, which can provide an accurate estimation of trends in time series with strong autocorrelation (in detail seeYue et al., 2002).Differences in spatiotemporal patterns in the IAV of each productivity proxy among different Köppen-Geiger climate zones are investigated(Kottek et al., 2006).We consider mainly three dominant climate zones in our study including the warm temperate humid zone (Cf, grouping of Cfa, Cfb and Cfc), the warm temperate arid zone (Cs, cluster of Csa and Csb) and the snow humid zone (Df, grouping of Dfa, Dfb and Dfc).

Table 1 .
Generalized linear model analyses for the relationships between productivity proxies and climate in different climate zones.Climate zones used in this study are based on the Köppen-Geiger climate classification.c GLM: generalized linear model.TEM, WAI, CFD, and MMT are mean growing-season temperature, total growing-season water availability index, total spring (March-May) consecutive frost days and mean growing-season maximum of daily maximum temperature, respectively.The growing seasons are different for different crops.We fixed the growing season as March-June for barley and wheat and as March-September for potato and grain maize.d The scaling factor for the GLM coefficients of CFD is 10 −3 . b