Future climate change significantly alters interannual wheat yield variability over half of harvested areas

Climate change affects the spatial and temporal distribution of crop yields, which can critically impair food security across scales. A number of previous studies have assessed the impact of climate change on mean crop yield and future food availability, but much less is known about potential future changes in interannual yield variability. Here, we evaluate future changes in relative interannual global wheat yield variability (the coefficient of variation (CV)) at 0.25° spatial resolution for two representative concentration pathways (RCP4.5 and RCP8.5). A multi-model ensemble of crop model emulators based on global process-based models is used to evaluate responses to changes in temperature, precipitation, and CO2. The results indicate that over 60% of harvested areas could experience significant changes in interannual yield variability under a high-emission scenario by the end of the 21st century (2066–2095). About 31% and 44% of harvested areas are projected to undergo significant reductions of relative yield variability under RCP4.5 and RCP8.5, respectively. In turn, wheat yield is projected to become more unstable across 23% (RCP4.5) and 18% (RCP8.5) of global harvested areas—mostly in hot or low fertilizer input regions, including some of the major breadbasket countries. The major driver of increasing yield CV change is the increase in yield standard deviation, whereas declining yield CV is mostly caused by stronger increases in mean yield than in the standard deviation. Changes in temperature are the dominant cause of change in wheat yield CVs, having a greater influence than changes in precipitation in 53% and 72% of global harvested areas by the end of the century under RCP4.5 and RCP8.5, respectively. This research highlights the potential challenges posed by increased yield variability and the need for tailored regional adaptation strategies.


Introduction
Interannual crop yield variability is one of the primary drivers of food system instability (IPCC 2019).Assessing the effects of climate change on yield variability is critical to understanding the impact of climate change on food security (FAO 2019).Due to trends in global warming (Lobell et al 2011) and the changing frequency and intensity of climate extremes (Trnka et al 2014), potential decreases in the mean yields of crops and an increase in the interannual yield variability could adversely affect the livelihoods of producers, create spikes in food prices, lead to hunger (IPCC 2014), and even cause political instabilities at a regional level (Sternberg 2011).Previously, the impact of climate change on mean crop yield (Lobell et al 2011, Rosenzweig et al 2014) has been investigated with a focus on food availability (Campbell et al 2016).From a climate risk perspective, the concept of time of climate impact emergence has recently been introduced, linking mean yield changes with historical yield variability (Jägermeyr et al 2021).Yet, the impact of climate change on future interannual yield variability has not received sufficient attention (Wheeler et al 2013, Challinor et al 2014).
Interannual yield variability has always been one of the key risk indicators of crop production.Early studies have either assumed a stationary process without considering variability changes (Ray et al 2015, Ceglar et al 2016, Tao et al 2016, Matiu et al 2017) or linked changes in variability to non-climatic factors (Kucharik and Ramankutty 2005, Döring and Reckling 2018, Knapp and Van Der Heijden 2018, Müller et al 2018).Recent studies have provided evidence for changes in the interannual yield variability of major cereal crops and identified significant impacts of climate change at the global scale 0.5 • grid level or at the country level (Osborne and Wheeler 2013, Iizumi and Ramankutty 2016).These studies have been followed up by regional, county-level analyses of the interannual yield variability of maize (Hawkins et al 2013, Lobell et al 2014, Leng 2017).Efforts have also been devoted to projecting the impact of future climate change on interannual yield variability, focusing on wheat and maize at global and regional scales, using process-based crop models (Moriondo et al 2011, Liu et al 2019) and statistical models (Urban et al 2012, Ben-Ari et al 2018, Tigchelaar et al 2018).Results from these studies have indicated substantial changes in interannual yield variability as a result of climate change, and that the sign and magnitude of change varies by production region.
Climate-related risk assessment on crop yield requires reflecting the spatial heterogeneity of both agricultural systems and climate change effects relevant for interannual yield variability (Benami et al 2021).There are still major research gaps in our understanding of these linkages across regions.In terms of major staple crops, only changes in yield coefficient of variation (CV) of wheat (Liu et al 2019) and maize have been analysed (Tigchelaar et al 2018) at the global scale.As these studies have used either site-based simulation or globally homogeneous warming perturbations, it is difficult to deduce robust conclusions on changes in interannual yield variability, reflecting the spatial heterogeneity of climate projections (Leng and Hall 2020).In addition, although the mechanism of impact and the mean yield response to change in climate drivers (e.g.temperature, precipitation, and CO 2 ) have been intensively discussed (Schlenker andRoberts 2009, Zhu et al 2019), the response of interannual yield variability to changes in the various climate drivers is not well understood.
The aim of this study is to evaluate potential changes in interannual wheat yield variability under two climate change scenarios globally, and to attribute individual contributions of temperature, precipitation, and CO 2 .The main research questions are: (a) How could climate change affect interannual wheat yield variability on current wheat-growing areas by the end of the century?(b) How much of these changes can be attributed to changes in temperature, precipitation, and their interaction, respectively?(c) To what extent can elevated CO 2 concentrations mitigate potential increases in yield variability? Answers to these questions will provide crucial information for climate risk assessment and effective adaptation measures.
We address these questions by conducting multimodel ensemble simulations with crop model emulators forced with global climate projections at high spatial resolution (0.25   Müller et al 2021).Emulators substantially improve computational efficiency and reduce data-processing requirements compared to running the original models, without sacrificing much prediction performance (Blanc and Sultan 2015, Blanc 2017, Folberth et al 2019, Franke et al 2020a, Ringeval et al 2021).The use of a large ensemble of GCM projections in combination with the ensemble of crop yield emulators allows for comprehensively evaluating changes in future yield variability and the associated distribution of extreme yield levels.

