Wheat yield loss attributable to heat waves, drought and water excess at the global, national and subnational scales

Wheat, with about a 2.1 million km² total harvested area, is the most abundant crop in the world: it is the first rain-fed crop after maize and the second irrigated crop after rice (Portmann et al 2010). With a total production that surpassed 700 million tons (MTons) in year 2010, it is contributing to about the 20% of the total dietary calories and proteins worldwide (Lobell and Gourdji 2012, Shiferaw et al 2013).

Field studies of wheat yield variability have shown that higher than optimal temperatures during the growing season are generally accelerating the progress of the plant phenological stages, affecting photosynthesis and its balance with respiration, and reducing the final yield (Lobell and Gourdji 2012, Rezaei et al 2015). Higher temperatures also increase the atmospheric demand for water and reduce the crop water-use efficiency (Ray et al 2002). Exposure to extremely high temperatures (i.e. heat stress) leads to plant damages by inducing perturbations in cellular structures and metabolic processes (Nakamoto and Hiyama 1999). Isolated occurrences of extreme high temperatures around a sensitive stage of crop development, such as flowering and grain-filling, can reduce grain yield considerably (Tashiro and Wardlaw 1990, Ferris et al 1998, Porter and Gawith 1999, Luo 2011), while a prolonged period of extreme high temperatures might result in almost total yield loss (Semenov and Shewry 2011).

The detrimental effect of the heat stress on wheat yield may worsen when coinciding with drought (Pradhan et al 2012). Irrigation, to alleviate local water stress, can also mitigate the impact of the highest temperatures, by evaporative cooling, even at larger spatial scales (Tack et al 2015, Troy et al 2015, Mueller et al 2016). On the other hand, extreme amounts of precipitation and water excess in the soil can also be responsible for wheat loss due to proliferation of pests and diseases, leakage of nutrients, inhibition of oxygen uptake by roots, and interference with agronomical practices (e.g waterlogging during harvest). This may have happened during the 2016 cropping season in France, which showed a huge wheat yield anomaly that has not been officially confirmed yet (MARS 2016, Bloomberg 2016, France24 2016, Arvalis 2016).

Compared to other important crops, the main wheat producing regions are characterized by 'close-to-average' yield variability (Ben-Ari and Makowski 2014). However, increasing unfavorable conditions under observed and projected climate change conditions will probably impact wheat production variability (Gourdji et al 2013, Deryng et al 2014, Siebert and Ewert 2014, Asseng et al 2015, Ray et al 2015), especially because of the exacerbated effects of heat stress on grain number and grain size (Lobell et al 2015).

