Research ArticleAGRICULTURE

# Spatial variations in crop growing seasons pivotal to reproduce global fluctuations in maize and wheat yields

See allHide authors and affiliations

Vol. 4, no. 11, eaat4517

## Abstract

Testing our understanding of crop yield responses to weather fluctuations at global scale is notoriously hampered by limited information about underlying management conditions, such as cultivar selection or fertilizer application. Here, we demonstrate that accounting for observed spatial variations in growing seasons increases the variance in reported national maize and wheat yield anomalies that can be explained by process-based model simulations from 34 to 58% and 47 to 54% across the 10 most weather-sensitive main producers, respectively. For maize, the increase in explanatory power is similar to the increase achieved by accounting for water stress, as compared to simulations assuming perfect water supply in both rainfed and irrigated agriculture. Representing water availability constraints in irrigation is of second-order importance. We improve the model’s explanatory power by better representing crops’ exposure to observed weather conditions, without modifying the weather response itself. This growing season adjustment now allows for a close reproduction of heat wave and drought impacts on crop yields.

## INTRODUCTION

Year-to-year fluctuations in crop yields and weather-induced shocks pose a considerable risk for food security, via both direct effects on local supply and prices, and remote effects via the global crop market (13), which can have severe impacts for producers and consumers (4, 5). Weather-induced agricultural losses belong to the most costly natural disasters (6, 7). For maize, the most important global staple in terms of production, the 2003 European heat wave resulted in yield reductions of 13% (8), and total direct economic losses in the agricultural sector are estimated at more than US$12 billion (9, 10). In 2012, drought conditions caused a decline in U.S. maize yields of >22% (11, 12), resulting in direct economic damages through crop losses alone of almost US$30 billion (7, 12). On average (1964–2007), global heat waves and droughts [EM-DAT record (13)] have reduced national reported maize yields by 12 and 7% compared to undisturbed reference conditions, respectively. Wheat yields, in turn, have been less sensitive to heat waves (4% reduction), but slightly more to droughts (7% reduction) (14).

Weather variability in general (1, 15, 16) and the frequency and intensity of heat waves and droughts are expected to increase under global warming (1719). It is thus essential to develop a detailed understanding of the associated yield responses to assess future risks to food production. Empirical models forced by observed weather fluctuations can, however, explain only about 30% of the observed yield variability at global average levels (2023). Process-based models, incorporating the mechanistic understanding of driving factors, achieve similar levels globally but are shown to explain >50% of the national average yield variability in individual countries (24, 25).

Both the process-based and empirical approaches leave substantial fractions of observed variances unexplained, which are likely due to changes in management decisions and impacts from weeds, diseases, and pests not represented in the models. But they may also be associated with an inappropriate or incomplete representation of the actual weather signal, the yield response to weather forcing, or reporting errors. A fundamental prerequisite to disentangle these confounding factors in global gridded crop models is a realistic representation of plant development dates and phenology, that is, the timing of planting, anthesis, and maturity (11, 26). Detailed spatially and temporally resolved knowledge about the exact growing season window would ensure the correct exposure to actual weather conditions, ruling out inadequate forcing that might cause models to obscure. However, global benchmarking of crop model responses to weather forcing is limited particularly because of missing information about management and other non–weather-related drivers (24, 25).

Here, we demonstrate to what degree the explanatory power of process-based yield simulations can be improved by comprehensively accounting for available local information about the average month of sowing and harvest for different crops—representative of average climate conditions around the year 2000 (27). We particularly test to what degree an adequate representation of spatial variations in growing seasons improves the reproducibility of observed global yield losses due to heat waves and droughts [derived from EM-DAT (13)], without modifying the model’s weather response itself. To this end, spatial variations in phenological heat units (PHU; in °C day), a cumulative measure of heat required to reach maturity, are used to implement information on the timing of maturity for different crops into a global gridded crop model (LPJmL; see Materials and Methods). Geographic variations in PHU can be considered as proxy for cultivar choices of the same crop in the model. We note that cultivars are not fully specified by heat sums but are characterized by additional properties (28), which are only partly resolved in global gridded crop models and not modified in our simulations. In general, however, cultivars grown in colder climates require smaller heat sums to reach maturity compared with those grown in warmer climates (29).

