Bivariate return periods of temperature and precipitation explain a large fraction of European crop yields

This study seeks the linkage between multivariate climate conditions and crop yields. This paper is one of the first which employs bivariate return periods of temperature and precipitation as the indicator of climate variability to explain crop yield variability. It is clearly written and obtains the interesting finding that the combination of temperature and precipitation can explain more crop yield variability than models relying directly on temperature and precipitation as predictors on average in Europe. The result also reveals different sensitivities of crops to climate conditions. A need to incorporate the nonlinear impacts into the climate-crop yield assessment is highlighted. For all these reasons, I recommend publication after addressing a few comments regarding the statistical examinations.

The authors examined 6 types of copulas in order to represent different combined effects of temperature and precipitation, i.e. dry and hot, dry and cold, wet and hot, wet and cold. My understanding is that the impact of these combinations can be due to a single variable or both variables being in an extreme state. The examined copulas, however, do not include extreme-value copulas, which are usually considered more appropriate for reproducing the interrelationship/interdependence structure between rare events. Did the authors compare the result using extreme-value copulas?
The t-copula and some of the Archimedean copulas allow tail dependence, i.e. dependence in the extremes, see e.g. Schoelzel and Friederichs (2008) Sections 4.1.2 and 4.2. The t-copula is symmetric such that if the upper tails are dependent, the lower tails are as well. In contrast, the Gumbel and Joe copulas can model only upper tail dependence whereas the Clayton copula can model lower tail dependence (all three are asymmetric). Note that in this application we aim for modeling the whole distribution to also capture events with high return frequency and not only the extremes. Hence, only using extreme value copulas might not be appropriate. We have added a note in Section 2.4 to highlight that the copulas that we use are able to capture dependence in the extremes. "The t copula and some of the Archimedean copulas allow tail dependence, that is, dependence in the extremes, (see e.g. Schoelzel and Friederichs, (2008)). The t copula is symmetric such that if the upper tails are dependent, the lower tails are as well. In contrast, the Gumbel and Joe copulas can model only upper tail dependence, whereas the Clayton copula can model lower tail dependence (all three are asymmetric). Note that in this application, we aim for modeling the whole distribution to also capture events with low return periods and not only the extremes. Hence, only using extreme value copulas might not be appropriate." Here the authors applied copulas to the seasonal and 2-month averaged climate variables. One of the prerequisites to apply copula is the assumption of temporal independence of variables, e.g. by examining the autocorrelation. Did the input variables meet this requirement?
Thank you for this comment. We have now tested for autocorrelation in the yearly time series of temperature and precipitation averaged over the different time scales use in the study (66 years). 8.8% of the temperature time series and 4.4% of the precipitation time series are significantly autocorrelated for lag 1 at the 5% level. Since we look at bivariate distributions, we conclude that autocorrelation should not have a significant effect on the conclusions of the paper. We have add this information in the revised manuscript in Section 2.2: "We tested for autocorrelation in the yearly time series of temperature and precipitation averaged over these different time scales. 8.8% of the temperature time series and 4.4% of the precipitation time series are significantly autocorrelated for lag 1 at the 5% level. Because our focus is on bivariate distributions, we conclude that autocorrelation does not have a significant effect on our conclusions."