Gridded crop model data for emulator construction
The input and output data for the simulation of global gridded wheat yield were obtained from the GGCMI phase 2 experiment dataset (Franke et al 2020b).The spatial resolution of this dataset is 0.5  , including daily maximum and minimum temperatures, precipitation, and solar radiation (Ruane et al 2015).Based on these baseline reference simulations, the GGCMI phase 2 experiment used systematic perturbations in each grid cell with seven temperature levels (from −1 K to +6 K in 1 K interval, with +5 K skipped), nine precipitation levels (from −50% to +30%, in 10% interval, with −40% skipped), four CO 2 -concentration levels (360, 510, 660, and 810 ppm), and three nitrogen levels (10, 60, and 200 kg ha −1 ) (table S2; Franke et al 2020b).Twelve GGCMs were then forced with each of these perturbations of the original reanalysis weather data.The GGCMs used a national and subnational crop calendar for wheat that is based on Sacks et al (2010), Portmann et al (2010), and environmentbased extrapolations (Elliott et al 2015).
The output data contained irrigated and rainfed yield simulations from 1980 to 2010 for each of the different perturbation levels.In this study, we selected 8 out of the 12 crop models in the GGCMI phase 2 experiment for constructing the emulators.These were APSIM-UGOE, EPIC-IIASA, EPIC-TAMU, GEPIC, LPJ-GUESS, LPJmL, pDSSAT, and PEPIC.CARAIB was not included as it did not consider nitrogen stress.ORCHIDEE-crop was not included as it did not provide simulation results for spring wheat.PROMET and JULES were not included as they used different climate inputs.Although these eight crop models differed in their representation of crop phenology, leaf-area development, yield formation, root expansion, and nutrient assimilation, all accounted for the effects of water and heat stress and assumed no technological change (Blanc 2017).All input and output data sets were provided by GGCMI at the standardized spatial resolution of 0.5 • .More detailed descriptions of the individual crop models and the input and output data characteristics are available in the supplementary material (SM).

Data for emulator-based yield projections
To project a high spatial resolution global wheat yield, the Earth Exchange Global Daily Downscaled Projections (NEX-GDDP) dataset (Thrasher 2012), with a spatial resolution of 0.25 • , was obtained from the National Aeronautics Space Administration (NASA).This database contains the global daily maximum/minimum near-surface air temperature and precipitation data from 21 General Circulation Models (GCMs) from the Coupled Model Intercomparison Project phase 5 (CMIP5, Taylor et al 2012) under two representative concentration pathways (RCP4.5 and RCP8.5), covering the years 1950-2100.Other RCPs are not available through NEX-GDDP.
The emulator-based projections used a national and subnational crop calendar for wheat from MIRCA2000 (Portmann et al 2010).Given that MIRCA2000 has only monthly resolution, it was assumed that the first day of the month was the date of planting, and the last day of the month was the date of harvesting (Elliott et al 2015).The calendar we used to project yield was only MIRCA2000 because if we used the calendar of the GGCMI phase 2, we would be troubled with the mismatch between the separated spring and winter wheat calendar and only wheat harvested areas in spatial production allocation model (SPAM).Global wheat harvested area distribution around the year 2005 was obtained from the SPAM for rainfed and irrigated systems at five arc-minute resolution (You et al 2014).

Methods
The methodologies for evaluating changes in wheat yield variability under future climate scenarios includes the following steps (figure 1): (a) Develop annual yield emulators for the process-based GGCMI crop models.(b) Conduct emulator-based yield projections based on the NEX-GDDP climate model ensemble.(c) Summarize the future changes in wheat yield variability relative to the baseline; decompose the changes in yield variability into changes in mean yield and yield standard deviation.And (d) Separate the contributions from the changes in climatic drivers to the changes in the yield variability.