We compare the effect of improved growing season timing with model improvements regarding refined representation of water stress on irrigated land through accounting for irrigation water constraints, starting from the reference model version assuming unlimited water supply in irrigated systems. We also consider a basic version neglecting water stress through unlimited irrigation on all crop land. Simulation performances are evaluated on the basis of (i) the variance in reported national maize and wheat yield anomalies [Food and Agriculture Organization (FAO) (8)] that can be explained (“explanatory power”), (ii) the performance in reproducing reported composite yield impacts of heat waves and droughts at global scale, and (iii) the reproduction of the yield declines induced by the European heat wave 2003 as a case study. Analyses are focused on the group of main producer countries that collectively provide 90% of the respective global maize and wheat production (2000–2011 average; hereinafter “main producers”). Of these 20 countries for maize and 26 countries for wheat, we particularly focus on the 10 countries that exhibit highest weather sensitivity in yield anomalies, that is, where the maximum of the explained variances provided by the considered model versions solely forced by observed weather information is highest (hereinafter “most weather-sensitive main producers”).

While reported sowing dates have already been used to harmonize crop model simulations, for example, within the Global Gridded Crop Model Intercomparison (GGCMI) (24, 30) and other studies (11, 12, 31), individual models are rarely constrained to meet observed harvest dates and often assume static crop cultivars globally (24) or across large regions (32). Here, we make use of the spatially explicit information about sowing and harvest dates by crop to constrain global gridded crop model simulations and systematically quantify for the first time associated gains in model performance regarding (i) representing annual fluctuations in national maize and wheat yields and (ii) reproducing observed yield responses to heat waves and droughts.

## RESULTS

We use the global gridded crop model LPJmL (33) to calculate historical annual crop yields at 0.5° spatial resolution, forced by three observational climate datasets (3436). To study the contribution of both refined growing season timing and irrigation water constraints to increased explanatory power of crop model simulations, we evaluate the following three model versions (see Table 1 for overview and Material and Methods for details).

Table 1 Experimental design.

Global gridded crop model versions used in this study are characterized regarding assumptions on irrigation, phenological heat units (PHU), and sowing dates; see Materials and Methods for details.

View this table:

The reference model [1. LPJmL–Ref (37), as it is used for example, in GGCMI (24) or in the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) (25)] internally derives sowing dates from climate conditions, and maturity is reached at a simplistic PHU estimate with limited spatial variation (semistatic). Neither sowing nor harvest is adjusted to meet reported data. Irrigation is assumed to occur on areas currently equipped for irrigation (27) without accounting for system efficiencies and water constrains.

In a two-step approach, the model is first updated with an advanced representation of local surface-water constraints based on a detailed representation of irrigation systems and their efficiencies [2. LPJmL–WaterLimIrr (38)]. In a second step, growing seasons are additionally adjusted by prescribing grid cell–level sowing dates and PHU requirements to meet local average harvest dates, individually for rainfed and irrigated systems (3. LPJmL–PHU). Crop calendar information are derived from observational MIRCA2000 data (27), complemented by LPJmL-derived seasons in countries, where MIRCA2000 information leads to lower agreement with reported yield data.

To quantify the fundamental contribution of water stress to crop yield fluctuations, we introduce a fourth model run in which water stress is eliminated through unrestricted irrigation in both currently rainfed and irrigated systems (4. LPJmL–NoWaterStress). Assumptions herein on sowing dates and PHU are identical to the LPJmL–Ref simulation.

### Reproducing national crop yield variability

Accounting for water constraints in irrigated farming systems, and especially the representation of spatial crop cultivar variations, substantially increases the fraction of the variance in observed national crop yield anomalies that can be explained by the process-based model simulations forced solely by weather fluctuations (Fig. 1). For maize, the explained variance (R2) roughly doubles—from 29% (LPJmL–Ref) to 58% (LPJmL–PHU) averaged across the most weather-sensitive main producers (Fig. 1; 20 to 37% across all main producers). Differences [root mean square error (RMSE)] between observed and simulated standardized crop yield anomalies decrease from 0.95 to 0.67 (fig. S1; 1.06 to 0.87 all main producers). Respective time series are shown in fig. S2.

The effect of the improved representation of water constraints in irrigation—R2 increases from 29% (LPJmL–Ref) to 34% (LPJmL–WaterLimIrr) across the most weather-sensitive main producers (Fig. 1; 20 to 23% all main producers)—is consistently surpassed by contributions from improved timing of the growing season in LPJmL–PHU (R2 reaches 58% across most weather-sensitive main producers and 37% for all main producers).

Eliminating water stress by assuming full irrigation in LPJmL–NoWaterStress reduces the explained variances to only 10% (most weather-sensitive main producers, 6% all main producers), which highlights the fundamental importance of water stress for explaining interannual yield fluctuations. In turn, while the water stress representation (difference between LPJmL–NoWaterStress and LPJmL–WaterLimIrr) can increase the explained variance to 34% for the most weather-sensitive main producers (and 23% for all main producers), the growing season adjustment improves simulation performance to the same extent (to 58% for the most weather-sensitive main producers and 37% for all main producers; see Fig. 1 and fig. S3). Our results thus suggest that using more reliable information on the timing of the growing season improves maize-yield simulations as much as the fundamental representation of water stress itself, which is one of the main growth-limiting constraints in large-scale crop models.

