Uncovering the Past and Future Climate Drivers of Wheat Yield Shocks in Europe With Machine Learning

Recently, yield shocks due to extreme weather events and their consequences for food security have become a major concern. Although long yield time series are available in Europe, few studies have been conducted to analyze them in order to investigate the impact of adverse climate events on yield shocks under current and future climate conditions. Here we designated the lowest 10th percentile of the relative yield anomaly as yield shock and analyzed subnational wheat yield shocks across Europe during the last four decades. We applied a data‐driven attribution framework to quantify primary climate drivers of wheat yield shock probability based on machine learning and game theory, and used this framework to infer the most critical climate variables that will contribute to yield shocks in the future, under two climate change scenarios. During the period 1980–2018, our attribution analysis showed that 32% of the observed wheat yield shocks were primarily driven by water limitation, making it the leading climate driver. Projection to future climate scenarios RCP4.5 and RCP8.5 suggested an increased risk of yield shock and a paradigm shift from water limitation dominated yield shock to extreme warming induced shocks over 2070–2099: 46% and 54% of areas were primarily driven by extreme warming under RCP4.5 and RCP8.5, respectively. A similar analysis conducted on yields simulated by an ensemble of crop models showed that models can capture the negative impact of low water supply but missed the impact of excess water. These discrepancies between observed and simulated yield data call for improvement in crop models.

This study focused on yield shock, which is often accompanied by food price spikes and can have dramatic influence on food security and social conflicts (Chen & Villoria, 2019;Crost et al., 2018). Yield shock is normally interrelated with single or compound extreme weather events, pest and disease outbreaks or geopolitical events (Beillouin et al., 2020;Ben-Ari et al., 2018;Ciais et al., 2005;Cottrell et al., 2017;Webber et al., 2020). Yet, the influential climate variables are diverse and not fully identified, and interactions between them matter as well. For example, heat stress accelerates crop leaf senescence (Lobell et al., 2012); excessive precipitation causes lodging or increased risk of pests and diseases (Zampieri et al., 2017). Several studies investigated the impact of specific climate variation or extreme climate events on crop yields (Beillouin et al., 2020;Jin et al., 2017;Lesk et al., 2016;Rötter et al., 2018) but no study has been conducted to attribute yield shocks to different types of weather events in Europe over several decades. This information is essential to identify the most influential climatic factors and the most sensitive periods of crop development. Considering that compound extreme events are more frequent in warmer climate (AghaKouchak et al., 2020;Zscheischler et al., 2018) a better understanding of the leading drivers of yield shock will help design a targeted climate adaptation policy for agricultural production systems.
Winter wheat is widely planted in Europe and highly sensitive to climatic warming (Asseng et al., 2015;Porter & Semenov, 2005). Compared to spring crops, winter wheat has a longer growing season and thus higher probability of being exposed to extreme events, making it more difficult to be simulated by crop models (Guarin et al., 2020). For example, until approaching the time of harvest, all forecasting systems failed to predict the 2016 yield shock that occurred in the breadbasket of France. In contrast, a statistical model was able to attribute this extreme yield shock to a compound extreme of abnormal autumn warming and spring wet conditions (Ben-Ari et al., 2018). Therefore, a systematic analysis of the climate drivers underlying different types of shocks beyond well-studied drought events and an evaluation of the performances of current crop models to capture the drivers of those shocks is required to improve model credibility and help agricultural decision making (Rötter et al., 2011;Wang et al., 2017).
Here we built a data-driven framework with long term sub-national European winter wheat yield survey during 1980-2018 to attribute the yield shocks to different climatic drivers. We first used the random forest (RF) machine-learning algorithm to associate climate stresses at different wheat growth stages with wheat yield shocks. We then used an attribution analysis based on game theory (Shapley additive explanations) to assess the relative contribution of extreme weather to each yield shock event ( Figure 1). Specifically, in the following content, we try to address the following questions: (1) How have the wheat yield shocks evolved over the historic period? (2) What were the primary climate drivers for the identified yield shocks? (3) How will the yield shock occurrence probability change under future climate scenarios? (4) How well are the current crop models able to predict wheat yield shocks under current and future climate? ZHU ET AL.