Development of annual GGCM emulators by extreme gradient boosting (XGB)
A previous study developed emulators of climatological-mean yield based on GGCMI phase 2 experiment data (Franke et al 2020b).We, however, develop an emulator capable of capturing year-toyear variability in yield.A machine-learning (ML) approach was used in this study for its flexibility for data-driven development of models with high accuracy (Folberth et al 2019) and its associated computational efficiency.
Development of the emulator consists of training-via a ML algorithm-on specific GGCM input and output datasets, so that the emulator replicates the complex process of yield simulation within the crop model.Variables that have been frequently reported to significantly influence wheat yield were prepared as the predicting variables, including climate, soil type, length of growing season, and management practices (table S3).All the data for training were computed/adapted from the GGCMs' input and output datasets.
All prediction variables were computed/obtained from the GGCMI phase 2 data archive.The climate data are supplied as daily values and were, in a first step, aggregated to monthly sums or averages (MON).For each grid, the month of planting was defined as month 1 to harmonize, on a global basis, the order of months from planting.Subsequently, prediction variables were calculated for each month in the growingseason months and the entire growing season (GS, based on the planting and harvesting dates for the GGCMs).Soil properties were adopted primarily for the topsoil class.Additional characteristics like length of growing season were regarded as a cultivar characteristic.The total amount and fraction of the nitrogen fertilizer application and CO 2 concentration were uniform for each grid.
In total, 32 different emulators were trained for the eight GGCMs, each with two water management modalities (rainfed and irrigation) and two wheat types (spring wheat and winter wheat).An XGB algorithm was used due to its better performance in terms of goodness-of-fit, cross-validation errors, and computation efficiency compared with a random forest algorithm (Folberth et al 2019).The predicting variables and the simulated yield in the GGCMs were randomly split into training and validation sets, which contained 75% and 25% of the samples (Yue et al 2019), respectively.Depending on the size of the dataset supplied by each GGCM, 1.7 × 10 6 -2.3 × 10 7 (0.9 × 10 7 -1.9 × 10 8 ) samples were used for model training and 0.6 × 10 6 -0.8 × 10 7 (0.3 × 10 7 -0.6 × 10 8 ) samples were used for validation of irrigation (rainfed) conditions, covering the period of 1981-2010.More details on emulator training, validation, and performance evaluation are available in SM and figures S1 and S2.

Emulator-based wheat yield projections
The emulators were then used to project global wheat yield by using future GCM projections.It is important to ensure that emulator-based projections do not exceed the range of training samples to avoid unrealistic extrapolation effects (Folberth et al 2019, Franke et al 2020b).The GGCMI phase 2 perturbations were designed to accommodate high-end warming scenarios (RCP8.5-2080s).For growing season average maximum temperature, the range of training data (interannual and spatial variability in AgMERRA + GGCMI perturbation) covered the entire range of the GCM projections.CO 2 concentrations were averaged with a 30 year moving window and the highest CO 2 concentration under RCP8.5-2080s is 760 ppm.

Climate drivers Descriptions 'T'
Using future scenarios of temperature, other drivers taken from baseline 'P' Using future scenarios of precipitation, other drivers taken from baseline 'T + P' Using future scenarios of temperature and precipitation, holding CO2 constant at 360 ppm 'T + P + CO2' Using future scenarios of temperature, precipitation, and CO2 Ensemble yield projections were conducted at the global level for grids with a spatial resolution of 0.25 • for the years 1976-2005 (baseline), 2006-2035  (2030s), 2036-2065 (2050s), and 2066-2095 (2080s) under RCPs of 4.5 and 8.5.If the spring and winter wheat are grown in parallel at national or subnational level, we determined the wheat type with larger harvested areas according to MIRCA2000 (Portmann et al 2010).The multi-model ensemble approach improves the robustness of future climate-change impact estimates and allows for analyses of spatial heterogeneity and inter-model uncertainty (Martre et al 2015).There were 336 future wheat yield estimates (21 GCMs × 8 emulators × 2 RCPs), each simulated for 4 × 30 year periods.Throughout all simulations, planting dates, cultivar selection, soil properties, and management practices were assumed to remain constant over time, which is consistent with the GGCMI Phase 2 experimental design A0 (Franke et al 2020b), which is used for training the emulators here.Final estimates of future yield responses are based on the median across the crop model emulators and GCMs.

Measuring the change in yield variability
Rainfed and irrigated yield were first aggregated to grid and national levels using an area-weighted average (Müller et al 2017), as described in the following equation: where i is the index of any grid cell assigned to the spatial unit in year t, n is the number of grid cells in that spatial unit, y i,firr,t is the emulator-projected yield under fully irrigated conditions in grid cell i, and y i,noirr,t is the emulator-projected yield for rainfed conditions in grid cell i; area is the harvested area in grid cell i, either due to fully irrigated or rainfed, obtained from SPAM.
We used the CV a measure of interannual yield variability, where CV = σ/µ, in which σ and µ are the standard deviation and mean, respectively, over a reference period.We compare the baseline period  with six future scenarioperiods: RCP4.5-2030s,RCP4.5-2050s,RCP4.5-2080s,RCP8.5-2030s,RCP8.5-2050s, and RCP8.5-2080s.The percentage change in yield CV in one of the six future scenario-periods relative to the baseline period is then measured by: (2)

