Climate variation explains a third of global crop yield variability

Many studies have examined the role of mean climate change in agriculture, but an understanding of the influence of inter-annual climate variations on crop yields in different regions remains elusive. We use detailed crop statistics time series for ~13,500 political units to examine how recent climate variability led to variations in maize, rice, wheat and soybean crop yields worldwide. While some areas show no significant influence of climate variability, in substantial areas of the global breadbaskets, >60% of the yield variability can be explained by climate variability. Globally, climate variability accounts for roughly a third (~32–39%) of the observed yield variability. Our study uniquely illustrates spatial patterns in the relationship between climate variability and crop yield variability, highlighting where variations in temperature, precipitation or their interaction explain yield variability. We discuss key drivers for the observed variations to target further research and policy interventions geared towards buffering future crop production from climate variability.

. Model bias. The bias is scaled and thus the legend reflects a bias of -5% to +5%. White areas indicate where crop is not harvested or analysed.
Supplementary Figure 9. Model significance levels as determined from conducting F-tests against a random climate model. Regions in green had p=0.01 (sample size study period of 3 decades per political unit; F-test) or models at 99% significance; yellow areas had significance between p=0.01 and p=0.05 or 99% to 95% significance (study period of 3 decades per political unit; F-test). Most of the regions had models at p=0.05 level. We additionally had included some more areas with p values between p=0.05 and p=0.1 or 90% to 95% confidence. White areas indicate where crop is not harvested or analysed.  [12] found similar results using either the CRU or the UDEL global climate dataset; hence we felt that using only the CRU TS3.10 was sufficient to conduct this study.
The spatial resolution of our gridded crop data (see next section) is at 5 minutes globally. Thus, within each half degree grid cell of the CRU TS3.10 data are several crop grid cells. We assigned each of the 5 minute crop grid cell that was contained within a half degree CRU grid cell the identical reported weather value. Those that lay partly within more than a single CRU grid cell were assigned the value of the nearest CRU grid cell. However note that the analysis was done at the political unit level which can contain numerous CRU grid cells. What we computed was the harvested-area-weighted temperature and precipitation information for all of the grid cells within a political unit each year; this was the political unit's climate value per month for any specific year.

Crop data
The crop yield and area harvested information came from the newly developed dataset noted in three of our previous reports [13][14][15]. In view of the importance of the Russian Federation especially for wheat and maize production we updated Russian maize and wheat statistics at the oblast / county level using an approach similar to those described in Ray et al. [2012] [14] specifically for this study; otherwise the rest of the crop data used is similar to the ones reported in [14,15].
The crop statistics database provides information on crop yield and harvested areas for each of the four crops: maize, rice, wheat, and soybean at the 5 minute spatial resolution. The yield and harvested information does not distinguish between the primary crop and secondary crops such as rice or soybean grown in double cropping systems. Similarly, crops such as wheat that are grown as winter or spring wheat are not distinguished. There are several other agronomic distinctions in cropping systems, such as irrigated crops, different cultivars, upland and lowland crops that would likely have different levels of sensitivity to climate variability. We did not make these distinctions in our crop database. The original dataset covers the period of 1961 to 2008 and contains ~2.5 million unique statistics of maize, rice, wheat, and soybean harvested and yield information. We used the period 1979 to 2008 for this study which contained ~2 million reported statistics on these four crops including the updated information for the Russian Federation [16-18].