For wheat, model improvements due to a refined representation of irrigation processes and adjusted growing seasons follow a similar pattern: Gains from refined growing-season timing exceed those from refined irrigation processes and surface water availability. Average explained variances increase from 45% (LPJmL–Ref) to 54% (LPJmL–PHU) for the most weather-sensitive main producers (Fig. 1; 25 to 34% for all main producers). LPJmL–PHU wheat simulations achieve largest relative gains in countries at the lower end of simulation performance (fig. S3), turning nonsignificant correlations in five of eight countries statistically significant (P <0.1; table S2). Differences (RMSE) between simulated and observed standardized yield anomalies decline from 0.8 to 0.72, respectively (fig. S1; 1.02 to 0.94 for all main producers). The LPJmL–NoWaterStress simulation does not exceed 9% explained variance on average.

In most countries, the gain in model performance due to the growing-season adjustment is dominated by the effect of spatially resolved representation of PHUs, while the selection of prescribed sowing dates only plays a minor role (fig. S4). Thus, LPJmL–PHU simulations achieve substantially better results compared to LPJmL–WaterLimIrr (that is, LPJmL–Ref phenology but water constraints as in LPJmL–PHU), even when the growing season is not based on observational data, but derived by the LPJmL–Ref model (R2 of 53% for LPJmL–PHU with model-derived growing season and 34% for LPJmL–WaterLimIrr, respectively, for maize across most weather-sensitive main producers; fig. S4). This is explained by excessive PHU requirements in the LPJmL–Ref phenology, which often lead to premature harvests when the maximal growing season length is reached (Materials and Methods). In LPJmL–PHU, this is avoided by associated spatial variations in PHU requirements that are based on local climate conditions. LPJmL–PHU based on MIRCA2000 growing seasons alone reaches R2 of 54% for maize across the most weather-sensitive main producers (fig. S4). The MIRCA2000 crop calendar outperforms the LPJmL-calculated seasons on 80% of the maize and 59% of the wheat cropland area, respectively (fig. S5)—predominantly in countries in which MIRCA2000 provides subnational data (for example, USA, China, India, and Brazil). Explained variances depend on the considered climate forcing dataset, but associated uncertainties are minor yet somewhat larger for maize than for wheat (fig. S6). The mean of simulations forced by different climate inputs achieves higher explanatory power than any of the individual simulations in many countries.

### Quantitative understanding of crop yield responses to heat waves and droughts

Worldwide heat waves between 1964 and 2007 reduced observed national maize and wheat yields on average by 12.4% (quartile range, −1 to −20%, referring to the 75th and 25th percentile of the normalized distribution of all 65 events) and 4.1% (+3 to −9%, 81 events), respectively, compared to the mean of the three before and following nonevent years (Fig. 2, see Materials and Methods for details on this compositing analysis). Throughout the same period, droughts led to mean yield declines of 7% (+2 to −16%, 175 events) for maize and 7.1% (+3 to −16%, 146 events) for wheat (table S3 lists all considered events). Neither heat waves nor droughts show a lagged yield level response in the year following reported events (Fig. 2).

Except for simulations neglecting water stress, all model versions agree on the sign of heat wave and drought yield influences on maize and wheat. Even the versions without growing-season adjustment (LPJmL–Ref and LPJmL–WaterLimIrr) are able to simulate global average effects of extreme events for wheat. However, they only reproduce about half of the observed decline in maize yields in response to heat waves. For droughts, LPJmL–Ref only reproduces half of the observed maize response, while water constraints in irrigation (LPJmL–WaterLimIrr) already lead to a close agreement with observations. Yet, this model version overestimates the drought effects for wheat yields (Fig. 2). Throughout all cases, the adjustment of growing seasons (LPJmL–PHU) leads to close reproduction of the observed losses across reported events, that is, in both maize and wheat yields and for heat waves and droughts. Improvements are particularly relevant for maize, with remarkably accurate replications of the yield signal also in the years following the extreme event (Fig. 2, A and B).