The effects of changes in temperature, precipitation, and CO 2
The effects of changes in temperature, precipitation, and CO 2 were separated by using individual climate driver perturbed simulations, with one climate factor at a time taken from a climate scenario and the rest from the baseline.Four such climate driver sensitivity simulations (table 1) were conducted to isolate the effects of changes in temperature, precipitation, their interaction effects, and the CO 2 fertilization effect.The climate driver sensitivity simulations listed in table 1 allow for addressing the following: (a) The effects of temperature and precipitation changes can be derived by comparing the results of groups 'T' and 'P' with the baseline simulations, respectively.(b) The interaction between temperature and precipitation changes can be evaluated using the difference between groups 'T + P' and 'T' + 'P' .(c) The effect of CO 2 fertilization can be evaluated using the difference between groups, 'T + P + CO 2 ' and 'T + P' .

Global patterns of future change in wheat yield interannual variability
By the end of the century, model simulations indicate an overall decrease in wheat yield CV, but in some regions, including major producing countries, there would be more unstable wheat yield (figure 2).The spatial patterns of CV changes intensify towards the end of the century, indicating a more polarized pattern under the long-term scenarios RCP4.5-2080s and RCP8.5-2080s (figure S3).Under RCP8.5-2080s, with the CO 2 fertilization effect ('T + P + CO 2 '), the yield CV increases significantly in 18% of harvested areas (p < 0.05; see figure S4 for significance test), while 44% of harvested areas experience significant decrease of the yield CV (p < 0.05).Under RCP4.5-2080s, with the CO 2 fertilization effect ('T + P + CO 2 '), 23% of harvested areas undergoes significant increase of yield CV (p < 0.05), while yield becomes more stable in 31% of harvested areas significantly (p < 0.05).Western Europe, northern Australia, central US, South Asia, Southwest China, and Myanmar are found to experience a small increase in yield CV (<40%).In eastern Europe, southern Australia, and central India yield CV is indicated to decrease by >20% under RCP8.5-2080s(figure 2).The spatial patterns of changes are consistent across different scenarios and time periods, but the size of changes varies (figure S5).The uncertainty across crop yield projections (standard deviation of CVs across all 8 emulators and 21 GCMs) ranged between 17% and 119% with the CO 2 fertilization effect, with a global mean of 39% under RCP8.5-2080s.Uncertainty was most pronounced in central Europe to eastern Russia, and in the northern Indian production regions (figure S6).We further break the total uncertainty to those associated with the emulators and those with the GCMs, by analysis of variance.Disagreement across the emulators explained less than 50% of the total variance in 47% of the harvested areas (figure S7).
Changes in yield CV are linked to changes in mean yield and yield standard deviation.Under RCP8.5-2080s, mean yield levels increase in 92.1% of harvested areas and the yield standard deviation increases in 95.3% of harvested areas (figure S8).About 30.8% of the areas in which CV is found to increase, CV changes are dominated by increases in yield standard deviation (|SD+| > |MY+|).In regions where CV is decreasing, 59.3% of the areas are dominated by mean yield increases (|SD+| < |MY+|) (figure 3, table 2).Under RCP4.5-2080s, mean yield levels and the yield standard deviation increase in 92.8% and 94.5% of the harvested areas, respectively (figure S8).About 42.7% of the areas in which CV is found to increase, CV changes are dominated by increases in yield standard deviation (|SD+| > |MY+|).In regions where CV is decreasing, 47.6% of the areas are dominated by mean yield increases (|SD+| < |MY+|) (figure 3, table 2).

Changes in the yield CV across different climatic regions
Changes in the wheat yield CV exhibited a clear relationship with the baseline regional temperature, precipitation, and nitrogen fertilizer application rate.In general, regions with hotter growing seasons (growing season average temperature > 20 • C) or with lower nitrogen fertilizer application rates (nitrogen application rate < 200 kg ha −1 ), experienced the largest relative increase in wheat yield CV (figure 4).
The increases in yield CV tend to be greater in regions with hotter growing seasons under both RCP4.5 and RCP8.5, including sub-Saharan Africa, India, Australia's wheat belt, South East US, and southern Brazil and Argentina.These are regions in which mean wheat yields are expected to decrease under high-emission climate change scenarios, whereas at higher latitudes with lower growing season temperatures mean wheat yields are generally projected to increase (Jägermeyr et al 2021).The change in yield CV undergoes smaller decline under RCP8.5 and even experiences subtle increase under RCP4.5 in regions with wetter growing seasons, which can be attributed to stronger variability of precipitation in   wetter regions.Underperforming wheat production system regions, like Brazil, sub-Saharan Africa, and South East Asia, with lower levels of nitrogen application, are likely to experience a greater increase in yield CV under both RCP4.5 and RCP8.5.

