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

2.1.1. Gridded crop model data for emulator construction

2.1.2. Data for emulator-based yield projections

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.

Zoom In Zoom Out Reset image size

2.2.1. Development of annual GGCM emulators by extreme gradient boosting (XGB)

2.2.2. Emulator-based wheat yield projections

2.2.3. 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:

$$\begin{equation}{y_t} = \frac{{\sum\limits_{i = 1}^n {{y_{i,{\text{firr}},t}} \cdot } are{a_{{\text{firr}}}} + \sum\limits_{i = 1}^n {{y_{i,{\text{noirr}},t}} \cdot } are{a_{{\text{noirr}}}}}}{{\sum\limits_{i = 1}^n {\left( {are{a_{{\text{firr}}}} + are{a_{{\text{noirr}}}}} \right)} }}.\end{equation} \tag{ 1 }$$

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,{\text{firr}},t}}$ is the emulator-projected yield under fully irrigated conditions in grid cell i, and ${y_{i,{\text{noirr}},t}}$ is the emulator-projected yield for rainfed conditions in grid cell i; ${\text{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 ${\text{CV}} = \sigma /\mu$, in which $\sigma$ and $\mu$ are the standard deviation and mean, respectively, over a reference period. We compare the baseline period (1976–2005) with six future scenario-periods: 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:

$$\begin{equation}{\delta _{{\text{scenario}}}} = \frac{{{\text{C}}{{\text{V}}_{{\text{scenario}}}} - {\text{C}}{{\text{V}}_{{\text{baseline}}}}}}{{{\text{C}}{{\text{V}}_{{\text{baseline}}}}}} \times 100\% .\end{equation} \tag{ 2 }$$

2.2.4. The effects of changes in temperature, precipitation, and CO₂

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₂ fertilization effect ('T + P + CO₂'), 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₂ fertilization effect ('T + P + CO₂'), 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₂ 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).

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

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⁻¹), experienced the largest relative increase in wheat yield CV (figure 4).

Zoom In Zoom Out Reset image size

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.

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).

Zoom In Zoom Out Reset image size

After separating the relative contributions of climate drivers (without elevated CO₂ 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).

Zoom In Zoom Out Reset image size

The elevated CO₂ 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).

Zoom In Zoom Out Reset image size

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. 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_gsinterV, figure 8(c)), intra-seasonal variability (T_gsintraV, 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_gsmean > 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_gsmean < 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_gsmean, figure 8(f)), but positively correlated with interannual variability of precipitation (P_gsinterCV, figure 8(g)), and drought intensity (consecutive drought days, CDD_gs, figure 8(h)), all statistically significant.

Zoom In Zoom Out Reset image size

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₂ 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 of change for other sites. Spatial heterogenous yield variability changes are the main cause of inconsistence between our projections and Liu et al (2019). This is the case in Glen and Bloemfontein in South Africa, and Dharwar in India for the four inconsistent sites, as well as 11 out of the other 21 inconsistent sites, including those in northwest US, around the western or northern coast of the Black Sea, in central southern Russia, North China, and south-eastern Australia. Besides, another cause may be the choice of the crop model ensemble and underlying uncertainties. We examined our results for each crop model emulator. The direction of each site-based change in yield CV reported by Liu et al (2019) can be found in the result of at least one of our emulators, indicating GCMs-crop models ensemble combination is critical to yield projection.

Spatial heterogeneity of crop yield variability creates a huge challenge for agricultural risk management (Benami et al 2021). The spatially heterogeneous yield CV changes are also found in earlier reports that yield variability changes in rice and wheat are sensitive to spatial resolution (Iizumi et al 2018). Previous yield variability projections conducted with site-based, process-based models have found that regional yield variability changes are not consistent across different sites (Liu et al 2021). Thus, the gridded process-based crop models can provide an overview of global or regional changes in yield variability (Ostberg et al 2018, Parkes et al 2018). The ability to represent this spatial heterogeneity in yield variability in light-weight emulators allows for more comprehensive assessments of the risk of changes in 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₂ 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₂. Wheat as a C3 crop is known to have a high capacity to benefit from elevated CO₂ levels, which has been confirmed by various previous experiment-based evidences (Kimball 2016, Toreti et al 2020). Negative effects of global warming on future wheat yield could potentially be fully compensated by yield-amplifying effects of elevated atmospheric CO₂ concentrations (Ye et al 2020).

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₂ 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, process-based 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 non-climatic 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).

The spatial scale of our estimates reached the sub-province 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 mean-increasing 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 inter-temporal 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.