Simulated yield declines in the wake of heat waves and drought are primarily induced by water stress. Eliminating soil water deficits by assuming full irrigation, even in currently rainfed systems (LPJmL–NoWaterStress), leads to increases in simulated maize yields during heat wave and drought years (Fig. 2), which indicates that additionally available radiation in such years is, on average, beneficial for crop growth, as long as the water demand is fulfilled. Reported national yields (Fig. 2) include yields achieved under irrigation, because FAO yield statistics (8) are not available separately for rainfed systems. In most countries, the fraction of cropland under irrigation is low, that is, national yields are generally dominated by production on rainfed land. Further LPJmL–PHU model simulations confirm that at the global scale, heat wave, and drought events predominantly affect rainfed yields, with marginal influences in irrigated systems (heat wave effects on maize yields: −19.7% rainfed only, −12.4% mixed, and −0.6% irrigated only; fig. S7). In general, if water availability permits, then irrigation not only reduces water stress but also alleviates direct adverse extreme temperature impacts due to surface cooling. The latter effect, however, is not represented in the model and seems to be of minor importance for reproducing observed yield effects at national level (see Discussion).

### Refined crop exposure to rainfall deficit

Improved representation of growing season timing in the model changes the crop’s simulated exposure to rainfall and thus water stress, which is a key driver of yield fluctuations in all gridded crop models. For maize, the median rainfall deficit (growing season rainfall relative to long-term average growing season rainfall; growing seasons depend on respective model version) aggregated across heat wave years is considerably lower in the model with improved growing season timing (LPJmL–PHU) compared to the reference model (Fig. 3A). This is particularly pronounced in the case of the Europe 2003 heat wave, where the median deficit across 10 European countries decreases from −11 to −44%, respectively, between the two models (Fig. 3A). Maize yields declined by 12.7% during the 2003 heat wave and LPJmL–PHU simulations (10.8% decline) closely reproduce the observed impacts, whereas previous model versions fail to show any impact (0.3% decline, fig. S8). The differences in rainfall deficits between model versions are minor for wheat growing seasons, where both models reproduce observed yield anomalies (fig. S9).

The LPJmL model does not account for potential permanent physiological crop damages or mortality caused by adverse weather conditions. While plant growth is reduced and carbon assimilation is altered under water stress, grain filling continues when conditions allow. This model simplification has the potential to cause fundamental differences in simulated crop yields for different growing seasons (see Fig. 3 for the example of Germany).

In summer 2003, an anticyclone anchored over western Europe in June and July and held back the arrival of depressions, leading to low precipitation, soil moisture depletion, and above-average temperatures, which, in turn, accelerated heat-sum accumulation and thus crop development by about 20 days (39). As for Germany, LPJmL–PHU reproduces the observed shortening of the growing season, leading to a rainfall deficit of >1200 mm. Since LPJmL–Ref simulates a later (assuming later planting) and longer growing season (assuming excessive PHU requirements), the timing of maturity stretches into considerable rainfall events in October (Fig. 3, C and D)—when the harvest in LPJmL–PHU was completed already. In LPJml–Ref, late rainfall compensates early-season rainfall deficits and eventually leads to above-average yield levels in 2003, which is in sharp contrast to the observed decline. As a result of the refined weather forcing in LPJmL–PHU, simulated maize-yield anomalies are in close agreement not only with observations in Germany (R2 of 93 and 18%, respectively, for the two models; Fig. 3E) but also with the reported yield impacts across other affected countries during the 2003 heat wave in Europe (fig. S8).

## DISCUSSION

We derive location-specific PHUs from spatial variations in observed [MIRCA2000 crop calendar (27)] and simulated (LPJmL-derived; see fig. S5) sowing and harvest dates. These variations are used as proxy to implement geographic differences among crop cultivars into a global gridded crop model to deliberately match spatially explicit crop calendar dates with global crop model simulations. Constraining the growing season means improving the coincidence of weather forcing with phenological stages of crops, which substantially increases the agreement between reported and simulated crop yield fluctuations without actually changing the representation of weather responses in the model. The adjusted timing of growing seasons is particularly important for the quantitative reproduction of global average maize-yield responses to heat waves. Improved representation of irrigation system functioning and surface water constraints also improves the reproduction of interannual yield variability, yet to a lower degree compared to spatially resolved growing seasons. The updated model simulations reproduce average heat wave and drought effects on national yield levels, including higher sensitivities of maize than wheat yields to extreme heat (11, 14, 40), while providing new evidence that these effects are dominantly affecting rainfed farming.

Results indicate that observed yield fluctuations, especially in response to droughts and heat waves (which may as well imply agrometeorological droughts), are induced primarily by water stress. This is in line with in-depth analyses of observed effects of heat on U.S. maize yields (41, 42), which suggest that yield declines occurring around 30°C are less related to direct heat stress on reproductive organs than to temperature-related water stress due to increases in water demand and reductions in water supply. The sensitivity of photosynthesis rates to water stress is generally well represented in global gridded crop models that are designed to capture temperature effects on soil water depletion and increased atmospheric water demand (40, 41), as well as differences in water-use efficiencies in C3 and C4 crops due to differences in the CO2 fixation pathway (43). Water stress is also simulated to have lasting effects on the development of leaf area index, harvest index, and thus grain yield in the models.