Climatic drivers of and their relative contributions to the change in yield CV
In simulations based on individual climate drivers, temperature changes alone increase the yield CV for 55% and 56% of the harvested areas under RCP4.5-2080s and RCP8.5-2080s, respectively.Under RCP8.5-2080s the magnitude of increased yield CV with temperature change alone is larger than that with precipitation change alone, but the extent of the area affected by increasing yield CV is smaller (figure 5).The yield CV increases in 64% and 60% of harvested areas when only precipitation change is assumed under RCP4.5-2080s and RCP8.5-2080s, respectively (figure 5).After separating the relative contributions of climate drivers (without elevated CO 2 concentration) under RCP4.5-2080s,precipitation was the dominant driver to increase the yield CV in 33% of harvested areas, even if the temperature change plays a more important role in yield CV change in over half of harvested areas (53%).The interaction between temperature and precipitation change played a dominant role in changes in yield CV in 10% of harvested areas.Under RCP8.5-2080s, temperature becomes a more important factor and was found to be the dominant driver in 72% of global wheat harvested areas, of which, yield CV increased in 41% of harvested areas.Precipitation was found to be the dominant driver in 21% of harvested areas, of which, yield CV increased in 17% of harvested areas.The interaction between temperature and precipitation played a dominant role in the change in yield CV in only 8% harvested areas (figure 6).
The elevated CO 2 concentration reduced the increase in yield CV, which was greatest in RCP8.5-2080s.The effect was strongest (>15%) in central Europe, south Asia, North and Southwest China, and North America.The mitigation effect was weaker under the other RCP4.5-2080s, but the spatial patterns were largely consistent with RCP8.5-2080s (figure 7).
To elucidate the link between changes in yield CV and climate factors, we further examined the change in yield CV with climatic factors changes in mean, variability, and extremes of temperature and precipitation by using perturbation 'T' and 'P' results.

Changes in future wheat yield variability
Our results indicate that wheat yield CV might increase significantly in 18% of the global harvested area under a high-emission climate change scenario.In turn, yield variability is found to decrease in 44% of currently cultivated areas, regions in which mean yields are projected to increase under climate change.Globally, our findings are consistent with those of earlier studies indicating that declining yield variability is wide-spread but increasing yield variability is found across important breadbasket regions (Iizumi and Ramankutty 2016, Leng 2017).Site-based simulation results for a 2 • C warming scenario (Liu et al 2019) have provided a more pessimistic estimation, with wheat yield CV increases in 36 out of the 60 sites, including the CO 2 fertilization effect.Our results confirm higher yield variability in hot regions as reported by (Liu et al 2019) and in regions with low nitrogen fertilizer application rates as reported by (Han et al 2020).Similarly, the low yield CV in high nitrogen fertilizer application rates regions is consistent with the findings of nutrients-driven intensification that additional nutrient inputs raise mean crop yields and thus decrease yield CV (Müller et al 2018).
A detailed comparison of yield CV changes between site-based projections (Liu et al 2019) and our gridded projections demonstrates the importance of revealing spatial heterogeneity of yield variability changes.Changes in yield CV were identified as significantly increasing at all 14 sites for the 2 • C warming scenario (Liu et al 2019).Among the 14 sites, our estimates were consistent with 10 of the 14 stations.For the other 46 sites, our results are largely consistent with the 25 sites across the central U.S., South America, the Middle East, the western European coastline, and Southern Russia; but are different in the direction

Climatic drivers of changes in future wheat yield variability
The present results indicate strong links between changes in the wheat yield CV and changes in temperature and precipitation.Previous reports have suggested that changes in yield interannual variability are closely related to changes in the variability (both interannual and intra-seasonal) (Iizumi et al 2013, Peng et al 2018) and extremes of climate factors (Iizumi and Ramankutty 2016, Chen et al 2018).In addition, due to the non-linear relationship between yield and temperature, changes in the mean temperature, away from the optimal range, will increase the interannual yield CV (Urban et al 2012, Tigchelaar et al 2018).The response of the interannual yield variability to changes in precipitation is more complex than for temperature.In general, changes in precipitation have smaller effects on irrigated yield than on rainfed yield (Tubiello et al 2002, Kothari et al 2019).In rainfed systems, yield interannual variability has been known to be closely related to interannual variability of rainfall, as well as frequency and intensity of drought (Webber et al 2018).The effect of total precipitation change largely depends on the baseline humidity of the production region.For drylands, increasing total precipitation increases mean yield (Fronzek et al 2018) and consequently reduces CV.Also, the interaction between temperature and precipitation changes can mitigate the increase in yield CV, although the magnitude of the interaction effect on change in yield CV is modest, within 10%.This is similar to the mitigation effect of irrigation on heat stress (Zaveri and Lobell 2019).However, the interaction effect cannot be explicitly explained, depending on the timing, intensity, and volume of rainfall (Tack et al 2017).
Higher atmospheric CO 2 concentrations mitigate variability changes in crop yield (Urban et al 2015), a consistent finding across different scenarios and time periods (figure 7).The mitigation effect is mainly attributed to increases in mean crop yield under elevated CO 2 .Wheat as a C3 crop is known to have a high capacity to benefit from elevated CO 2 levels, which has been confirmed by various previous experimentbased evidences (Kimball 2016, Toreti et al 2020).Negative effects of global warming on future wheat yield could potentially be fully compensated by yieldamplifying effects of elevated atmospheric CO 2 concentrations (Ye et al 2020).