Planting and harvesting dates
We used the main growing season determined by the average planting and harvest date for each crop as given in [19]. The crop planting and harvest dates have varied in the last three decades due to climate as well as technological and socioeconomic changes at each political unit [20]. In a perfect modeling system we would have the planting and harvesting dates as a function of time. At this point we have not built a dataset reflecting temporal crop planting and harvest dates. Moreover, the climate information is monthly, reducing sensitivity of the analysis to planting and harvesting dates. However, we did test the sensitivity of results at the top maize producing U.S. counties using daily station data and found statistically significant correlation (at p = 0.0) between monthly and daily climate data based analysis. In most global locations a few days advancement or delay in crop planting and harvesting still leads to the same month. Thus we used fixed crop planting and harvested dates as used in other studies [12,21].
We however updated the global maps of crop planting and harvest dates reported in [16] based on more recent information for some of the important agricultural countries such as China [22], India [23], Pakistan [24], Bangladesh [25], Nepal [26], Australia [27], South Africa [27], Brazil [17,28] and Argentina [17,29] (Supplementary Figure 6). We updated the planting and harvest dates for only maize and wheat. Rice crop appeared generally more accurate globally but note that in several global locations there is a significant second rice growing season. Our utilization of weather conditions 1 year prior to harvest of the main rice crops provides information regarding rice growing conditions for the secondary season. In future studies a more precise delineation of the main and the secondary growing season yields and weather conditions needs to be modeled.

Setting up of the statistical models
We set up the statistical models in three steps. We first computed the harvestedarea-weighted average seasonal and average annual (12 months prior to harvest month including the harvest month) temperature and precipitation information at each of the political units that harvested a crop (maize, rice, wheat, and soybean). In the next second step we detrended the climate and yield information by taking out the linear trends for climate and the best fit yield model following [12 and 14] for yields: where Y, T, and P are the observed yield, temperature and precipitation, and the subscript 'd' denotes the detrended values, 'l' denotes the linear fit, and 'b' denotes the best fit from among constant, linear, quadratic, and cubic fit. The best fit in turn was determined using AIC [30]. The * denotes the initial value similar to detrending noted in [12].
We used the crop yield trend modeling described in Ray et al., 2012 [i.e. reference 14] but other investigations into crop yield modeling are also present. Specifically in [14] comparisons with the statistical method presented in another study [e.g. 30], which in turn were further investigations from other recently published methods [31,32] was shown. The yield modeling method presented in [14] was applied to national level data to compare against the national scale yield modeling presented in [30] and gave similar results. Yield modeling however continues to be an area of active research [33,34], but their performance at our level of high spatial resolution globally remains unknown, as other studies have been either regional studies or national scale studies. Future studies could investigate these developing techniques of yield modeling in global analysis of yield variability.
Only those political units that harvested a crop (maize, rice, wheat, and soybean) at least half of the years of the study period entered into our analysis. This means that a political unit needed to harvest maize crop more than 15 out of the maximum 30 years to enter into our maize analysis and so on.
In the third step we built linear regressions connecting with the single and squared forms of and . Note, that there are both seasonal and annual forms of the climate variables. Thus, eight types of climate variables are possible and 27 different forms of: were generated. We had a sufficiently large number of permutations to cover a large set of statistical relations between the observed crop yields variability and climate variability.
Since we tracked ~13,500 political units this meant ~1.5 million statistical equations were generated for the four crops. We identified the best fit equation using the Akaike Information Criterion (AIC) [35] from the 27 equations at each political unit which resulted in using the best ~50,000 equations to draw on our results and conclusions for the four crops. Supplementary Figure 7 gives an example of this process from Shaanxi province, China. Supplementary Figure 1 provides the average distribution of crop yields.

Determining yield variability explained by climate variability
From the 27 regression equations the equation with the lowest AIC was further analyzed. Even though the model chosen fit the data best it may still be a poor model and no better than a model that related the variability in crop yields to random climate variations or a null model. Thus we conducted an ANOVA or F-test at the p = 0.10 level to determine whether the chosen model was statistically significantly better than a random climate model. Only in those political units where the chosen model was statistically significant at the p = 0.10 or less level we further reported how much of the yield variability was explained by climate variability. Political units where the p value was greater than 0.10 indicated that there was no statistically significant relationship in crop yield variability and climate variability. We classified these political units as those where we could not detect statistically significant climate variability impacts on crop yields as "no effect" in Figure 2 in the main text.
The coefficient of determination, or explained variance for the chosen model, i.e. the one with the lowest AIC (equation 4), is statistically the 'complete model' explained variance or R 2 which is plotted as categorical maps in Figure 2 in the main text. Reduced models of temperature and precipitation were also generated with only temperature and precipitation terms. If the selected model had only temperature terms then only the reduced model for temperature was generated and vice versa. The reduced models give an indication of how important temperature and precipitation variability was in explaining yield variability and the R 2 of these models are given in Supplementary Figures 2 and 3.