Extreme temperatures also have the potential to directly damage enzymes, tissues, and reproductive organs and therefore hamper plant growth, especially near anthesis and until early grain filling—mostly in tropical and subtropical regions, where average maximum temperatures during the growing period are exceeding 30°C (40, 44). However, these direct negative effects on grain number and on the duration of grain filling (45) are rarely implemented in global gridded crop models (40, 46). While photosynthetic activity is down-regulated at temperatures exceeding the crop-specific optimal range and warmer seasons lead to faster phenological development and thus early senescence (see Materials and Methods), the crop model considered in our study neither accounts for permanent yield-reducing effects of heat nor does it represent lethal effects of extreme temperatures (33). These direct heat stressors are expected to be partially mediated by ample soil moisture and irrigation through reduced leaf temperatures, a mechanism also not explicitly implemented in the model. Without accounting for irreversible damages and crop mortality, photosynthesis rates and biomass accumulation are reduced within the model during periods of water and heat stress, while unimpeded plant growth continues as soon as conditions permit (see, for example, LPJmL–Ref in Fig. 3C). Despite the lack of irreparable heat damage in current global process-based crop models, they prove capable of reproducing observed county-level yield declines at extreme temperatures (40), as well as heat wave impacts on national yield levels (this study). This contributes to a debate on the relevance of heat stress in crop modeling at larger spatial scales (40, 41, 45, 47, 48) and contradicts recent findings from empirical modeling, suggesting that heat stress is often the most important predictor (23). Other empirical findings suggest that high temperature might affect global yields only when concurring with drought conditions (49).

Relative improvements in explained variance of yield simulations achieved in this study are consistent across most main producer countries (fig. S3), but low explanatory power remains in some countries such as Mexico, Indonesia, or the Philippines for maize (table S1) and Mexico, Kazakhstan, or Denmark for wheat (table S2). These countries are also located at the lower end of explained variances provided by ensembles of different crop models (24, 25). Yield fluctuations in general do not show a statistically significant weather dependency in about a third of all main producer regions (20, 25). Additional reasons are difficult to differentiate and may include not only model shortcomings (for example, ozone, flood, pest, and diseases) and missing representation of agronomic decisions (for example, fertilizer use, crop rotation, planting density, and general socioeconomic condition), but also low-quality yield reporting data in some cases (20, 25).

In addition, MIRCA2000 sowing and harvest data quality is not consistent and considered less reliable in large countries, where data are reported at national level only (for example, Russia, Spain, France, Ukraine, Nigeria, Pakistan, and Argentina) (27). In these countries, the LPJmL-internally calculated crop calendar often outperforms MIRCA2000 information (fig. S5). Crop calendar improvements at high spatial resolution, such as recently published for rice (50), are expected to further improve crop yield simulations—especially in regions where growing seasons are closely linked to the onset of the rainy season (51, 52). There is no global information about temporal variations in growing season, but recent studies show that the implementation of high-quality crop development dates reported annually for the United States even allows for a historical reproduction of observed temporal trends in U.S. yields (11, 12).

The implementation of spatially resolved PHUs improves the simulation of maize yields more than that for wheat. This difference may be due not only to (i) a higher ab initio agreement between the observed and simulated wheat yield fluctuations already achieved by the reference model, associated with a different heat unit approach compared to maize (see Materials and Methods), and (ii) the fact that wheat simulations include winter crops, which are generally harvested earlier and thus less vulnerable to summer extreme heat, but also to (iii) the fact that MIRCA2000 growing season information is known as being more robust for maize than for wheat (27). In only 14 of 26 main wheat producers (59% of the cropland area), MIRCA2000 crop calendar information leads to superior results compared to LPJmL–internally derived seasons (fig. S5).

Our analysis exclusively addresses weather effects on crop yields, but not the effects on harvested areas. Management decisions related to the extent of the harvested area can compensate or compound yield impacts from extreme events (14). Yet, more and, most of all, spatially explicit data are needed to better understand the relation of interannual climate variability and sowing dates and harvested area for specific crops (53). While these data would further benefit crop simulation if available globally, this study provides evidence that the multiyear average growing season timing at the grid-cell level already provides one critical tech-factor needed to reproduce annual national yield variability.