Uncertainties
The spatial pattern of uncertainty in our results is consistent with different uncertainty distribution between high and low latitudes provided by crop model simulation (Xiong et al 2020).Including the CO 2 fertilization effect would further increase the total size of uncertainty in the projected yield.This is in agreement with a recent analysis on sources of uncertainty regarding GCMs and GGCM statistical emulators (Müller et al 2021).
The use of emulator ensemble simulation enabled the estimation of wheat yield variability change driven by climate change.Nevertheless, our approach has two limitations.First, crop damage from climate extremes is a major driving force of interannual yield variability (Trnka et al 2014), but the capability of most crop models in reproducing extreme climate damage to crops is still limited (Rötter et al 2018).For instance, process-based crop models of the GGCMI phase 1 experiment fail to reproduce yield impact from too wet conditions (Li et al 2019).Also, processbased crop models underestimate the extremeness of the 2003 heat-drought (Schewe et al 2019).We employ newly developed crop model emulators to project future wheat yield and these emulators are capable of capturing the direction of yield anomalies due to climate extremes, indicating the type of extreme event-induced yield variability that is captured by the models (heat, drought) will increase yield variability in a fair share of current cropland.Second, interannual yield variability driven by non-climatic factors is not considered in our analysis.These nonclimatic factors can strongly affect yield variability (Albers et al 2017) and changes in management can also strongly affect yield levels under climate change (Minoli et al 2019, Zabel et al 2021).

Implications
The spatial scale of our estimates reached the subprovince scale in China and the sub-state scale in the US and thus provided more insight than previous global estimates.First, gridded estimated yield variability change could provide more detail on spatial heterogeneity in local areas.Such local spatial differences were pronounced in South Africa, eastern Africa, and Central Russia.Second, when there is a need to estimate regional or country-level aggregated yield variability change, our gridded estimates could enable straightforward aggregation rather than upscaling from site-based estimates-these estimates rely heavily on the representativeness of sites.
High-spatial resolution gridded estimates of future yield variability change enabled global estimates of future change in yield CV. Globally, changes in yield CV tend to decrease in 44% of global harvested areas; but still yields would become more unstable in 18% of global harvested areas under RCP8.5-2080s,including several major production regions and countries.This indicates potential challenges to the stability of grain supply, market price, and consequently, the whole food system in the context of future climate change.It is important for local and regional economies to proactively implement adaptation measures and policy support (Iizumi and Ramankutty 2016).In light of this, our results can provide details of spatial heterogeneity in local areas and identify regions with urgent needs, including those hot and low-fertilizer application regions.The predominant climate driver is also identified, so that adaptation strategies can be tailored for regional or local challenges.
To face the challenge of increased yield interannual variability, adaptations including meanincreasing and variance-reducing strategies (Mehrabi and Ramankutty 2019), are needed because the changes in relative yield variability (CV) are sourced from both changes in mean yield and yield standard deviation (figure S9).The focus on the relative yield variability (CV) rather than absolute (e.g.SD) reflects the producer perspective, where the variability around the mean is relevant (storage and financial buffers) even if the mean is increasing in the long-term (Hasegawa et al 2021).Shifting cultivars (Liu et al 2010, Olmstead and Rhode 2011), changing sowing regions (Iizumi et al 2021) and adjusting planting dates (Lobell 2014, Huang et al 2020) have been recognized as effective adaptation options to address heat stress.Likewise, reinforcing irrigation equipment and adjusting irrigation strategies could relieve water shortages (Zhao et al 2020).Additionally, increasing nitrogen application rates in underperforming wheat production system regions could mitigate the increase in yield CV (Han et al 2020).From a risk management perspective, the risk in increased yield CV requires better domestic intertemporal reserves of wheat grain to smooth fluctuations in interannual production, market supply, and commodity price and better financial buffers at the producer level, to mitigate financial losses from local less-than-average yields.