Global maps
Even though the global maps are plotted at the grid cell the analysis was conducted only per political unit and the result for the political unit mapped onto all the grid cells. The bias of the selected model at each political unit, and for each crop, is plotted in Supplementary Figure 8.
One could also choose Bayesian Information Criterion (BIC) as the criterion to identify the best equation from the 27 regression equations [36], but we continued to use AIC similar to our previous studies [14,15]. Political units where crop yields completely failed would not report any harvested areas and production. These areas evidently would be undergoing an extreme event and these situations of total crop failure would not enter into the current analysis as harvested areas reported would then be zero. The numbers reported in the summary data table for each country and crop are weighted by harvested area. Further Supplementary Figure 9 provides the spatial locations of p values of the selected models for two increasingly higher confidence limits (p=0.05 and p=0.01). The models selected at the less restricted pvalue of 0.1 generally included models at p=0.05 or less.

Sensitivity to time period analyzed
Our study period spanned the thirty years from 1979 to 2008. We also conducted separate analyses to determine how sensitive the explained crop yield variability was to the length of time studied. As expected, results were slightly different.
Simulations for a twenty-five year period from 1984 to 2008 resulted in climate variability explaining 40%-45% of the overall yield variability. Over a thirty-five year period (1974 to 2008) the numbers were 29%-34%. Our second set of sensitivity analysis involved conducting thirty year analyses but starting the simulations from one and two years earlier to 1979, i.e. 1978 and 1977. In both cases the results were nearly identical to numbers for analysis starting in 1979. When simulations started off from 1978, we found that a nearly similar 32%-39% of crop yield variability was explained by climate variability and for the analysis from 1977 to 2006 the numbers were 34%-39%.

Use of gridded weather versus station data
Since the yield is measured at each political unit, the weather information should appropriately represent all crop-harvested areas within the political unit and valid for each political unit similar to the procedure used in [12]. Gridded weather data such as the CRU [1] use station data, which have sparse coverage (see for example the coverage in the top maize counties of the United States -Supplementary Figure  4) to build a complete coverage through statistical analyses. Direct use of station data, which essentially represents only a short fetch around the site of the ground station ( Supplementary Figure 4), should not be used to analyze the relationship between crop yields when measured at the political unit (covering a vast area) and weather. On the other hand if the yield is measured at a specific location, for example a farmer's field, then gridded weather data should not be directly used without appropriate downscaling of the data on account of relatively coarser gridded data at this scale. A better approach would then be to use a nearby sited station weather information. Our analysis being at a political unit, the unit of yield measurement, we upscaled the gridded weather data (as described in section 1.1, above) to match the yield measurement spatial unit, a method similarly used in [12].
We however conducted an alternate study to show how different our results would have been if spatial mismatch using direct station data was used. Since we did not build the CRU data, the analysis also serves as a rough check of the CRU data.
We selected the top ~100 maize producing counties in the United States (representing ~25% of total United States maize production) to perform this additional experiment. Over the study period we accessed all the archived weather station data for all the ~100 counties [37] and used the average of all, essentially point measurements, within a county boundary to generate a seasonal and annual weather for the county. The stations in general had very poor coverage of the county maize harvested areas (Supplementary Figure 4). Note that the CRU gridded data is essentially station data with statistically estimates of the weather in non-covered areas resulting in complete coverage over all land areas [1].
Next we conducted our statistical modeling directly using the station data and determined the year-to-year yield variability explained by year-to-year weather station derived climate variability (using the coefficient of determination R 2 metric). Finally we correlated the station analysis R 2 with the CRU weather derived R 2 for all counties considered. The scatter plot (Supplementary Figure 5) shows both these R 2 per county. The correlation was significant (0.54 at p<0.001).