Better knowledge about the current relationship between local climate conditions and heat unit requirements may also be used to project future cultivar selection [for example, faster growing cultivars (28)] in response to regional climate change. Thus, by providing spatially explicit heat unit requirements, our study provides a basis for investigating the potential of management adaptation strategies in view of adverse climate change.

This is the first study to our knowledge to highlight that a spatially explicit implementation of cultivar phenology allows for an adequate reproduction of extreme weather impacts on the variability of global maize and wheat yields and that better agreement between simulated and observed crop yield variations may depend less critically on advanced model responses to weather forcing than on the correct representation of exposure to weather forcing. In view of the advanced model’s capacity to reproduce observed historical interannual yield variability and extreme event impacts, we highlight its role for refined impact analyses of future climate-change projections. It may lead to more reliable crop yield simulations and quantitative understanding of adaptation potentials.

## MATERIALS AND METHODS

### Experimental design

The agrohydrological model LPJmL was used to show how (i) better representation of water constraints in irrigation and (ii) the implementation of spatially resolved crop cultivars based on observed and modeled growing seasons contribute to GGCM’s ability to explain observed interannual yield variability at the global scale (Table 1). An additional model run eliminating crop-water limitations in both currently rainfed and irrigated cropland areas was used to quantify the underlying contributions of water stress to observed crop yield variabilities.

By improving on a reference model version that does not account for water constraints in irrigation and based on a semistatic heat-unit parameterization not constrained by observations [LPJmL–Ref (33); details below], we first introduce an advanced representation of irrigation system functioning and water constraints (38) and then adjust PHU requirements to match spatially explicit information on sowing and harvesting dates.

Model versions were evaluated on the basis of national standardized yield anomalies (1980–2009) reported by the United Nations’ FAO (8) using the Pearson correlation (R2) and RMSE. In addition, we tested the ability to reproduce global average effects of heat waves and droughts on national crop yields between 1964 and 2007 (compositing) (14), where extreme events were derived from the International Disaster Database EM-DAT (13) (see the “Statistical analysis” section for more details).

Model simulations were forced by the three bias-corrected reanalysis climate datasets GSWP3 (34), PGFv2.1 (35), and WATCH + WFDEI (36). If not stated otherwise, results refer to the average of three crop model simulations forced by each climate dataset, respectively (fig. S6).

### Sowing and harvest dates

Sowing dates and harvest dates in all model versions except for LPJmL–PHU were internally derived per crop on the basis of local climatic conditions (37). In contrast, LPJmL–PHU was additionally constrained by reported information about growing seasons under rainfed and irrigated conditions around the year 2000 at 0.5° spatial and monthly temporal resolution according to MIRCA2000 (27). The first day of a reported planting month was used as the planting date and the last day of a harvest month as the harvest date. MIRCA2000 differentiates between up to five different seasons per grid cell and crop. Grain and silage maize were distinguished neither in MIRCA2000 nor in LPJmL. MIRCA2000 was originally released at the 5–arc min spatial resolution, aggregated to 0.5° for the use in global crop models. However, the underlying information about growing seasons has a much lower spatial resolution with 402 units, where, for example Russia, Spain, France, Ukraine, Nigeria, Mexico, and Argentina are single spatial units without distinguishing subnational regions. An alternative global crop calendar (52) has been published, but it relies on largely similar data sources and does not distinguish rainfed and irrigated crops and was therefore not considered in this study. Given the limited spatial resolution of reported data, we complemented the MIRCA2000 information with LPJmL–internally derived growing seasons based on local climatic conditions (37). For each therefore available growing season, we derived grid cell–specific heat units required to match targeted harvest dates (fig. S4). Per country, we selected the season that leads to the highest correlation between simulated and observed crop-yield anomalies when used as input in LPJmL–PHU. Figure S4 shows individual model performances of all seasons, for simulations solely based on LPJmL-derived growing seasons but with spatially resolved heat units, and also for simulations with growing seasons exclusively based on MIRCA2000 information. For final LPJmL–PHU simulations as described in Table 1, MIRCA2000 information was used across 80 and 59% of total main producers’ cropland area for maize and wheat, respectively (see the map in fig. S5).

### PHU requirements and phenological development

The simulation of crop phenological development from emergence through anthesis to maturity was controlled by temperature and modified by vernalization requirements. Spatially explicit accumulated thermal units—here PHU (in °C day)—are based on daily weather data averaged for the 1980–2010 period, individually derived from the three observational climate datasets mentioned above. For each climate dataset, crop, growing season (see above), and grid cell (i), we determined heat units required to reach maturity (PHUreq,i) by calculating the sum of daily (j) average temperature above a crop- and location-specific base temperature (Tb,i), weighted by a vernalization factor Vi,j from sowing (j = 0) to the harvest day (J)(1)where base temperatures Tb,i are calculated as(2)with minimum temperatures = 5°C for maize and 0°C for wheat, maximum base temperatures = 15°C for maize and 0°C for wheat (33, 34), and annual average temperatures .