Conclusions
This study presents one of the first projections of future wheat yield interannual variability change at high spatial resolution and disentangles the impacts from changes in temperature, precipitation, and CO 2 on those changes.Our results reveal that future climate change alters wheat yield interannual variability in over 60% of harvested areas.Wheat yield variability may decrease in over 40% of global wheat harvested areas under a high-emission climate change scenario (RCP8.5-2080s),while under RCP4.5-2080sonly 31% of harvested areas undergo the declined yield CV.However, 23% and 18% of harvested areas experience increased yield CV under RCP4.5-2080s and RCP8.5-2080s, respectively.Greater increase in yield standard deviation than that in the mean yield was the main reason for the increase yield variability under both RCP4.5 and RCP8.5.Yields in hotter or lower fertilizer regions are projected to become more unstable.Worldwide, changes in temperature have a stronger influence on changes in yield variability compared with precipitation in 72% of global harvested areas under RCP8.5-2080s,whereas under RCP4.5-2080s the areas controlled by temperature changes are smaller (predominant in 53% of harvested areas).The global mean of yield CV reduction due to rising CO 2 concentration across current harvested areas are 5% and 8% under RCP4.5-2080s and RCP8.5-2080s, respectively.High spatial resolution patterns of changes in wheat yield variability, as well as site-specific major driver identification results, have great implications for policy-making with regard to where food supply and farmer income need to be stabilized by additional measures in wheat production throughout the world.

Figure 1 .
Figure 1.Framework for evaluating changes in global wheat yield interannual variability (T: temperature, P: precipitation).

Figure 3 .
Figure 3. Factors of yield CV changes, including mean yield changes (MY) and yield standard deviation changes (SD).Positive changes are indicated with '+' and negative changes with '−' .The |MY+| and |SD+| denote the absolute increase of mean yield and yield standard deviation, respectively.
a '+' denotes positive changes, and the '−' denotes negative changes.'|SD+| > |MY+|' denotes the absolute value of increase in yield SD is greater than that in the mean yield.

Figure 4 .
Figure 4. Wheat yield CV change under RCP4.5 and RCP8.5, separated by different climatic and management bins: growing season mean temperature (a), growing season total precipitation (b), and nitrogen fertilizer application (c).The bin classification refers to baseline reference conditions.CV change is based on the T + P + CO2 simulations.Box-and-whisker plots show the distribution of yield CV changes across all cultivated grid cells in each class.The group divisions are based on approximately equal sample sizes.

Figure 5 .
Figure 5. Changes in yield CV under isolated temperature and precipitation perturbations (RCP4.5-2080sand RCP8.5-2080s).Yield CV changes are shown as the median of 21 GCMs and 8 crop model emulators.'T' is the effect of temperature change, and 'P' is the effect of precipitation change.

A
linear regression was conducted between median changes in yield CV and growing season climatic factors from food producing units (FPUs, Kummu et al 2010) (figure 8).The change in yield CV was positively correlated with change in interannual variability (T gs interV, figure 8(c)), intra-seasonal variability (T gs intraV, figure 8(d), and extreme degree day (EDD gs , figure 8(e)) of temperature.The relationship between mean temperature and yield CV varied by region.For regions with hotter growing seasons (T gs mean > 10 • C, figure 8(a)), a warming trend tended to increase the yield CV, and decrease the yield CV in regions with colder growing season (T gs mean < 10 • C, figure 8(b)).For the effect of precipitation change, results from grid cells with rainfed systems showed that change in yield CV was negatively correlated with change in total precipitation (P gs mean, figure 8(f)), but positively correlated with interannual variability of precipitation (P gs interCV, figure 8(g)), and drought intensity (consecutive drought days, CDD gs , figure 8(h)), all statistically significant.

Figure 6 .
Figure 6.Major contributors to the change in wheat yield CV (RCP4.5-2080sand RCP8.5-2080s).The dominant factors driving the change in yield CV are defined as the driving factors that contribute the most to the increase (or decrease) in the yield CV in each grid cell.The suffix '+' attached to the driving factor name indicates increase in the CV, whereas a '−' indicates a reduction in the CV.

Figure 8 .
Figure 8. Correlations of changes in yield CV with the changes in mean, variability, and extremes of temperature and precipitation at the FPU level relative to the baseline.In each panel, the changes in growing season climate factor are: (a) mean temperature (baseline value > 10 • C), (b) mean temperature (baseline value < 10 • C), (c) interannual variability (standard deviation) of mean temperature, (d) intra-season variation of daily temperature, (e) extreme degree days, (f) total precipitation, (g) interannual variability (coefficient of variation) of precipitation, and (h) consecutive drought days.
). Statistical crop model emulators are developed based on simulations from global gridded crop models (GGCMs), facilitated by AgMIP's Global Gridded Crop Model Intercomparison Project (GGCMI).Crop model emulators have recently gained popularity as a powerful tool for assessing the impact of climate change on crop yield (Lobell and Burke 2010, Holzkämper et al 2012, Oyebamiji et al 2015, Raimondo et al 2021, Franke et al 2020b).Baseline climate inputs were used from the AgMIP Modern-Era Retrospective Analysis for Research and Applications (AgMERRA) forcing dataset • .The input data included four different data types, i.e. climate, soil, atmospheric CO 2 concentration, and nitrogen fertilizer application rates (table S1 (available online at stacks.iop.org/ERL/16/094045/mmedia),