Introduction
Agriculture is essential to the well-being of humans and is directly affected by changes in climate (Lobell et al., 2011;Porter et al., 2014). Continuing climate change will likely increase the pressure on agriculture in the future with adverse impacts on food security (Porter et al., 2014;Wheeler and Von Braun, 2013;Rosenzweig et al., 2014). Understanding how crop yields vary with climate may help to ensure future food security through increased predictability of yields and adequate adaptation 5 measures. Trying to keep variability of yields low is a key objective as lower yield variability leads to more stable farmer incomes (Reidsma et al., 2010) and food supply (Slingo et al., 2005). Biological crop models exist at least since the 1970s (Loomis et al., 1979). However, such models are limited in their ability to quantify the impact of climate variability on crop yields over larger scales. For this purpose statistical models are applied. Globally, variations in precipitation and temperature account for about a third of the observed yield variability (Lobell and Field, 2007;Ray et al., 2015). Due to their sensitivity to 10 climate, crop yields are also strongly susceptible to climate extremes. In particular, disastrous droughts and heat waves often severely impact crop yields, reducing national cereal production by up to 9-10% (Lesk et al., 2016). How climate extremes impact crop production, however, varies across locations. Crops in the northern latitudes, for instance, might thrive under very warm conditions. In addition, impacts of climate extremes on vegetation strongly depends on when the climate extreme occurs and in which developing stage the underlying vegetation is (De Boeck et al., 2011;Frank et al., 2015;Sippel et al., 2016). 15 Because crops are often strongly adapted to specific climate conditions and covary with local climate (Osborne et al., 2009), it can be expected that the timing of climate anomalies is even more crucial for crops.
Crop yields can respond nonlinearly to changes in climate conditions and extremes (Porter and Semenov, 2005). Nonlinear impacts of extremely dry and hot conditions have been shown previously on terrestrial carbon fluxes. At the global scale, the impacts of concurrent dry and hot periods on carbon fluxes exceed the sum of impacts from either hot or dry conditions 20 (Zscheischler et al., 2014a). Studies investigating the variability of crop yields with climate usually rely on linear models (Osborne and Wheeler, 2013;Ray et al., 2015) though sometimes nonlinear transformations of climate variables are used as predictors (Ray et al., 2015). Building more complex nonlinear models is often not feasible because crop yield time series are too short (mostly not more then 60 years). However, climate extremes might have disproportionately large impacts on crops (Lesk et al., 2016). Due to the sensitivity of crop yields to both temperature and precipitation (Ray et al., 2015) and extremes 25 therein (Lesk et al., 2016), a variable incorporating the degree of extremeness of both variables at the same time might be more robustly related to crop yields. In the univariate setting, the magnitude of an extreme event is sometimes approximated by its return period. Here we use this analogy and derive bivariate return periods of temperature and precipitation which we then relate to crop yields. We explore how well bivariate return periods of temperature and precipitation are linked to crop yields.
Multivariate return periods have been studied increasingly in the recent past, mostly in the field of hydrology, for instance to 30 study floods (Grimaldi and Serinaldi, 2006;Salvadori et al., 2011Salvadori et al., , 2013 and droughts (Shiau and Modarres, 2009;Wong et al., 2009;De Michele et al., 2013). The introduction of the concept of copulas into meteorology and climate research has boosted the usage and modeling of multivariate distributions (Schoelzel and Friederichs, 2008). A copula is a multivariate probability distributions for which the marginal probability distribution of each variable is uniform. Copulas can be used to describe the 2 dependence between random variables. They can greatly simplify calculations involving multivariate distributions and have lead to suite of definitions of multivariate return periods (Salvadori and De Michele, 2004;Serinaldi, 2015). Most studies so far have focused on return periods of droughts and floods in a multivariate setting. AghaKouchak et al. (2014) were probably the first to compute return periods of events involving both precipitation and temperature.
Generally, there is no unique way to define multivariate return periods (Serinaldi, 2015). As there is no natural ordering 5 in higher dimensional spaces, one has to decide first in what direction one wants to "look". For instance in our case, we can compute return periods of hot and dry events, hot and wet events, cold and dry events, or cold and wet events. Furthermore, we have to decide how to set a threshold to compute threshold exceedances which will ultimately result in return periods of specific events. Various suggestions have been proposed in the recent years (Serinaldi, 2015). Here we use the definition motivated by the notion of "critical regions" (Salvadori et al., 2011) whose return periods can be computed using so-called survival functions 10 and survival copulas . ::: The ::::::: concept :: of :::::: critical :::::: regions :::: has ::::::: recently :::: been :::: well ::::::::::::: mathematically ::::::::: formalized :: via :::: the :::::: notions :: of ::::::: Hazard :::::::: Scenarios :::: and ::::: Upper ::::: Sets, :::::::: providing :: a ::::::::: consistent ::::::::::: mathematical :::::::::: framework. :::: The :::::::: interested :::::: reader : is ::::::: referred ::: to ::::::::::::::::::: Salvadori et al. (2016). : A survival function is the complementary cumulative distribution function, i.e., for a distribution function F (x) = Pr(X ≤ x) the survival function, denoted byF , is given byF (x) = Pr(X > x). Survival copulas are copulas based on survival functions. 15 Critical regions are separated from noncritical regions by a so-called critical layer {(x, y) ∈ R 2 :F (x, y) = Pr(X > x, Y > y) = t}, where all points (x, y) ∈ R 2 have the same probabilityF (x, y) = t in the joint survival distributionF . The critical region is then defined asR < t = {(x, y) ∈ R 2 :F (x, y) < t} (De Michele et al., 2013). A return period RP of an event can be defined as the inverse of the probability of falling in the critical region, where µ is the average interarrival time of events (1 year in our case). As we will see, from this definition it follows that events in which only one variable is extreme (e.g., either extreme heat but no drought, or extreme drought but no heat) and events in which both variables are moderately extreme (e.g., moderate heat and drought) may have similar return periods. Crop yield vary with climate, most importantly temperature and precipitation (Ray et al., 2015). If we assume that temperature and precipitation at potentially different time periods are of similar relevance for crop yields, a relationship between bivariate return 25 periods of precipitation and temperature can be expected. Here we apply this approach to climate data in Europe and investigate relationships between bivariate return periods of combinations of temperature and precipitation with crop yield variability. By "looking in all directions" in 30 the temperature-precipitation space we assess whether an intensification along a certain direction (dry and hot, dry and wet etc.) leads to an increase or decrease of crop yields. By computing return periods for different combinations of months (e.g. temperature in spring, precipitation in June and July) we estimate which combinations of months and climate conditions mostly affect crop yields. We compare this analysis with linear models fitted to the two predictors precipitation and temperature.  (Leff et al., 2004) and good spatiotemporal coverage, we selected the following crops (following the 5 official EUROSTATS terminology): Cereals (excluding rice), Wheat (including spelt), Maize (grain maize and corn-cob mix), Potatoes (including early potatoes and seed potatoes), and Sugar beet (excluding seed). Wheat does not include winter and spring varieties. Sugar beet is sown in spring and harvested in autumn in Europe. In the EUROSTAT database, yield is reported as the amount of dry matter suitable for consumption (Moors et al., 2010). Not every country grows every crop every year, accordingly some time series contain missing values. In total we collected 3552 crop yield years distributed over 19 countries 10 and 5 crops. Most data was available for cereals (781 crop yield years) and potatoes (766 crop yield years) followed by wheat (734), sugar beet (672), and maize (599 crop yield years).
We further used daily minimum (T min ), mean (T ) and maximum temperature (T max ), and precipitation (P ) data from E-OBS (Haylock et al., 2008) at 0.25 • spatial resolution and computed spatial averages for the 19 countries for which we have collected crop yields. Because crop area potentially changed over time, for each country we averaged climate data over the 15 entire country. In this way we investigate how climate variability at country scale is related to annual crop yields.

Data preprocessing
Long-term trends in crop yields are mostly caused by technical progress such as breeding, and changes in land use policy and management practices. Various corrections are typically used to adjust for this effect. Often the difference in values from one year to the next (Lobell and Field, 2007), or linear or quadratic detrending (Ray et al., 2015) are applied to obtain anomalies of 20 crop yields independent of long-term trends which are more closely related to climate variability. Here we rely on a different, more adaptive approach (Wu et al., 2014). We fitted cubic smoothing splines to each crop yield time series to capture the long-term trend. The original time series was then divided by this trend to obtain a dimensionless index. Values greater than 1 indicate positive yield anomalies whereas values below 1 indicate negative yield anomalies. The cubic smoothing splines were fitted with a frequency cutoff of 50% at 2/3 of the time series length to remove low frequency variability. This procedure 25 is analogous to standard dendrochronological procedures, where the aim is to filter out low-frequency trends related to tree geometry and not to climate (Babst et al., 2012;Cook and Peters, 1981). We used the R package dplR (Bunn, 2008) for these computations and denote the resulting yield anomalies by Y . We used the same normalization to subtract long-term trends in climate data (T min , T , T max , and P ) for each country, similar to (Wu et al., 2014). The resulting indices preserve high-frequency (i.e., interannual) variability, but not the trends or low-frequency signals (Babst et al., 2012;Cook and Peters, 30 1997).

Bivariate return periods
The computation of bivariate return periods was done based on copulas similar to (AghaKouchak et al., 2014). The joint 10 distribution of two variables X (e.g., temperature) and Y (e.g., precipitation) can be written as (Sklar, 1959;Salvadori and De Michele, 2004) :::::::::::::::::: F X (x) = Pr(X ≤ x), :::::::::::::::: and a copula C. Copulas are multivariate probability distributions with uniform marginal distributions designed to model the 15 dependence between multiple variables (Nelsen, 1999). As we are interested in exceedance probabilities, we rely on survival functions. The so-called joint survival distributionF (x, y) = Pr(X > x, Y > y) can be obtained using the concept of survival copulas (Salvadori et al., 2011 F with marginal survival functionsF X = 1 − F X andF Y = 1 − F Y , andĈ the survival copula. For any given (x, y) ∈ R 2 , 20 there exists a unique survival critical layer on which the set of realizations of X and Y share the same probability t ∈ (0, 1) (De Michele et al., 2013;AghaKouchak et al., 2014):LF t = {(x, y) ∈ R 2 :F (x, y) = t}. The survival critical layer partitions R 2 into a "safe" and a "critical" regionR < t = {(x, y) ∈ R 2 :F (x, y) < t}. A bivariate return period RP can now be defined as the probability of falling inside the critical regionR < t De Michele et al., 2013) through
For a given temperature-precipitation couple (x, y) we can now (i) determine the survival critical layerLF t associated with t =F (x, y) =Ĉ(F X (x),F Y (y)), (ii) compute the probability of an event being at least as extreme as (x, y) (i.e., falling inside the critical region defined by t) by p = 1 −K(t), and finally (iii) compute return periods via (4).

Results
We first present three examples for particular country-crop combinations illustrating our approach, and the comparison between

Examples
Yields often vary along a specific precipitation-temperature gradient. In figure 1 we present an example of our approach for maize yield anomalies in Germany. Maize yield anomalies are high when maximum temperatures in July and August are low and precipitation in June and July is high (lower left of figure 1a). Yields decrease with increasing temperatures and decreasing precipitation. Return periods of hot and dry events increase along the same gradient, as denoted by the curves of equal return periods. Accordingly, bivariate return periods along this gradient capture the variability of maize yield anomalies well, explaining 72% of their variance (figure 1b). This variability is less well captured if a linear model is directly fitted to the same predictors (M T P , R 2 = 0.62, figure 1c).
In a similar fashion, figure 2 presents an example of cereal yields in Lithuania. In this case, yields are high when maximum temperatures are low in June and July, and when precipitation is low in August and September (lower left of figure 2a).

10
They decrease with increasing temperatures and, in contrast to our first example, increasing precipitation most likely because Lithuania's climate is already wet (it has humid continental climate with an average annual rainfall of around 620 mm). Return periods of hot and wet conditions increase along the same gradient and thus they capture the variability of cereal yield anomalies well (R 2 = 0.56, figure 2b). Also in this example, the variability in yield anomalies is less well captured if a linear model is fitted to the same predictors (M T P , R 2 = 0.47, figure 2c).

Summary statistics 25
The models M rp perform significantly better than M T P (p < 0.05, figure 4). This is also evident for the individual cropcountry combinations (figure 5). On average the fraction of explained variance is 42% for M rp (colored bars) versus 36% for M T P (white bars). This corresponds to an increase in predictability of 17%. Crop yield anomalies mostly vary strongest along a hot-dry to cold-wet gradient with low yields during hot and dry conditions (41 out of 91, 45%, red bars). Crop-country combinations where yield anomalies increase along the hot-wet to cold-dry gradient (green bars) represent 29 of all 91 cases 30 (about 32%). The other two gradients appear less often, cold-wet to hot-dry 10 times (blue bars) and cold-dry to hot-wet 11 times (purple) out of 91.  Notably, specific precipitation sums over longer time periods (mostly 3 months and more) seems to be most relevant for crops.

5
In contrast, crop yields are sensitive to temperature conditions at shorter time scales (generally less than 3 months). Precipitation averaged over the whole spring and summer (13 incidences) or over one season (28 incidences) was selected in nearly 41% of the cases while averaging periods for temperatures were only 2 months in nearly 86% of the cases (78 incidences out of 91). In total, for precipitation 2.5 times more months are selected than for temperature (524 versus 211). Overall, climate conditions in summer (JJA) seem to be more relevant for crops and climate during other periods. For both temperature and precipitation, 10 63% of the selected months fall in the summer.

Discussion and Conclusions
Linear models based on return periods of specific climate conditions (M rp ) can better explain variability in annual crop yield anomalies than linear models based directly on precipitation and temperature (M T P , 42% versus 36% of explained variance in crop yield anomalies on average, p < 0.05). The reason for this might be that yields react nonlinearly to more extreme climate  13 conditions. Including more variables such as radiation, humidity, and soil moisture might further improve the prediction of crop yield anomalies. However, due to the shortness of yield time series (more than 50% of the yield series contain less than 30 years of data), robust models with many predictors are difficult to build. Our analysis demonstrates the different sensitivities of crops to climate conditions. Figures 5 and 6 can be used in combination to disentangle the sensitivities of crop yields to climate conditions at certain time periods. For example, in most countries Maize yields decrease for hot and dry conditions. However, 5 in countries such as Lithuania, Luxembourg and Great Britain, maize yields increase under hot and wet conditions ( figure 5).
For all three countries, maize yields are most sensitive to precipitation in summer ( figure 6). Overall, crops are most sensitive to climate variations along the hot-dry to cold-wet gradient. This in line with what can be expected from univariate assessments (Lesk et al., 2016), showing that crop yields can decrease substantially under heat or drought. Further, crops are most sensitive to maximum temperatures (in contrast to mean and minimum temperatures) whereas the most important time period largely 10 depends on the country and the crop. Because crops can exhibit a thresholds behavior at high temperatures resulting in crop damage (Porter and Semenov, 2005), maximum temperatures might be the most crucial variable for yield variability in most countries and crops.

Uncertainties in the estimation procedure
Estimating return periods is not free of uncertainties. Bivariate copulas are defined on the uniform marginal distributions, which 15 is achieved by ranking of the original marginal distributions. Copulas are thus related to the original bivariate distribution by (3). Naturally, the accuracy of this transformation improves with increasing sample size. After the data have been transformed to the unit square, copulas have to be fitted. The fitted copulas we used passed standard goodness-of-fit tests though longer climate time series certainly help to obtain more robust fits (E-OBS data currently covers . Uncertainties in the fitting procedure may be particularly large for extreme return periods (Serinaldi, 2015).

Return periods of climate events and their relation to impact-related variables
Despite the above-mentioned uncertainties, which are particularly relevant for extreme return periods, we argue that bivariate return periods of temperature and precipitation conditions can be used as a robust estimate of the climate variation along a specific temperature-precipitation gradient, outperforming conventional linear models. In particular if both temperature and precipitation are strongly affecting crop yields, bivariate return periods will be related to yield variation. This assumption holds 25 for many crops and regions in Europe (Peltonen-Sainio et al., 2010) but might be different for other regions in the world (Ray et al., 2015). The gradient which is most strongly related to crop yield variability or other impact variables can be chosen by testing different directions in the temperature-precipitation space. In addition, due to the nonlinear nature of return periods, extreme conditions might be captured more appropriately. In particular, extreme heat and drought can have extreme impacts on crops which might exceed the impacts predicted by linear models (Porter and Semenov, 2005;Zscheischler et al., 2014a;30 Lesk et al., 2016).
For most country-crop combinations, crops are sensitive to climate along the hot-dry to cold-wet gradient, with hotter and drier conditions leading to lower yields, though there is large variability among countries (figure 5). This is in line with previous research on the relationship between climate and other ecosystem variables. For instance, photosynthetic carbon uptake also generally increases along the hot-dry to cold wet gradient (Zscheischler et al., 2014a;Ahlström et al., 2015). Furthermore, most negative impacts on terrestrial carbon uptake are associated with hot and dry conditions (Zscheischler et al., 2014b).
We have demonstrated how bivariate return periods of climate conditions as an indicator for climate variability along a certain climate gradient can be used to predict crop yield anomalies. Our proposed approach explains significantly more yield 5 variability than linear models based on temperature and precipitation. However, climate time series should be long enough to robustly estimate bivariate return periods (preferably at least 50 years). Crops largely rely on precipitation and temperature at the right time periods, depending on the specific crop and the location. Our analysis illustrated that crops are often most sensitive to summer climate, maximum temperatures, and specific precipitation conditions over longer time periods (at least 3 months). Other ecological variables sensitive to impacts of climate variability such as tree rings and terrestrial carbon uptake 10 and release may be analyzed in a similar fashion. If other driver variables despite precipitation and temperature are crucially relevant for an impact variable (such as for instance radiation, humidity, and soil moisture) extensions to higher-dimensional return periods are possible, though the choice of the direction of a gradient becomes more difficult in higher dimensions. In conclusion, bivariate return periods allow to summarize bivariate climate conditions in a novel univariate measure which then can be used more easily to investigate associated impacts.