For maize, we assume no vernalization requirement (V = 1). Winter wheat is assumed to require exposure to low temperatures to reach anthesis (54, 55). It is distinguished from spring wheat if sowing takes place between the 243rd and 365th day of the calendar year on the Northern Hemisphere and between the 59th and 181st day of the calendar year on the Southern Hemisphere and if the mean temperature of the three coldest months within the growing season is below 10°C. Heat unit accumulation is down-regulated by V < 1 as long as the vernalization requirement Vreq is not yet fulfilled. Vreq is derived from cell-specific average winter temperatures of the five coldest months within the growing season.

For winter wheat, we assume a maximum vernalization requirement of 70 days. That being said, fewer vernalization days are required if monthly average temperatures of the coldest months i1,..., i5 are above 3°C(3)where if monthly average temperature is below 3°C and linearly decreases to zero at a monthly temperature of 10°C.

To reach the vernalization requirement, each day in the growing season is weighted according to their vernalization effectiveness and successively added up. Here, we assume that vernalization effectiveness is highest (=1) if daily temperature is between 3° and 10°C. Below 3°C, it linearly decreases and reaches 0 at −4°C, and above 10°C, it also decreases linearly to 0 at 17°C (55). As soon as Vreq is reached, the vernalization factor for daily heat sum accumulation V is set to 1. Before a minimum vernalization requirement of Vreq/5 is reached, V is set to 0. In between, it increases linearly from 0 to 1 with accumulated weighted vernalization days (29).

In the default phenology model (referred to as semistatic), LPJmL–Ref (and also LPJmL–WaterLimIrr and LPJmL–NoWaterStress), heat requirements are not derived from reported growing seasons, and sowing dates are internally derived on the basis of climate conditions. For maize, the harvest date is not based on local observations but on a global constant , while only the base temperature (Tb,i) is assumed to vary geographically (same as Eq. 2). For spring wheat, heat requirements are defined as a multiple of the mean annual temperature (, 1980–2010 mean) or a maximum value of and a minimum value of . For winter wheat, is calculated as(4)where Ds is the sowing date and Dk is 365 (Northern Hemisphere) or 181 (Southern Hemisphere). The default model also uses a simplified approach to the effect of vernalization with geographically constant vernalization requirements.

Phenological development is a function of daily accumulation of weighted heat sum increments over the growing season (PHUsum,j, calculated as in Eq. 1). Anthesis is assumed to occur as soon as PHUsum,j reaches 0.5*PHUreq for maize and 0.45*PHUreq for wheat (33). Physiological maturity is assumed to be reached as soon as PHUsum,j reaches either PHUreq or a crop-specific age limit of 240 days for maize, 334 days for spring wheat, and 364 days for winter wheat, which represents the maximum growing season length in the MIRCA2000 crop calendar. We further assume time requirements between crop maturity (green leaf area index reaches zero) and harvest, as suggested by Elliott et al. (30), that is, 21 days for maize and 7 days for wheat. As phenological development of specific crops can also be influenced by day length, PHU accumulation in LPJmL can be weighted by photoperiodism, as described in Schaphoff et al. (33), but results were outperformed by the simpler model without photosensitivity, which was thus presented herein. Additional physiological stresses were not considered to directly affect phenological development but explicitly influence the simulation of photosynthesis and carbon allocation (see next section).

### The global gridded crop model LPJmL

LPJmL mechanistically represents biogeochemical land-surface processes at the global scale, simulating daily water fluxes in direct coupling with the establishment, growth, and productivity of major natural and agricultural plant types at 0.5° resolution (33). Spatially explicit data on cropland extent were obtained from the MIRCA2000 land-use dataset (27), the extent of areas equipped for irrigation from Siebert et al. (56), and the distribution of irrigation systems (surface, sprinkler, and drip—applies only to simulation LPJmL–WaterLimIrr and LPJmL–PHU; see Table 1) from Jägermeyr et al. (38). Land-use patterns were held constant at year 2005 levels, and sowing dates were fixed throughout the simulation period. Crop yields were calibrated to match average (1998–2002) reported national management intensities (8).