Table 2 .
Attribution of wheat harvested area with yield CV changes to changes in mean yield (MY) and standard deviation (SD) under RCP4.5-2080s and RCP8.5-2080s.

W
Liu et alHuang M, Wang J, Wang B, Liu D L, Yu Q, He D, Wang N and Pan X 2020 Optimizing sowing window and cultivar choice can boost China's maize yield under 1.5 • C and 2 • C global warming Environ.Res.Lett.15 024015 Iizumi T, Ali-Babiker I A, Tsubo M, Tahir I S, Kurosaki Y, Kim W, Wang E, Yang X and Wang J 2010 Contributions of climatic and crop varietal changes to crop production in the North China Plain, since 1980s Glob.Change Biol.16 2287-99 Lobell D B 2014 Climate change adaptation in crop production: beware of illusions Glob.Food Sec. 3 72-76 Lobell D B and Burke M B 2010 On the use of statistical models to predict crop yield responses to climate change Agric.For.Meteorol.150 1443-52 Lobell D B, Roberts M J, Schlenker W, Braun N, Little B B, Rejesus R M and Hammer G L 2014 Greater sensitivity to drought accompanies maize yield increase in the U.S. Midwest Science 344 516-9 Lobell D B, Schlenker W and Costa-Roberts J 2011 Climate trends and global crop production since 1980 Science 333 616-20 Martre P et al 2015 Multimodel ensembles of wheat growth: many models are better than one Glob.Change Biol.21 911-25 Matiu M, Ankerst D P and Menzel A 2017 Interactions between temperature and drought in global and regional crop yield variability during 1961-2014 PLoS One 12 e0178339 Mehrabi Z and Ramankutty N 2019 Synchronized failure of global crop production Nat.Ecol.Evol. 3 780-6 Minoli S et al 2019 Global response patterns of major rainfed crops to adaptation by maintaining current growing periods and irrigation Earth's Future 7 1464-80 Moriondo M, Giannakopoulos C and Bindi M 2011 Climate change impact assessment: the role of climate extremes in crop yield simulation Clim.Change 104 679-701 Müller C et al 2017 Global gridded crop model evaluation: benchmarking, skills, deficiencies and implications Geosci.Model Dev. 10 1403-22 Müller C et al 2018 Global patterns of crop yield stability under additional nutrient and water inputs PLoS One 13 e0198748 Müller C et al 2021 Exploring uncertainties in global crop yield projections in a large ensemble of crop models and CMIP5 and CMIP6 climate scenarios Environ.Res.Lett.16 034040 Olmstead A L and Rhode P W 2011 Adapting North American wheat production to climatic challenges, 1839-2009 Proc.Natl Acad.Sci.108 480-5 Osborne T M and Wheeler T R 2013 Evidence for a climate signal in trends of global crop yield variability over the past 50 years Environ.Res.Lett.8 024001 Ostberg S, Schewe J, Childers K and Frieler K 2018 Changes in crop yields and their variability at different levels of global warming Earth Syst.Dyn. 9 479-96 Oyebamiji O K, Edwards N R, Holden P B, Garthwaite P H, Schaphoff S and Gerten D 2015 Emulating global climate change impacts on crop yields Stat.Model.15 499-525 Parkes B, Defrance D, Sultan B, Ciais P and Wang X 2018 Projected changes in crop yield mean and variability over West Africa in a world 1.5 K warmer than the pre-industrial era Earth Syst.Dyn. 9 119-34 Peng B, Guan K, Pan M and Li Y 2018 Benefits of seasonal climate prediction and satellite data for forecasting U.S. maize yield Geophys. Res.Lett.45 9662-71 Portmann F T, Siebert S and Döll P 2010 MIRCA2000-global monthly irrigated and rainfed crop areas around the year 2000: a new high-resolution data set for agricultural and hydrological modeling Glob.Biogeochem.Cycles 24 GB1011 Raimondo M, Nazzaro C, Marotta G and Caracciolo F 2021 Land degradation and climate change: global impact on wheat yields Land Degrad.Dev.32 387-98 Ray D K, Gerber J S, Macdonald G K and West P C 2015 Climate variation explains a third of global crop yield variability Nat.Commun.6 5989 Ringeval B, Müller C, Pugh T, Mueller N, Ciais P, Folberth C, Liu W, Debaeke P and Pellerin S 2021 Potential yield simulated by global gridded crop models: a process-based emulator to explain their differences Geosci.Model Dev.Discuss.14 1639-56 Rosenzweig C et al 2014 Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison Proc.Natl Acad.Sci. 111 3268-73 Rötter R P, Appiah M, Fichtler E, Kersebaum K C, Trnka M and Hoffmann M P 2018 Linking modelling and experimentation to better capture crop impacts of agroclimatic extremes-a review F. Crop.Res.221 142-56 Ruane A C, Goldberg R and Chryssanthacopoulos J 2015 Climate forcing datasets for agricultural modeling: merged products for gap-filling and historical climate series estimation Agric. 