Using an ensemble of crop model simulations, Lobell et al (2015) estimated for the period 1980–2008 a 5.5% reduction of global wheat production due to changes in temperature and precipitation. This reduction occurred together with an increase in wheat yield/production of about 50% over the same period, mostly due to improved management and higher yielding crop varieties. Therefore, considering that the human population increased proportionally from 4.4 billion in 1980 to 6.9 billion in 2010 (UN 2015), negative yield deviations from the mean trend can be a potentially increasingly threat to food security. Indeed, changes in wheat yield variability have already been suggested as one of the primary factors influencing global food prices, market stability and food security, especially in developing countries (FAO 2011, OXFAM 2012, Porter et al 2014, Deryng et al 2014, GFS 2015). Moreover, the current urbanization trend at the expense of cropland area extent would imply, to be sustainable, a further increase of crop yield (Bren d'Amour et al 2016), and presumably an increased sensitivity of food security to the climate induced yield interannual variability.

The supporting evidence for these assertions is still relatively weak, since only few recent studies have addressed the relationships of the observed inter-annual climate and wheat yield fluctuations at the global scale. Lizumi and Ramankutty (2016) found consistent increases in standard deviations of annual yields and of an agro-climatic indicator accounting for solar radiation, heat stress and drought. Ray et al (2015), through an automatic spatial selection of several agro-climatic indicators based on temperature and precipitation anomalies (and using high resolution yield data), developed a statistical model accounting for the sign of the fluctuations and explaining about one third of the observed yield variance. This is in line with the composite analysis by Lesk et al (2016), based on an a priori knowledge of extreme events, showing that droughts and extreme heat could have reduced national annual yields by 9%–10%. No effects of floods and extreme cold events were identified by Lesk et al (2016) on national data.

It is worth noting that none of aforementioned studies explicitly considers the possible impact on wheat yield of water excess events, which are not necessarily only related to floods. Moreover, the indicator used by Lizumi and Ramankutty (2016) accounts more for sub-optimal climatic conditions, rather than for the extremes. On the other hand, the indicators used by Ray et al (2015) are computed using the seasonal mean variables, which are not capturing the extremes occurring on a sub-monthly scale as the heat waves. Both of them are based on the exceedance of prescribed thresholds that are not explicitly mentioned. This might not be justified at the global scale, because of the different climates and different agronomical practices (such as irrigation and variety selection).

Here, we aim to model explicitly the effects of temperature and soil moisture anomalies (either positive and negative) on global and regional wheat production anomalies. The proposed approach does not require an a priori knowledge of extreme events (as in Lesk et al 2016) and it has the clear advantage of disentangling the different contributions of heat and water anomalies, potentially resulting in stress for the crop. Furthermore, it helps to better identify and characterize the role of the major producers in driving global yield anomalies, showing explicitly the effects of the individual climatic anomalies and extremes for each year in the period 1980–2010. Compared to Ray et al (2015) and Lizumi and Ramankutty (2016), we use the most up-to-date heat waves and water balance indicators currently available. These are, namely, the Heat Wave Magnitude Index daily (HWMId; Russo et al 2015) and the Standardized Precipitation Evaporation Index (SPEI; Vicente-Serrano et al 2010, 2013). Compared to previous studies, both indexes have been developed to be generally applicable under a wide range of conditions. They bear the great advantage of allowing the diagnosis and the comparison of extreme events in different regions of the world and time of the year on the basis of globally defined thresholds (Zampieri et al 2016, Vicente-Serrano et al 2010, 2013, see Data and Methods section).

In this study, the HWMId has been slightly modified in order to account for temperature anomalies and heat waves events during the stages of wheat growth that are sensitive to heat stress (i.e. three months before harvesting, roughly corresponding to the period including anthesis and grain filling). Therefore, the new indicator bears some analogy with the Growing Degree Day (DGG) index, but instead of summing the daily temperature anomaly w.r.t. a base value, it sums up the heat wave magnitude function defined by Russo et al 2015 (see equation (1) in the Data and Methods section). We refer to it as Heat Magnitude Day (HMD). While the HWMId selects the maximum value of the summations computed on consecutive periods with maximum daily temperature over a certain threshold (i.e. a heat wave), the HMD integrates all of them similarly to the GDD. Therefore, the HMD is able to capture not only heat stress events, but also sub-optimal conditions.

The SPEI is a multi-temporal-scale index that quantifies persistent anomalies in the soil moisture balance (i.e. precipitation minus potential evaporation) over different time- periods (1–24 months). We use the SPEI6 computed during the previous six months before harvesting, and we refer to it simply as SPEI. The SPEI6 provides an estimate of the soil moisture state during the sensitive stage of wheat, as for the heat stress, but also accounting for a preconditioning period of three months before flowering consistently with the time-scale of the root soil moisture dynamics (see e.g. Zampieri et al 2009). Large negative values of the SPEI indicate drought, while positive SPEI indicate wetter-than-normal soil conditions.

In order to account for concurrent heat and water stress events, we define a combined index as a simple linear superposition of the standardized HMD and SPEI (Combined Stress Index, CSI). We calibrate the CSI by determining the coefficients that maximize the CSI explanatory power with a bilinear ridge regression of the observed yield fluctuations at national level (FAOSTAT data), accounting for the covariance of the two agro-climatic indexes (see the Data and Methods section). Thus, the resulting multiplicative coefficients combining the HMD and the SPEI into the combined index depend on the country. The comparison between their absolute values can be interpreted as our best statistical estimate of the relative contribution of heat and water stresses on the national yield anomalies. Moreover, the sign of the parameter multiplying the SPEI indicates, for each country, larger yield sensitivity either to drought or to water excess.

We test the robustness of our proposed approach by performing a specific case study at the subnational level for France, where wheat is resulting to be more sensitive to over-wet conditions. This country is indeed characterized by pronounced heterogeneity of agro-climatic zones and by the availability of both high-quality yield data (from 1989 to 2014) and meteorological data (updated in near real-time). This allows us to discuss also the possible compensation effects while aggregating the indexes over large countries to compare with the FAOSTAT national yield statistics. Finally, we demonstrate the usefulness of the CSI to characterize the impact of climatic extremes on wheat yield anomaly observed during the 2016 wheat-cropping season in France.

Most physiological processes responsible for the final crop yield are negatively affected by maximum daily temperatures above 30 °C (Porter et al 2014). Above 34 °C–35 °C, wheat reaches its limits for survival (Porter and Gawith 1999, Hatfield et al 2011, Rezaei et al 2015), the so-called 'cardinal temperatures' for wheat (Porter and Gawith 1999).

To provide realistic estimates of fixed based threshold estimators often used to characterize heat stress on wheat, biases in temperatures need to be corrected, especially when gridded temperature data and/or reanalysis are used instead of local field measurements. Furthermore, wheat physiology is sensitive to the temperature of its canopy, which might be higher than the 2 m temperature (the commonly available variable in the observational data sets) in unstably stratified surface layers, especially above arid soils (Siebert et al 2014). The thermal vertical structure of the atmospheric surface layer can also be altered by evaporative cooling due to irrigation (Tack et al 2015, Troy et al 2015, Mueller et al 2016). These issues, together with the lack of information about the varieties and the detailed agronomical practices that are implemented in the different regions, raise some concerns about the uniform global applicability of heat stress indicators based on the cardinal temperature thresholds used in previous studies (Lizumi and Ramankutty 2016, Ray et al 2015).

Here, we adopt indicators based on statistically determined climate thresholds aiming to account for the physiological effects of heat and soil moisture anomalies (and, in principle, other correlated factors) with respect to the local climate variability.

The heat analysis is inspired by the HWMId (Russo et al 2015, Zampieri et al 2016). A heat wave is defined as a consecutive period when T_max exceeds the 90th percentile of T_max, computed for all days of the year considering a 2 w time window (T₉₀). Over these periods, the HMWId is then defined as the sum of the difference between T_max and the 25th percentile (T₂₅) normalized by the interquantile temperature range (i.e. the magnitude function M_d, equation (1)), computed in a similar manner to the 90th percentile.

$$\begin{equation} \begin{array}{*{20}c} {M_d (T)} \hfill & = \hfill & {(T - T_{25} )/(T_{75} - T_{25} )\;\;\text{if}} \hfill \\ {} \hfill & {} \hfill & {T > T_{90} ;\;\text{zero otherwise}} \hfill \\ \end{array} \end{equation} \tag{ 1 }$$

This non-parametric approach (i.e. without prior assumption on the statistical distribution of T_max) takes both the length and the intensity of the heat waves into account. Whether a severe or a 'close-to-normal' event is diagnosed depends on the numerical value that the HWMId assumes. A great advantage is that the severity of the events can be defined on the basis of global thresholds, no more depending on the region (Russo et al 2014, 2015, Zampieri et al 2016).

This definition of the HWMId is the same described in Russo et al (2015). Here, we slightly modify it to diagnose heat stress and temperature anomalies that are relevant for crops. While the original formulation of the HWMId is based on the maximum event over the whole year, we adapt it for wheat crop by cumulating the magnitude function over the last three months of the growing season before harvesting similarly to the way the Growing Degree Days (GDD) is defined. Therefore, the modified index accounts for all events occurring during the sensitive stage for the crop. We refer to this modified index as HMD. Figure S9 in the online supplementary information provides a visual representation of the HMD computation, with respect to the T₉₀ and the interquartile temperature range computed during 2003 growing season in central France.

Daily T_max data are retrieved from the AgMERRA climate forcing dataset (Ruane et al 2015) that was created to support the Agricultural Model Intercomparison and Improvement Project (AgMIP). The AgMERRA dataset is available at a 0.25° × 0.25° resolution from 1980 to 2010 globally. This dataset is based on the NASA Modern-Era Retrospective Analysis for Research and Applications (MERRA, Rienecker et al 2011), blending models and observational data through a mathematical procedure called 'data assimilation', but it corrects biases in temperature using CRU (Climate Research Unit) data. CRU data are amongst the most widely used climatic datasets. They consist of a gridded datasets derived from quality checked station measurements, excluding time-series that do not pass several quality and homogeneity tests (see Osborn and Jones 2014, and references herein). CRU data are defined on a monthly time-scale, while the day-to-day variability is provided from the MERRA reanalysis. We evaluated the HMD considering also AgCFSR, i.e. the other main agricultural forcings dataset adopted in AgMIP (which is also correcting the bias using CRU data, Ruane et al 2015), finding almost equal results (see figure S7 of the supplementary material, available at stacks.iop.org/ERL/12/064008/mmedia).

The soil moisture analysis is tackled with the SPEI (Vicente-Serrano et al 2010, 2013). The SPEI is a multi-temporal-scale water balance index integrating the precipitation and potential evapotraspiration budget on the basis of CRU data, available at 0.5° × 0.5° resolution from 1901 to 2014 globally (http://spei.csic.es/). It can be used for determining the onset, duration and magnitude of drought conditions with respect to normal conditions in a variety of natural and managed systems such as crops. SPEI anomalies are considered large, or extreme when specific thresholds (globally defined) are surpassed.

To compute the HMD and evaluate the SPEI over the wheat growing season and production areas, we use the global data set of monthly irrigated and rainfed crop areas around the year 2000 (MIRCA2000, Portmann et al 2010). MIRCA2000 integrates data from several different sources, with the specific aim to enhance the consistency with subnational statistics collected by countries and by the FAO. It provides both irrigated and rainfed crop areas of 26 crop classes for each month of the year, and includes information on the occurrence of multiple cropping systems allowing to account for winter and spring sown wheat, or irrigated and rainfed wheat when they coexist in the same grid cell. Portmann et al (2010) also provide an estimate of the uncertainties of MIRCA2000 in comparison with other sources. MIRCA2000 is available at various resolutions. Here, we use the 0.5 × 0.5 degrees resolution, which is compatible with the meteorological data used in this study. MIRCA2000 is consistent with the only other global agricultural calendar data available to our knowledge, as it is mostly based on the same sources (Sacks et al 2010).

MIRCA2000 (online supplementary figure S1) provides sowing dates but not the varieties. The majority of the harvested areas correspond to rainfed winter sown wheat (118 million hectares), especially in the middle latitudes (in both hemispheres). Rainfed spring sown wheat (30 million hectares) is dominant mainly in the northwestern US and Canada, and it is present in other regions (Spain, England, Poland and India) as secondary crop. Irrigated winter sown wheat (50 million hectares) is found mainly in India and China, and it is characterized by the largest local concentrations. Irrigated spring sown wheat (8 million hectares) harvested areas are almost negligible with respect to the other categories. Online supplementary figures S2 and S3 show the sowing and harvesting months, respectively, derived from the MIRCA2000 dataset.

Wheat is mostly sensitive to climatic extremes during the three months period before harvesting (Lobell and Tebaldi 2014). We define these periods for each grid cell of the climatic data using MIRCA2000 (see figure S3). Accordingly, we compute the indexes separately for each of the four wheat cropping categories (i.e. irrigated and non-irrigated crops, and double cropping), since the sowing and (most importantly) the harvesting dates can be different. In the grid boxes where different categories coexist, we average the values of the indexes using the harvested area of each wheat category as a weighting factor in order to obtain one value of each index for each year and each grid cell.

The HMD is computed for period 1980–2010 by summing equation (1) over all relevant events during the three months before harvesting. The SPEI, which integrates the precipitation and evaporation anomalies occurred earlier, is evaluated at the harvesting month. The timescale (i.e. the number of months) defining it has to be specified. Assuming that the crop is sensitive to soil moisture anomalies in the same period as for heat, but also accounting for the effects of the previous precipitation and evaporation anomalies in determining the current soil moisture values, we used SPEI6. In this way, we account for the water balance anomalies during the three months sensitive period before harvesting and also during an earlier preconditioning period of preceding three months, accordingly with the typical time scales of soil moisture dynamics (see e.g. Zampieri et al 2009).

Online supplementary figures S4 and S5 show the daily maximum temperature 90th percentile (T₉₀), used as a threshold in the heat stress estimation, corresponding to the beginning and to the end of the three months period from harvesting. Notably, there are regions where the T₉₀ is above 30 °C already at the beginning of the growing season, especially for irrigated winter wheat in India, and for rainfed spring wheat (figure S4). On the other hand, there are regions as Canada, for instance, where T₉₀ is relatively low (this is confirmed by Jarvis et al 2008). T₉₀ is above 30 °C in most regions at harvesting time (figure S5).

Global, national and subnational spatial aggregation are computed by averaging grid-point data, using again the MIRCA2000 harvested area data as weighting factors, so that the indexes can be compared with FAOSTAT3 national statistics of wheat production and yield (www.faostat3.fao.org/). The quality of yield data represents a main issue in this kind of studies, especially for countries where there are no alternative estimations available (Desiere et al 2016). The subnational analysis for France is based on winter wheat yields from 92 French départements, provided by AGRESTE Ministère de l'Agriculture (AGRESTE, 2015) and covering the period between 1989 and 2014. Subnational crop yields have been quality controlled (Ceglar et al 2016). Weather data for the computation of the CSI has been obtained from the high-resolution EOBS gridded weather dataset over Europe (Haylock et al 2008).

Since most of the tendencies affecting yield statistics are mostly due to improvement in the agricultural practices, we adopt a non-linear detrending procedure, i.e. the locally weighted scatterplot smoothing (LOESS, Cleveland and Devlin 1988) with the 'lambda' (span) parameter equal to 0.75. This value was optimized by trial-and-error procedure. However, the linear correlation coefficient we obtain are only weakly sensitive to changes of the span parameter around this value. In fact, the final results are robust to changes of the lambda parameter within a reasonable range of values. Moreover, this technique yields comparable results with respect to other analog procedures (see e.g. Ben-Ari et al 2016). This processing has been conducted with the R software build-in loess function.

Having removed in this way also the part of the signal that is potentially related to climate, we are compelled to apply the same procedure also to the HMD and SPEI time-series. This is performed only when the time-series exhibit either a significant trend (with the Mann–Kendall test), or a change point (Bai and Perron 1998), using 'Kendall' and 'strucchange' R libraries; otherwise, we just remove the mean. In such a way, we can isolate the effects of climate anomalies and extremes on the year-to-year wheat production variability. Online supplementary figure S8 shows an example of LOESS detrending (and mean removal) of yield and of the indexes time-series.

Our analysis is based on the anomalies obtained after such procedure. The CSI is then defined as a simple linear combination of the standardized HMD and SPEI (i.e. SPEI6) anomalies as expressed by equation (2):

$$\begin{equation} \text{CSI} = a \cdot \text{HMD}_{\text{detrend},\text{std}} + b \cdot \text{SPCI}_{\text{detrend},\text{std}} \end{equation} \tag{ 2 }$$

where the subscripts indicate the fact that the time-series are detrended and standardized. Standardization of the SPEI is advisable for using it in period 1980–2010, consistently with the reference period used to compute the HMD.

In equation (2), the a and b parameters are estimated by ridge linear regression with the country level FAOSTAT data (again after detrending, standardization and inversion [i.e. multiplication by −1]), accounting for the co-linearities between heat stress and drought. Therefore, the CSI is an estimate of the standardized yield anomalies. The numerical values (and sign) of the a and b parameters depend on the country. The sign of a is always resulting positive, as the CSI is defined such that positive values are corresponding to negative yield anomalies, which is the reason of the inversion of yield data. On the other hand, the sign of the b parameter indicates if the national yield is more affected by negative or positive soil moisture anomalies. Since all time-series are standardized, the comparison between the values of the parameters provide a qualitative indication of the relative importance of the two regressors in determining the yield anomalies. A quantitative indication could be provided if no co-linearities between the two regressors were present. In any case, the value of the a parameter provides an estimate of the effects of heat stress alone in determining the yield anomalies (i.e. by imposing no soil moisture anomaly in equation (2)). We refer to the a parameter as the heat sensitivity parameter, and to the b parameter as to the soil moisture sensitivity parameter.

Linear correlations between the individual regressors and national yield data are computed as well, for comparison with the CSI skill and to test the robustness and consistency of the procedure. In all three cases, Pearson linear correlations and tests (10%, two-sided) are applied in order to compute the statistical significance of the correlations coefficients.

Global wheat production approximately doubled from 1970 to 2015 (figure 1(a)), and increased by about 50% in the period 1980–2010 (the one here considered for the climatic analysis). The largest contributor to the global inter-annual variability is the Russian rainfed winter wheat production (see online supplementary figure S1), characterized by a standard deviation of the production anomalies of about 12 MTons corresponding to the 25% of the mean national production of 48 MTons per year (2000–2015 average).

Zoom In Zoom Out Reset image size

Wheat production interannual variabilities of EU (rainfed, mainly winter wheat), China (mainly rainfed and irrigated winter wheat) and USA (rainfed winter wheat) are following in decreasing order. Standard deviations consist of about 7 MTons that correspond to the 5%, 6% and 11% of the averaged national productions of 140, 107 and 57 MTons, respectively. Ukraine (rainfed winter wheat), Australia (rainfed winter wheat), Canada (rainfed spring wheat) and India (all agro-management types) are characterized by decreasing standard deviations (between 5 and 3 MTons) corresponding to the 27%, 19%, 14% and 4% of the averaged productions of 22, 24, 25 and 80 MTons, respectively.

These eight countries and regional aggregations together produced, on average, 495 MTons of wheat, corresponding to the 77% of the global amount of 640 MTons. The associated standard deviation almost equals the global value of 22 MTons, which is 3.4% of the total production (and the 4.4% of the sum of the 8 producers). Quantitative analysis of the fluctuations attributable to the individual regions displayed in figure 1(a) reveals that concurrent anomalies of two of these producers are able to explain more than 50% of the global inter-annual fluctuations in almost all years considered, if no compensating anomalies of the opposite sign are present.

Figure 1(a) also shows the occurrence of climatic extremes in the wheat cropping regions, as represented by the CSI. During the period 1980–2010, significant negative production departures from the mean global tendency are often captured by the CSI, especially in the years 1987–1988, 1994, 2003, 2006–2007 and 2010, coinciding with the observed spikes in cereal prices (see figures 7-3 in IPCC-WG2). We observed a significant level of consistency also in the years characterized by positive production anomalies and less frequent damaging climatic extreme events, like in 1986, 1992–1993, 1997, 2004–2005 and 2008–2009. The abundant wheat production in 1990, being characterized by moderate climate conditions, represents an outsider in our analysis. That was the last year of intensive fertilization in Russia (Lioubimtseva et al 2015), which is the main contribution of the recorded production anomaly.

Overall, most of the largest anomalies that could be related to climate variability are induced by the major producers, while smaller and/or compensating positive and negative anomalies characterize the 'normal' years. This is partly expected because the single CSI anomalies are weighted by the national productions in figure 1(a). However, the ranking of the anomalies contributions by region slightly differs from the one of production: the largest contributor to the variance of the CSI is Europe, followed by the US, China, Russia and India.

This difference could be a consequence of the spatially variable performance of the CSI in capturing the stress events for wheat, or to the different relative importance of climate compared to other factors in determining the production anomalies in the different regions, or it could be an effect of the spatial aggregations in large heterogeneous regions. Higher resolution yield and production data at the global scale, together with detailed agronomical information, would very likely improve the general understanding and ease the modeling of the yield anomalies.

The overall performance of the CSI is quantified by the linear correlation coefficients with yield anomalies, plotted in figure 1(b). We estimated significant linear correlation coefficients for most countries. Among the largest producers, US, China, India and Argentina display smaller but statistically significant values. As for the US, this result is likely to be caused by the mitigation effect of large-scale irrigation that tends to decouple yield anomalies from climate variability also in regions where it is not applied directly (Tack et al 2015, Troy et al 2015, Mueller et al 2016).

The linear correlation between the total production and the global averaged CSI computed for all wheat cropping areas is −0.65, which corresponds to the 42% of the total variability, strengthening the result of Ray et al (2015).

The effects of heat stress are comparable or larger than the effects of water stress in most countries (see figure 1(c)). Water stress effects tend to be more substantial in the Mediterranean region. Negative WSI can be associated to drought and related to negative yield anomalies, especially in arid and semi-arid regions. In the Tropics and the middle/high latitudes, our statistical model suggests an opposite relation with positive WSI being associated to negative yield anomalies (see figure 1, bottom-right panel). This is partly consistent with the findings of Vicente-Serrano et al (2013) who showed a reduced sensitivity of natural vegetation dynamics to drought in regions characterized by climatological water excess. Our study suggests that wheat may be, on average, more sensitive to water excess extreme events than to drought in these regions, and that dryer-than-normal conditions could have, on the other hand, beneficial effects.

Figure 2 shows the spatial and temporal characteristics of the individual regressors forming the CSI, as well as their direct relationships with the FAOSTAT yield data.

Zoom In Zoom Out Reset image size

At the global level, heat stress alone (HMD, figure 2(a)) has an explanatory power similar to the CSI, and displays a larger trend. Linear correlation between total wheat production and globally averaged HMD anomalies is −0.57, slightly lower than the value of −0.65 obtained for the CSI. HMD has the advantage, w.r.t. to the CSI, that it does not need calibration with yields data. However, including the SPEI integrating the water balance in the six months prior to harvesting (i.e. SPEI6, see Data and Methods) is necessary to significantly explain the production variability in India, China and France, and to improve the correlation coefficients in several other countries (compare figure 1(b) and figure 2(c)).

Given the dual role of the SPEI in the CSI (figure 1(d)), the globally averaged SPEI (figure 2(b)) is not significantly correlated with the global wheat production anomalies. With respect to the HMD, SPEI6 time-series are less often significantly correlated with the national productions (compare figure 2(c) and figure 2(d)). In most of the countries the sign of the correlation is positive, i.e. negative cumulated water balance is associated with negative yield anomalies. This is correctly accounted for in the combined index. The CSI is always outperforming the individual indexes such that it can be significantly correlated with the national yield anomalies even if the HMD and/or SPEI are not.

It is worth noting that heat and water anomalies are significantly linearly correlated in several countries (figure 2(e)), but there are many regions where heat and soil moisture anomalies and extreme events seem to be more independent. While several authors suggested a tight linkage between these two types of extreme events in water limited regions (see e.g. Zampieri et al 2009, Seneviratne et al 2010), our analysis suggests that this is not often the case for most of the wheat producing regions that are found to be more sensitive to water excess than to drought (figure 1(d)).

France, with about 37 MTons of rainfed winter wheat produced per year, is the largest producer of the European Union (25%), and consequently it is responsible for a corresponding fraction of the inter-annual standard deviation. Therefore, EU and France production variation are significantly correlated, with a linear correlation coefficient of 0.49. However, European climate is characterized by pronounced spatial gradients and prominent production anomalies of the individual countries (see figures 1(b)–(d)) that can compensate at the European level. This compensation effect has been observed to moderate relatively large anomalies preferably in the first part of the time series, in particular in 1980, 1985, 1986, 1988. It appears to be reduced towards the end of the period. Consistently, the linear correlation of the French and EU productions computed until the end of the 1980s is reduced to 0.18. This implies that there is a current tendency towards more coherent anomalies, and, therefore, towards an increased threat to food security. This shifting behavior might also be manifested by other large countries and spatial aggregations. Indeed, it is consistent with the emerging global scale heat waves pattern observed recently (Zampieri et al 2016).

According to our analysis at the national scale, the linkage between climate and yield variability in France is controversial. Although there is a significant correlation between heat waves and drought (see figure 2(e)), the heat and water stress indexes are not directly correlated with the production anomalies and the French production is more negatively affected by water excess than drought. This motivated us to conduct a deeper analysis for France.

Figure 3 represents the main results of the analysis conducted for France at the subnational level. Yield data are available from 1989 to 2014, while the meteorological data are updated in near real-time, allowing the evaluation of the indexes for more recent years.

Zoom In Zoom Out Reset image size

The CSI is able to explain yield variations significantly in a considerable part of the country (figure 3(a)), with a linear correlation coefficient of about −0.5, especially in the central-northern part where the main wheat producing areas are located. In western France, facing the Atlantic, no significant correlations are found. The analysis conducted for heat (HMD) and water (SPEI) indexes alone (figures 3(b) and (c)) shows interesting spatial variability. The central part of France is significantly affected by heat waves. Mediterranean France is more sensitive to drought while the northern part results to be more sensitive to water excess (consistently with Ceglar et al 2016). The CSI computed at the subnational level indicates higher sensitivity to water excess than to drought (i.e. negative 'b' parameters, see Data and Methods) in the northern and central parts of the country, where most of the wheat is grown. This is consistent with the results obtained at country level.

We have used the CSI as a diagnostic tool to assess the impact of weather events on French wheat yield in 2016, characterized by an extremely negative anomaly. The indexes computed for 2016 (figure 3(e) and (f)) are diagnosing an excess of water in central, northern and eastern parts of France (mainly due to heavy precipitation in May), while drought occurred in the southern part of the country. In addition, a moderate heat wave influenced the western and southwestern regions (towards the end of the growing season). Therefore, the combined index computed for 2016 would have indeed suggested a negative yield anomaly (figure 3(d)).

By using the CSI to estimate the production anomaly at country level in 2016, we obtain a reduction of about 1 MTon lower than the average, corresponding to the 40% of the national production standard deviation. Therefore, the CSI is probably underestimating the actual anomaly. However, it suggests that the adverse weather conditions responsible for the water excess (possibly around the wheat flowering stage) might have played a central role.

It is worth noting, however, that the CSI does not capture the single precipitation events, nor their timing, but it rather accounts for their cumulative effect on the soil water balance. Therefore, while the CSI would have predicted a negative yield anomaly, the precise causes of the 2016 French yield loss still have to be determined.

We have analyzed the historical global wheat production annual anomalies computed from the FAOSTAT data set from 1980 to 2010. During this period, few of the major producing regions were responsible for significant anomalies at the global level. Concurrent anomalies of two of the major producers suffice to explain more than 50% of the global inter-annual fluctuations in almost all years, if no significant compensating anomalies were present. Russia is the largest contributor to the global annual variability.

The combined stress index (CSI), defined as a linear superposition of the Standardized Precipitation Evaporation Index (SPEI, Vicente-Serrano et al 2010) and of the Heat Magnitude Day (HMD, this study), allows explaining the 42% of the variability globally (confirming and strenghtening the result of Ray et al 2015), and often at country level as well. Therefore, accurate yield forecasts in the major producing countries could provide useful insights in global production, potentially increasing food security. In fact, we are currently working on the application of the CSI in a seasonal forecasts context (Ceglar et al submitted).

Our result points to a clear role of climate anomalies and extreme events on wheat yield anomalies and enables identifying the relevant contributors and the associated effects. Heat stress is often the most important predictor, consistently with field studies and future projections, and it is in general as important as drought. As a prominent exception, we have found that in the Mediterranean countries drought carries a larger detrimental effect on wheat yield than heat stress.

Heat stress over wheat cropping regions increased significantly in the period 1980–2010, especially since the mid-1990s. This produced less compensating and more concurrent yield anomalies, motivating the general concern about food security.

Our results also point to the sensitivity of wheat yield to the water excess issue, rather than to drought, especially in tropical regions and in some regions of the mid/high latitudes. While heat stress globally remains the most important factor determining the yield anomalies (compare figure 2(c) and figure 2(d)), our findings show that water stress, and in particular water excess, is essential to explain yield anomalies of important wheat producers such as China and India.

The results obtained at the subnational level for France (the main European wheat producer) support and strengthen our main findings. The combined indicator was, indeed, able to capture sub-regional yield anomalies, and provided consistent results with the national scale analysis. This case study also showed the importance of accounting for water excess to gain some predictive skill for extreme events such as the large wheat production reduction occurred in France of 2016.

We are indebted with Maurits van den Berg for his help during the revision process and two anonymous Reviewers that greatly contributed to the improvement of the quality and clearness of the message we are conveying.