Photosynthesis is simulated in dynamic coupling to absorbed photosynthetically active radiation, temperature (optimum for photosynthesis, maize 21° to 26°C and wheat 12° to 17°C), day length, canopy conductance, and water stress (33). Daily carbon assimilated through photosynthesis is allocated to the crop organs roots, leaves, harvestable storage organ (for example, cereal grain), and stem. Allocation to each compartment is a function of the phenological development stage and biomass increment. LPJmL computes daily water stress as the ratio of soil water supply to the crop water demand. Increased vapor pressure deficit or water stress therefore affects the daily rate of carbon assimilation through reduced stomatal conductance. The daily increment of leaf area index is also down-regulated by water stress, with lasting effects over the growing season, affecting the harvest index and thus grain yield.

Model simulations followed a 900-year (no land use) and 120-year (land use) spin-up, recycling the first 20 years of input climatology (mentioned above). LPJmL–PHU outputs were corrected for unintended multiple harvests potentially occurring within a single calendar year when harvest dates oscillate around the end of the year. In cases where multiple harvest events are detected, yields were reported for the year in which most of the growing season occurred. The full irrigation setup (LPJmL–NoWaterStress) assumes unrestricted access to irrigation water to fulfill crops’ water demand (any soil water deficit is balanced the next day).

### Statistical analysis

Explained variances (R2, in percent) are based on Pearson correlation coefficients derived from simulated and observed national time series of standardized crop yield anomalies (1980–2010), calculated as detrended (first quadratic polynom subtracted) and normalized (mean subtracted) yields, divided by the SD (hereinafter “yield anomalies”). To quantify residuals between the observed and simulated yield anomalies, we calculated RMSEs (unitless, as standardized anomalies are without unit). Note that the Pearson correlation coefficient is unaffected by the standardization of anomalies.

The composite analysis (Fig. 2 and fig. S8) was constructed by extracting a 7-year time window from historical national yield time series, centered on the respective “extreme year.” Extreme years were defined according to EM-DAT (13), which reports the country and year for various extreme event types from 1961 to 2010, if at least 10 people died, 100 or more people were injured, made homeless, or required immediate assistance, or a country declared a state of emergency, or called for international assistance. In this study, we considered all heat waves and droughts recorded between 1964 and 2007 (3 years before and after an extreme event required for the construction of the composite analysis), if the respective crop contributes more than 5% to the total cropland area in the affected country [based on Portmann et al. (27)]. For multiyear extreme events, we averaged consecutive extreme years into a single disaster year signal, so that the time window always consisted of seven entries, centered on the event signal.

This method creates a subset of 65 heat waves and 175 droughts for maize and 81 heat waves and 146 droughts for wheat (table S3). The extracted 7-year time series were divided by the average of the 3 years preceding and following the event to remove the absolute magnitude of national data from the signal. Any other data entry co-occurring with another extreme event was excluded from calculating the mean. Last, any linear trend was removed from the composite time series [in contrast to the study by Lesk et al. (14)], as model simulations are not expected to reproduce observed trends in yields driven by technological progress. The detrending was applied after compositing as it maintains the absolute event signal, which is thus directly comparable to previous studies (14). A discussion on potential type I and II errors associated with compositing sample size was presented in Lesk et al. (14).

Rainfall deficits in Fig. 3 were calculated as the relative difference of cumulative growing season precipitation in the respective year and the long-term average cumulative precipitation during the MIRCA2000 growing season.

## SUPPLEMENTARY MATERIALS

Fig. S1. Root mean square error of standardized country-level yield anomalies for maize and wheat.

Fig. S2. Observed and simulated historical maize and wheat yield anomaly time series.

Fig. S3. Sensitivity of mean explained variance to number of countries considered.

Fig. S4. Evaluation of available growing season inputs.

Fig. S5. Best-performing crop calendar per country.

Fig. S6. Evaluation of different climate inputs.

Fig. S7. Influences of heat waves and droughts on rainfed and irrigated yields.

Fig. S8. Observed and simulated influences of the 2003 European heat wave on maize yields.

Fig. S9. Effects of growing season adjustment on wheat rainfall deficit.

Table S1. Explained variances and RMSE of maize country-level yield anomalies.

Table S2. Explained variances and RMSE of wheat country-level yield anomalies.

Table S3. List of extreme events considered in this study.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

Acknowledgments: Funding: This study was funded within the framework of the Leibniz Competition (SAW-2013-PIK-5) and with partial support from the U. Chicago Center for Robust Decision-making on Climate and Energy Policy (NSF grant #SES-146364). We thank C. Müller and S. Schaphoff for helpful comments on the study design and two anonymous reviewers for their excellent comments that helped improve the manuscript. Author contributions: J.J. and K.F. designed the study. J.J. developed the model code and performed the model simulations and analysis. J.J. and K.F. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper and the model code may be requested from the corresponding author.
View Abstract