Wheat and Climate Data
Winter wheat yield and cropping areas data in Europe over 1980-2018 were obtained based on survey data collected in a previous study (Beillouin et al., 2020) which covered 17 European countries at subnational scale (i.e., at nomenclature of Territorial Units for Statistics from EUROSTAT 2 and 3) and 1,435 administrative units across the 17 countries. To fill the missing yield data in some units, we employed a gridded yield data set at a spatial resolution of 0.5° (Iizumi and Sakai 2020). This global-scale gridded yield data was created through disaggregating national level yield data to finer spatial resolution based on grid cell crop harvest area and satellite estimated crop productivity. Correlation analysis suggests a good consistency between the two yield data sets ( Figure S1) in most units. Given the spatial disparity in climate and potentially locally adapted farmer management practices in wheat production system, the whole Europe yield data set was divided into four regions: North, South, West, and East Europe, as done by Beillouin et al. (2020) (details in Table S1).
Yield shocks were identified by fitting long-term yield trends through locally weighted scatterplot smoothing ("loess," Cleveland, 1979) and then estimating relative yield anomalies  , i t as: i t Y is the expected yield value in the i administrative unit at year t. The expected yield value, , i t Y corresponds to the long-term yield trend estimated using loess. Finally, yield shocks were defined as the lowest 10th percentile of relative yield anomalies in each region ( Figure S2). For those administrative units having missing values in the yield survey, relative yield anomalies based on yield survey and gridded yield by Iizumi and Sakai (2020) were estimated separately. Then yield shock was identified based on the two combined yield anomalies. The yield shock events identified in Iizumi and Sakai data constitute 18% total yield shock events. We also calculated the area fraction of yield shock for each year across all administrative units as: where Fraction t is the area fraction of yield shock for year t; , Gridded climate data during 1979-2018 were obtained from the Joint Research Center's MARS meteorological database (Toreti, 2014). This database was selected to be consistent with the climate forcing data used in the following crop model simulations. This climate data included daily minimum (T min ) and maximum (T max ) surface air temperature, precipitation, 10-m wind speed, global radiation, vapor pressure deficit and potential evapotranspiration from a crop canopy at 25 km resolution covering Europe and neighboring countries.
The winter wheat growing season in Europe generally starts with sowing in autumn followed by a dormancy period in winter and plant development in spring, and ends with maturity in the following late spring or early summer. Considering the different response of the plant to climate stresses in different growing stages (Duncan et al., 2015;Siebers et al., 2017;Zhu et al.,  stress, cold stress, water supply, and water demand. Specifically, we used two variables (mean precipitation and fraction of rainy days) to characterize the water availability accounting for precipitation distribution, as a recent study suggested that uneven precipitation distribution will offset the yield benefit of precipitation (Fishman, 2016). The 20 climate variables (five types of variable for each of the four growth stages considered) were used as the predictors of yield shock probability.

Crop Models and Future Crop Model Projections
Crop models have been used to study the climate change impacts on crop productivity, given their potential to simulate the large-scale spatiotemporal variations of crop yield (Mueller et al., 2017;Webber et al., 2018).
Here we used simulated yield data produced by eight-crop models (4M, FASSET, HERMES, MONICA, SIM-PLACE-Lintul5, SIRIUS 2015, SiriusQuality v3, and SSM-Wheat) which simulate wheat yield at 25 km resolution across Europe in the historical period 1980-2010 and future climate scenarios (RCP4.5 and RCP8.5). From the multiple model experiments, we selected the model simulation with responsiveness to mean temperatures, drought, and heat stress with simulated canopy temperature. Although other global-scale crop model simulation experiments were available (Mueller et al., 2017), the eight-crop models selected here were designed specifically for Europe and assimilated observed crop phenology data from Eurostat (http:// ec.europa.eu/eurostat/web/main) to better match crop growth seasonality. Additionally, the simulated yield showed a good performance in reproducing the interannual variation of wheat yield at national level.
For the historic period simulation , all eight crop models were forced by the climate data obtained from the Joint Research Center's meteorological database (Toreti, 2014). For the climate projections under the two future scenarios (RCP4.5 and RCP8.5), 6 out of 8 crop models reported the outputs forced with elevated CO 2 and future climate models (GFDL-CM3, GISS-E2-R, HadGEM2-ES, MIROC5, and MPI-ESM-MR) during the period 2070-2099. The climate projections were downscaled to the resolution of the baseline data by assigning values to each 25 km grid. The climate projections were created using the enhanced delta change method by applying changes in simulated temperature and precipitation variability to historic mean climate (Ruane et al., 2015). This climate data can be accessed at: http://open-research-datazalf.ext.zalf.de/ResearchData/DK_59.html. The same climate variables in Table 1 were derived for each climate model under RCP4.5 and RCP8.5, and were also used as predictors to predict future yield shock probability. In future climate scenarios, PET was calculated using Penman-Monteith equation (Allen et al., 1998).
To identify yield shock in the crop model simulation during both the historic period  and future scenarios (RCP4.5 and RCP8.5), simulated yield data were aggregated to the administrative unit and grouped into the same four regions as those considered in the yield survey data. Similarly, the relative yield anomaly was estimated with Equation 1 and then yield shocks (Shock sim ) were defined as the lowest 10th percentile of relative yield anomalies in each of the four regions.

Data-Driven Framework of Yield Shock Occurrence Modeling and Attribution
Yield shock events identified in the historical yield survey data were treated as a binary variable. We used a RF machine-learning algorithm to predict the shock occurrence probability using the 20 climate variables defined above (Table 1) as predictors. RF algorithm was selected as it was able to take into account non-linear responses and complex interactions between multiple predictor variables. The output of RF was expressed as a predicted probability of shock occurrence computed by averaging across an ensemble of ZHU ET AL.

10.1029/2020EF001815
4 of 13 f wd Fraction of warm days (days with T max >30°C or its 90th percentile in each growth stage) f cd Fraction of cold days (days with T min < its 10th percentile in each growth stage)

Prec
Mean daily precipitation in each stage (mm/day) f rd Fraction of rainy day (rainy day is a day with precipitation > 1 mm) PET Mean daily potential evapotranspiration amount in each stage (mm/day) decision trees built with a training data set (Breiman, 2001;Cutler et al., 2012). One separate RF model was built for each of the four regions to account for the potential disparity in the response of yield to climate variation. The two hyperparameters-"mtry" and "ntree"-in the RF model, representing the numbers of predictors sampled for splitting at each node and number of trees grown in the RF model, were tuned using the available set of data. The "ntree" was set equal to 1,000, since we found that out-of-sample model accuracy plateaued after this threshold. The "mtry" was tuned through 5-fold cross-validation and set equal to 11, 14, 12, and 11 for North, West, East, and South Europe, respectively. In order to assess the RF model performance in discriminating yield shock and non-yield shock, we computed the area under receiver operating characteristic curve (AUC) value (Fawcett, 2006) for each of the four regions by five-fold cross validation (AUC equal to 1 means a perfect classification, while a value of 0.5 suggests a performance equal to a random classification). For each RF model, we used partial dependence plots (PDP) to analyze the relationship between the yield shock occurrence probability and each predictor variable, averaging over all other predictors (Friedman, 2001). Similarly, eight RF models were built to fit the yield shock derived from the eight sets of crop model simulations with the same climate predictors during 1980-2010, respectively. RF models were implemented using the "ranger" R package (Wright et al., 2017).
Although PDP provides a graphical description of how yield shock probability responds to climate stressors in their multi-dimensional space, this approach does not quantify how much each climate predictor has contributed to the predicted yield shock probability. An attribution analysis based on game theory (Shapley additive explanations, SHAP) allowed us to evaluate the relative contribution of each climate variable to the predicted yield shock probability and then to identify the most adverse types of stress. SHAP is a method of estimating the expected marginal contribution of a covariate across all possible combinations and has been used in recent study to interpret machine learning model (Padarian et al., 2020). For each observed yield shock event, we used SHAP to break down the RF model predicted yield shock probability into individual predictor contribution. We then selected the predictor with the largest SHAP value as the primary driver for a yield shock event. Based on the category of climate variables in Table 1 and the reported prevalent causes of wheat yield shock events (Webber et al., 2020;Zheng et al., 2012;Zampieri et al., 2017), we focused on the following five drivers: extreme warming, high water demand, excessive water supply, low water supply, and cold stress. The selection criteria for each primary driver is detailed in Table 2. SHAP value based attribution analysis was implemented using the "iml" R package (Molnar et al., 2018).
Based on the criteria listed in Table 2, all yield shocks associated with the same primary climate driver were grouped together, and the five resulting groups were used to calculate the area fraction of yield shocks with the same primary driver for the four regions and the whole of Europe. We also divided the area fraction of yield shocks with the same primary driver in 1980-2018 into two periods (1980-1999 and 2000-2018) to investigate how a certain climate stress driven yield shocks have changed over time.
To predict the probability of yield shock in the future, we applied the RF model trained over historic yield survey data to the two climate scenarios considered (RCP4.5 and RCP8.5). Then, shock probability above its top 90th percentile in each region was identified as yield shock. We then applied the SHAP method to decompose the yield shock probability into contributions of individual climate predictors and finally determined the primary driver for each shock in future climate scenarios. Future yield shocks predicted by RF models (Shock RF ) were also used as pseudo true values to evaluate the performance of crop models in ZHU ET AL.  forecasting yield shock. A yield shock occurrence predicted by the RF model during both historic and future periods might be identified as yield shock by 0, 1, 2 ,…, 6 crop models (Shock sim ). The level of agreement between Shock RF and Shock sim was assessed as a function of the type of climate variables considered as the main driver by RF. For example, for N yield shocks (Shock RF ) primarily driven by extreme warming, we could get N0, N1 ,…, N5, N6 (N0 + N1 + … N5 + N6 = N, with N the total number of Shock RF for this driver) yield shocks also identified as shock (Shock sim ) by 0, 1, 2 ,…, 6 crop models, respectively. A high value of N4 + N5 + N6 + N7 suggested a strong agreement between RF and crop models.

Spatial and Temporal Distributions of Yield Shocks in Europe
Time series of wheat cropping area influenced by yield shock during 1980-2018 suggested that the year showing the highest area impacted by yield shock corresponded to year 2003 (Figure 2a). According to several previous studies, this year was characterized by a compound heat and drought stress (Ciais et al., 2005;van der Velde et al., 2010). When looking into the four regions, North Europe and East Europe show increasing trend in the area fractions of yield shock ( Figure S3). Spatially, the yield shock frequency (yield shock occurrence numbers/total years) is higher in southern areas characterized by warmer conditions, especially in Romania, Spain, and Southern France (Figure 2b).

Response of Yield Shock to Climate Stresses
AUC values estimated with 5-fold cross validation for the four regions were all greater than 0.85, suggesting that the trained RF model had a high capacity to discriminate between yield shock and non-yield shock, with North Europe scoring the highest AUC (0.95, see Table S2). With the trained RF model, we also derived the PDP to show the marginal effect of each predictor on yield shock probability. Overall, the same predictor might have quite different PDPs across the four regions, suggesting spatially varying sensitivities of yield shocks to climate and also justifying our partitioning of Europe into several regions. Based on PDP, ZHU ET AL.
10.1029/2020EF001815 6 of 13 Year Area fraction of yield shock 2010 1990 we found a greater f wd in a given region generally corresponded to a higher probability of yield shock, especially for extreme temperature in the reproductive period ( Figure 3). As f wd collectively accounted for >30°C heat stress or >90th percentile abnormal warming, the positive relationship suggested both heat stress and abnormal warming should be carefully considered in crop models to predict yield shock. For water supply (precipitation or f rd ), RF models suggested both low and excessive water supply can result in higher yield shock probability, especially for those in winter and in the reproductive period in North and West Europe (Figure 3). In terms of water demand (PET), higher PET generally corresponded to higher yield shock probability, especially in the vegetative and reproductive periods. We also noted that although water demand was usually highly correlated with warming, yield shock was less responsive to water demand relative to warming (Figure 3). The yield shock responded only weakly to cold stress, except for the winter cold stress in North and East Europe.
On the other hand, compared with the PDP of RF models based on official yield surveys, the crop model derived PDP suggested an overall positive but a weaker response of yield shock probability to f wd ( Figure S5). For water supply, lower water supply resulted in higher yield shock probability, especially for precipitation ( Figure S5). However, the higher yield shock probability with excessive water supply as observed in Figure

Attribution of Yield Shock Probability to Individual Climate Drivers
With yield shock probability being attributed to individual predictors (Figure 4a), the yield shocks primarily driven by the first four climate drivers in Table 2 comprised 83% of total yield shocks across Europe (cold stress was not presented, as cold stress driven yield shocks comprised less than 5% of total yield shocks). Across all studied countries, the attribution analysis suggested that low water supply was the most pervasive primary climate driver of shocks, followed with extreme warming, high water demand and excessive water supply (Figure 4a, although the dominant climate driver varied for each region). There was also a distinct spatial disparity about low water supply driven yield shock; the relatively drier and warmer regions (South and East Europe) were more dominated by low water supply than North and West Europe, and vice versa for the excessive water supply (Figure 4a). Seasonally, reproductive period was the critical stage for extreme warming driven yield shocks; however, for shocks driven by low water supply and high water demand the most important season was vegetative period, which was probably related with the lasting effect of water supply (Zhang & Oweis, 1999).
When the total time period was divided into two periods 1980-1999 and 2000-2018, for the whole Europe there was a distinct shift, such that the most pervasive emerging primary driver changed from low water supply to extreme warming from the earlier period to the later period (Figures 4b  and 4c). However, low water supply in South Europe remained as the most pervasive primary driver. Seasonally, the two time periods had a similar pattern to the whole time period; reproductive period was the most important stage for extreme warming and vegetative period was the most important for low water supply and high water demand. Further, the crop models overrepresent the effect of low water supply on yield during the reproductive period compared to the yield survey data Figure 5.

Change in Climate Drivers of Wheat Yield Shocks due to Climate Change
We projected yield shock probability to two climate scenarios (RCP4.5 and RCP8.5) using the RF models trained with historical survey yield data. As expected, future warming resulted in higher yield shock probability and RCP8.5 showed even higher shock probability with RCP8.5 representing the higher emission scenario and warmer climate than RCP4.5 ( Figure 6). Spatially, East European countries like Hungary and Romania were more prone to be influenced by yield shock in both scenarios. When the top 90th percentile of projected yield shock probability was identified as a yield shock event, the attribution analysis suggested that extreme warming was the most pervasive primary driver in both RCP4.5 and RCP8.5, while shocks driven by low water supply decreased significantly (Figure 7). This shift in yield shock drivers could be regarded as the extension of climate driver shift since 2000 as identified in Figure 4. Out of the 4 regions, East Europe had the most area with extreme warming as the primary driver in both scenarios (Figure 7), which might explain why East Europe had a higher yield shock probability in Figure 6.
ZHU ET AL.

Yield Shock Prediction Agreement Between RF Models and Crop Models
Generally, crop models had a satisfactory performance for predicting yield shocks primarily driven by low water supply during 1980-2010. This category of yield shock constituted 32% of total yield shocks and more than 60% of them could be predicted by more than 3 out of 6 crop models during 1980-2010 ( Figure 8a). However, there was a lower model agreement for other climate stress driven yield shocks, especially for the yield shock driven by autumn and winter (S1) stresses. For example, 73% yield shocks primarily driven by vegetative and reproductive period (S2) extreme warming could be predicted by more than 3 out of 6 crop model, while only 17% for yield shocks primarily driven by autumn and winter (S1) extreme warming ( Figure 8a). Yield shock primarily driven by excessive water supply had a low model agreement for both stages (Figure 8a), which could be inferred from crop model response curves where yield shock probability was unresponsive to excessive water supply ( Figure S5). Relative to the historic period, the model agreement became lower in future climate for all kinds of climate driven shock events (Figures 8b and 8c). This implies that in the future it might be even more challenging to forecast yield shock and design adaptation strategies accordingly using existing crop models as predictive tools.

Discussion and Conclusions
Here we built a data-driven attribution framework to investigate the primary climate drivers that may contribute to yield shocks in both historical and future climate scenarios. The attribution analysis suggests that for all studied areas low water supply was the dominant driver for yield shock. However, the most dominant climate drivers varied for each region: North Europe was dominated by excessive water, West Europe by excessive water, East Europe by high water demand, and South Europe by low water supply. This suggests that agricultural adaptation strategies should be tailored to meet the regional environment. Future warmer ZHU ET AL.
10.1029/2020EF001815 9 of 13  Unlike most studies where a machine learning algorithm was used to directly model yield variation (Crane-Droesch, 2018), here we classified yield data into yield shock and non-yield shock with the lowest tenth percentile yield as threshold. Then yield shock and non-yield shock was treated as a binary variable and a RF model was trained to classify them. This treatment allowed the RF model to focus on the extreme low yield cases while minimizing the influence of small yield variation on the estimation of yield response to climate. High AUC (>0.85) based on out-of-sample prediction suggested that the RF model produced a reliable classification. We also showed that the attribution analysis was robust to an alternative yield time series detrending method using Savitzky-Golay smoothing function ( Figures S8-S9). We should note that although we divided Europe into 4 groups, our model was unable to capture the finer spatial variation of farmer-adapted management practices in each administrative unit, like shift in sowing date and irrigation implementation. For example, irrigation is expected to expand to maintain crop yield at future warmer climate (Peña-Arancibia et al., 2020); however, our projection did not consider the potential expansion of irrigation area in the future, which might overestimate the risk of future yield shock.
On the other hand, our selection of four regions is based on the relatively homogeneous weather and geographic patterns within each region. Dividing the whole data set into more regions might better account for the spatial heterogeneity but will also reduce the training sample size for each regional RF model. We also conducted a sensitivity test through using an alternative spatial aggregation method (Table S3). The attribution analysis based on this second spatial aggregation ( Figure S10) also suggested that the major driver of yield shock during historical period is water limitation, although the climate drivers in each growing period differ between the two aggregation methods. For example, the extreme warming stress in Eastern Europe becomes more dominant yield shock driver in the alternative aggregation method ( Figure S10). We looked into this and find it is probably because yield shock probability shows a steeper response to warming stress after removing Austria in the second spatial aggregation ( Figure S11). This sensitivity test suggests that more attention need to be paid to spatial aggregation if more detailed administrative units become available in the future.
Since we did not account for the dynamics of crop phenology, which might reflect the farmer adaptation practices to climate change, we conducted a sensitivity test to see whether our conclusion will change with shifted crop phenology. When we shifted the crop phenological date by two weeks (±14 days), our main conclusion still holds, although a few differences were noticed in certain growing periods of climate drivers ( Figure S12). Although our results are robust, we suggest that the dynamic information of crop phenology Figure 8. Yield shock prediction agreement between RF models and crop models. There are totally 6 crop models reported future yield simulations. Therefore, for each Shock RF , it might be identified as yield shock by 0, 1, 2, …, 6 crop models. Shock RF was grouped based on its primary driver on different seasons (S1 means the early season: autumn and winter; S2 means the later season: vegetative period and reproductive period). The numbers on the grid represents the percentage of yield shocks also identified as shock (Shock sim ) by 0, 1, 2, …, 6 crop models, respectively. For example, for N yield shocks (Shock RF ) primarily driven by extreme warming in S1, there could be N0, N1, …, N5, N6 (N0 + N1 + … N5 + N6 = N) yield shocks also identified as shock (Shock sim ) by 0, 1, 2, …, 6 crop models, respectively. A higher value of N3 + N4 +N5 + N6 suggested a higher confidence to predict the yield shock. RP, reproductive period. S2 S1 S2 S1 S2 S1 % should be considered when available to better account for the influence of farmer management practices on crop yield.
Our analysis confirmed that accounting for the accurate timing information of climate stress is important, which can be inferred from the different responsive curves for the same climate variable in different seasons. On the other hand, lower crop model predictive power for autumn and winter stress relative to the later growing season suggested that crop models need to improve the simulation of abnormal weather effects during these seasons. Studies have revealed that abnormal warm winter resulted in insufficient vernalizing days and affected wheat leaf and tiller number or flowing timing (Kirby, 1992;Wu et al., 2017), which is likely related to the recent yield shock in the breadbasket of France (Ben-Ari et al., 2018). In addition, the excess water stress was always poorly represented by crop models, which calls for incorporating water excess effects in crop models to improve their predictive performance (Li et al., 2019;Rosenzweig et al., 2002).
We note that yield shock could be influenced by single or compound climate extremes as well as non-extreme variability (Anderson et al., 2019;Beillouin et al., 2020;Ben-Ari et al., 2018;Ciais et al., 2005). The traditional yield loss attribution analysis based on crop models requires a set of model experiments wherein certain crop stress schemes are turned on or off and then their outputs are compared (Webber et al., 2018). However, climate stresses are often interrelated and the interaction effect between water and heat stress could be also important for determining crop yield (Jin et al., 2017;Zandalinas et al., 2018). Omitting either of them could lead to unintended results. Our data-driven attribution framework provides a flexible alternative, as it does not require a priori knowledge on crop physiological response to climate stress and can handle nonlinear responses. The interaction effects of climate stresses were implicitly considered in the RF model through the hierarchical classification tree of predictors within the structure of the model. Then, the game theory based attribution gave a reasonable quantification of the relative importance of climate stresses for either single or compound extremes.
Overall, our study analyzed the spatial-temporal pattern of wheat yield shock and applied a data-driven attribution framework to quantify primary climate drivers of yield shock during the historical period. Our attribution analysis suggests that the overall dominant climate driver of yield shock during 1980-2018 is water limitation. However, it was surpassed by extreme warming in the late period (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). The future projection based on a calibrated RF model revealed enhanced yield shock probability and a paradigm shift from water limitation dominated yield shock to extreme warming dominated yield shock with warmer climate. However, crop models generally attributed low water supply as the primary climate driver of yield shock and were unresponsive to excessive water supply. The subsequent yield shock agreement analysis suggested that yield shocks driven by autumn and winter season extreme warming, high water demand, and excessive water supply were not as easily predicted by the crop models compared to the RF model. The future warmer climate resulted in even lower performance of crop models in predicting yield shocks. This analysis suggests that these climate stress effects during specific growth stages should be improved in future crop models to better develop climate adaptation strategies in cropping systems.