# Warming reduces global agricultural production by decreasing cropping frequency and yields

Peng Zhu, Jennifer Burney, Jinfeng Chang, Zhenong Jin, Nathaniel Mueller, Qinchuan Xin, Jialu Xu, Le Yu, David Makowski, Philippe Ciais

► To cite this version:

Peng Zhu, Jennifer Burney, Jinfeng Chang, Zhenong Jin, Nathaniel Mueller, et al.. Warming reduces global agricultural production by decreasing cropping frequency and yields. Nature Climate Change, 2022, 12 (11), pp.1016. 10.1038/s41558-022-01492-5 . hal-03810580

HAL Id: hal-03810580

https://hal.science/hal-03810580v1

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Warming reduces global agricultural production by decreasing cropping frequency and yields

Received: 4 September 2021

Accepted: 5 September 2022

Published online: 10 October 2022

<figure><img src="image_2.png" /></figure>

Check for updates

Peng Zhu Ⓓ1,2,3✉, Jennifer Burney Ⓓ4, Jinfeng Chang Ⓓ5, Zhenong Jin Ⓓ6, Nathaniel D. Mueller7,8, Qinchuan Xin Ⓓ3, Jialu Xu Ⓓ9, Le Yu Ⓓ10,11, David Makowski Ⓓ12 and Philippe Ciais Ⓓ1

Annual food caloric production is the product of caloric yield, cropping frequency (CF, number of production seasons per year) and cropland area. Existing studies have largely focused on crop yield, whereas how CF responds to climate change remains poorly understood. Here, we evaluate the global climate sensitivity of caloric yields and CF at national scale. We find a robust negative association between warming and both caloric yield and CF. By the 2050s, projected CF increases in cold regions are offset by larger decreases in warm regions, resulting in a net global CF reduction (−4.2 ± 2.5% in high emission scenario), suggesting that climate-driven decline in CF will exacerbate crop production loss and not provide climate adaptation alone. Although irrigation is effective in offsetting the projected production loss, irrigation areas have to be expanded by >5% in warm regions to fully offset climate-induced production losses by the 2050s.

Global food demand is expected to increase in the coming decades with growing population and shifting dietary patterns<fnref n="1" />,². In the past, global production has kept pace with rising demand through both cropland expansion and intensification. Cropland expansion into uncultivated areas is an expedient but unsustainable way of increasing crop production, since it has caused cascading environmental harms like increased soil erosion, loss of wildlife habitat and substantial carbon emissions³,⁴. Crop intensification includes within-season yield improvements through use of high-yielding varieties, fertilizer, pesticide and other inputs²,⁵, as well as increased cropping frequency (CF) or the number of crops cultivated per year on a given area of land. At the global scale, increases in caloric production (CalP) have been

driven mainly by improvements in crop caloric yield (CalY, crop yield in calories produced per area), whereas increases in CF and cropland expansion have played relatively minor roles (Fig. 1)⁶,⁷.

Anthropogenic climate change is expected to make efforts to increase total crop production more difficult due to rising average temperatures and more frequent extreme weather events⁷, directly threatening food security objectives codified in the United Nations Sustainable Development Goals. It has been noted that each component of crop production might be influenced by climate change through different pathways and over different time scales⁹,¹⁰. Crop yield is generally sensitive to in-season and interannual climate variability, and there is robust evidence that, without major adaptation, global staple crop

<footnote n="1">¹Laboratoire des Sciences du Climat et de l'Environnement, CEA CNRS UVSQ Orme des Merisiers, Gif-sur-Yvette, France. ²Department of Geography, The University of Hong Kong, Hong Kong SAR, China. ³Guangdong Key Laboratory for Urbanization and Geo-simulation, School of Geography and Planning, Sun Yat-Sen University, Guangzhou, China. ⁴School of Global Policy and Strategy, University of California, San Diego, La Jolla, CA, USA. ⁵College of Environmental and Resource Sciences, Zhejiang University, Hangzhou, China. ⁶Department of Bioproducts and Biosystems Engineering, University of Minnesota, St. Paul, MN, USA. ⁷Department of Ecosystem Science and Sustainability, Colorado State University, Fort Collins, CO, USA. ⁸Department of Soil and Crop Sciences, Colorado State University, Fort Collins, CO, USA. ⁹School of National Safety and Emergency Management, Beijing Normal University, Zhuhai, China. ¹⁰Department of Earth System Science, Ministry of Education Key Laboratory for Earth System Modeling, Institute for Global Change Studies, Tsinghua University, Beijing, China. ¹¹Ministry of Education Ecological Field Station for East Asian Migratory Birds, Beijing, China. ¹²Université Paris-Saclay, INRAE, AgroParisTech, UMR 518 MIA, Palaiseau, France. ✉e-mail: peng.zhu@lsce.ipsl.fr</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Changes in crop CalP, CF, CalY and cropland area during 1979–2018. Here, different regions are defined by FAO. Each time series is normalized using data points at year 1979. Note different y axis scale for South America.</figcaption></figure>

yield will be reduced under climate change<sup>11,12</sup>. Climate change further influences cropland area (CA) over longer time scales, with cropped areas migrating slowly to counteract environmental changes<sup>13</sup>. CF is potentially susceptible to climate changes over multiple time scales. Dry or excessively wet soil conditions before planting schedule might delay or prevent a second cropping season<sup>14,15</sup>. Warmer temperature might be beneficial for increasing CF in cold areas, as the frost-free period expands and shortens single-crop growth duration<sup>16–18</sup>. However, rising temperatures in warm regions, which have an already long growing season, might result in lower CF due to the higher risk of crop failure (total loss of cropping cycle) caused by more frequent heat or drought stresses<sup>19,20</sup>. These dynamics have been confirmed for the Brazilian corn-soybean cropping system using satellite-observed CF, suggesting that CF will be reduced by warming to a greater degree than will yields<sup>19</sup>. These divergent conclusions support the need for holistic quantification of how climate warming influences global-scale CF to manage climate warming impact on regional and global food security.

Studies reporting positive effects of warming on CF normally relied on prescribed crop phenology models driven by temperature or 'space-for-time' substitution to infer the potential change<sup>16,17,21,22</sup>. Such modelling strategies may overlook other environmental factors constraining the feasibility of boosting CF<sup>23,24</sup>, such as increased heat stress in warmer conditions<sup>8</sup>, decreases in water availability for irrigation<sup>25</sup>, drought stress<sup>26</sup>, soil fertility decreases and increased pests and pathogens in warmer climates<sup>27</sup>. Recent progress on using satellite data or statistics-derived regional-scale CF makes it possible to characterize regional disparities and explore how climate variation will influence CF with a data-driven approach<sup>19,28,29</sup>. However, most data-driven model studies have focused on crop yield<sup>12,30,31</sup>, harvested area<sup>32</sup> and total factor productivity<sup>33</sup>; such models have not been used to assess climate impacts on CF specifically.

Here, we address this gap by building an empirical model fed by the Food and Agriculture Organization (FAO) national statistical data (http://www.fao.org/faostat/en/#data/QC) to examine how and to what degree recent climate changes have altered CF and to assess the potential impact of future warmer climate on CF and total food production globally. Annual national mean CF is obtained by summing the harvested area of all crops produced within a country (up to 161

monitored products) and then dividing it by the national total CA–providing an estimation of the frequency cropland resources are used<sup>21</sup> (Supplementary Fig. 1). It is an aggregate of all cultivation patterns, including fallow croplands (CF = 0) used for conservation purposes or resulting from land abandonment, areas with single-cropping per year (CF = 1) and areas with multiple-cropping cycles per year (CF > 1). We also calculate the national annual total CalP by summing the caloric content of all national crop production and further obtain the national CalY by dividing CalP by the total harvested area of 161 crops in a nation (Methods). On the basis of these definitions, CalP = CalY × CF × CA. Although crop production additionally supplies critical protein, vitamins and fibre, calories from food crops are essential for meeting dietary energy requirements, and thus, calorie-weighted metrics are highly relevant for informing global food security. Taken together, evaluating the climate sensitivities of these four metrics provides a global picture of crop production change and adaptation potential under future climate warming.

## Temporal changes in production

Over the study period (1979–2018), global CalP increased by 125% and CF, CalY and CA increased by 18%, 73% and 11%, respectively (Fig. 1). There are substantial regional variations in growth; South America shows the greatest CalP increase (by 320%), followed by Africa (by 220%) and Asia (by 207%) (Fig. 1). However, unlike most continents, where CalY is the leading driver of CalP, the increase in African CalP is primarily driven by CA over time, corroborating other research<sup>1,2</sup> suggesting that there are large potentials to improve African crop production with yield intensification. In fact, CF is the second primary driver of increasing CalP in several regions, including North America, South America and Asia. By contrast, CA explains only a small share of the CalP increase in most regions. After 2010, the rate of CA increase has slowed down in all continents and CA even started to decrease in North America and Western Europe (probably due to land resource conservation policies<sup>6</sup>) and in Eastern Europe.

Although most continents show increasing trends in CF, CalY and CalP, the time series of these variables are characterized by strong year-to-year fluctuations, especially in Oceania, where the interannual variability of CalP seems primarily driven by the concordant temporal
<figure><img src="image_4.png" /></figure>

**Fig. 2 |** Response of CF, CalY and CalP to temperature. **a–c**, Response function of CF (**a**), CalY (**b**) and CalP (**c**) to Tmean with and without irrigation. Response functions are established on the basis of bootstrap. The 1,000 estimated regression coefficients are then used to determine the 95% CI (shaded areas) of model-estimated coefficients. The medians of the coefficients are used to determine the response curves. Two response curves corresponding to without and with irrigation are determined through setting irrigation fraction as zero and historical global mean, respectively. The 2.5th and 97.5th percentiles are used to define the lower and upper bounds of the 95% CI. Here, the curves are shifted

vertically so that the peak values of CF, CalY and CalP under no irrigation are 1. Triangles on the x axis indicate the optimal temperature with different levels of irrigation. **d**, Marginal effect of 1 °C warming on global average CF, CalY and CalP estimated with eight panel models (the numbers under the line correspond to the models in Methods). The global average of marginal effect of 1 °C warming is a country crop area weighted average of warming effects in each country. Error bars represent 95% CI of each estimation. The ensemble mean of eight panel models estimation is indicated by the horizontal line with shaded area as the 95% CI of ensemble mean.

variation in CalY (Fig. 1). Detrending CalY and CF with a spline function shows that the interannual variations of CF and CalY are positively correlated at the global scale and in several representative countries or regions (Supplementary Fig. 4). This synchronization suggests that both CF and CalY might be sensitive to climate anomalies. For example, adverse environmental conditions driving declined crop yield might also lower farmer's prospects for raising a second crop in the same year¹⁹,³⁴; extreme drought might result in a total crop failure¹⁴. However, we also note that CF and CalY in China and Brazil are not significantly (*p* > 0.05) correlated, perhaps because irrigation and other socioeconomic factors have decoupled multiple-cropping decisions from crop productivity there.

## Climate effects on CF, CalY, CA and CalP

We estimate best-fit parameters for panel regression models (baseline model M1; Methods) relating national CF, CalY or CalP to climate (annual mean temperature and precipitation) and major management practices (irrigation and nitrogen fertilizer application). With the median coefficients *β*, two response curves corresponding to 'without irrigation' and 'with irrigation' can be determined through setting irrigation fraction as zero and historical global mean (3.1%), respectively (Fig. 2a–c). Without irrigation, the estimated temperature response functions of CF, CalY and CalP peak at 7.5, 2.6 and 4.2 °C, respectively (Fig. 2a–c). As we use annual mean temperature rather than growing

season mean temperature, the estimated optimal temperature of CalY is lower than the optimal temperature of the crop yield recently estimated³¹. Additionally, due to significant interactions (*p* < 0.05) between temperature and irrigation area fraction (Supplementary Tables 1–3), more irrigation not only increases CF, CalY and CalP but also raises the optimal temperature (Fig. 2a–c), suggesting the heat stress mitigation potential of irrigation application. Similarly, the nonlinear responses of CF, CalY and CalP to annual precipitation suggest that higher precipitation promotes CF, CalY and CalP until a high threshold is reached (Supplementary Fig. 5). The significant negative interaction (*p* < 0.05) between irrigation area fraction and precipitation suggests that irrigation does indeed ease water availability constraints on CF, CalY and CalP (Supplementary Tables 1–3).

With this regression model, we also quantify the notable beneficial effect of fertilizer application (Supplementary Tables 1–3). Our empirical model suggests that there will be 15.4%, 6.64% and 7.61% increases in CF, CalY and CalP per 100 kg ha⁻¹ application of nitrogen fertilizer, respectively. In contrast, more fertilizer application is predicted to reduce CA (Supplementary Table 4), as higher crop yield reduces the need for land reclamation and slow cropland expansion³⁵.

Using the coefficients from baseline model M1, we show—under historical management practices—that +1 °C increase in annual mean temperature will reduce global average CalY, CF and CalP by 3.6 ± 1.5%, 1.8 ± 0.9% and 5.7 ± 1.9% (mean ± 95% confidence interval (CI)) (Fig. 2d),
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 3 | Effect of warming and irrigation on CF, CalY and CalP. a-c,** Marginal effect of 1 °C warming on CF (a), CalY (b) and CalP (c) for each country estimated with baseline model. The marginal effect of 1 °C warming was estimated as the difference between predicted CF, CalY or CalP with uniform 1 °C warming and the original CF, CalY or CalP, based on model M1 and historical irrigation fraction.
**d-f,** Sensitivity of CF (d), CalY (e) and CalP (f) to irrigation area fraction, estimated with the baseline model M1 ( $(\frac{\partial Y_{i,t}}{\partial \text{irrig}_{i,t}} = \beta_3 T \text{mean}_{i,t} + \beta_6 \text{Prcp}_{i,t})$ .</figcaption></figure>

Sensitivity for each country can be obtained with multiyear mean climate variables during 1979–2018. As $Y_{i,t}$ is logged, the estimated sensitivities indicate the percentage change in CF (CalY or CalP) with unit percentage change in irrigation area fraction. The insets in **d-f** show the sensitivity of CF, CalY and CalP to irrigation fraction in climate space, which is delineated by country-level annual precipitation (x axis) and annual mean temperature (y axis) during 1979–2018. Colours in each climate space indicate the sensitivity of CF, CalY and CalP to irrigation fraction.

respectively. This estimated value of CalY decline is quantitatively close to the estimation of 1 °C warming effects on staple food crops (maize, rice, soybean and wheat) on the basis of a similar statistical model[12,30,31], probably because staple food crops constitute a large proportion of global calorie production. Further, the reduced CF could be because higher temperatures increase heat extremes resulting in more crop failure, as well as reducing available soil water supply[8,36], making it hard to support multiple-cropping seasons. Warming effects on CalP are close to the sum of individual warming effects on CalY and CF, whereas the complementary components—warming effects on CA—are very small and not statistically significantly (p > 0.05) different from zero (Supplementary Fig. 6). This is also reflected by the statistically non-significant (p > 0.05) regression coefficients relating climate variables to CA (Supplementary Table 4).

Models based on spline function of climate variables (M2), quadratic function of 2-month climate variables with LASSO regression

(M3) and a growing degree day (GDD) model with LASSO regression (M4) produce similar estimations as the baseline model (Fig. 2d). In addition, we find the baseline model estimation still holds when using other climate datasets (M5 and M6), an alternative time period to aggregate climate variables (M7) and a comprehensive management index (total factor productivity input) to represent agricultural management practices (M8).

With the nonlinear temperature response functions driven by historical observations, the baseline model shows prevalent decline in CF, CalY and CalP for warm areas (for example, Africa) but lower declines or even small gains for cold areas per 1 °C of warming (Fig. 3a-c). For example, CF in Canada, Scandinavian countries, Mongolia and Russia increases with +1 °C warming (Fig. 3a-c). These increases, however, are insufficient to fully offset the declines in warm (and dry) regions, resulting in an overall decline in global average CF (Fig. 2d). Such spatial divergences could explain inconsistent findings of CF response
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

**Fig. 4 | Projected changes in CF, CalY and CalP for 2031–2070 relative to the reference period 1979–2018.** a,b, Boxplot of global and continental changes in CF, CalY and CalP under SSP 126 and SSP 585 estimated with baseline model. Boxplots indicate the median (horizontal line), 25–75th percentile (box) and

5–95th percentile (whiskers) of estimated change of all year and country combinations relative to historical period (1979–2018). **c–h**, Projected changes in country-level CF (**c,d**), CalY (**e,f**) and CalP (**g,h**) under SSP 126 and SSP 585 for 2031–2070 relative to the reference period 1979–2018.

to climate warming in different geographical areas¹⁷,¹⁹,³⁷. Overall, the estimated +1 °C warming effects on CF partially offset the declines in CalY found in Canada, Scandinavia and Russia but contribute to greater declines in CalP in warm areas like South America and Africa.

## Benefits of irrigation on CF, CalY and CalP

The temperature and precipitation response curves in different levels of irrigation (Fig. 2 and Supplementary Fig. 5) raise the question of how much irrigation can partially offset the negative effects of warming and drought stress on CF, CalY and CalP. With the interaction term between irrigation area fraction and climatic variables in the baseline model, we

are able to assess the sensitivity of CF, CalY and CalP to irrigation area fraction, which is expressed as a linear function of temperature and precipitation. We find a positive and statistically significant (*p* < 0.05) sensitivity of CF, CalY and CalP to irrigation area fraction with higher beneficial effect in warm and dry conditions (Fig. 3d–f). For example, in tropical and subtropical regions, the sensitivity of CalP to irrigation can reach 5%, meaning that a 1% increase in irrigation area fraction will lead to 5% increase in CalP. As these regions are characterized by a high decline of CF, CalY and CalP with temperature rise, irrigation expansion could be an effective way of mitigating production losses in this part of the world, provided that there are sufficient water supplies.
<figure><img src="image_19.png" /><figcaption>**Fig. 5** | Projected irrigation area fraction to offset climate change-induced decline in CalP. a,b, The irrigation area fraction required to offset climate change-induced decline in CalP is projected with baseline model M1 driven by ensemble mean of future (2031–2071) climate models outputs under SSP 126 (a) and SSP 585 (b) and historical reference period (1979–2018) climate dataset.</figcaption></figure>

## Changes in CF, CalY and CalP for 2031–2070

Projected changes in CF, CalY and CalP for 2031–2070, relative to the reference period 1979–2018, show in the low-warming scenario Shared Socioeconomic Pathway (SSP) 126, the projected global average changes in CF, CalY and CalP are equal to –2.9 ± 0.5%, –5.8 ± 1.1% and –9.0 ± 1.7% (mean ± 95% CI), respectively (Fig. 4a). In the higher warming scenario SSP 585, climate change leads to greater reductions in CF (–4.2 ± 0.6%), CalY (–8.1 ± 1.3%) and CalP (–12.7 ± 2.1%) (Fig. 4b). The projected impacts show substantial regional disparities (Fig. 4c–h; see Supplementary Figs. 7 and 8 for separate effects of temperature and precipitation). For example, CF is projected to increase in both scenarios in Canada, Scandinavian countries, Mongolia and Russia (representing –12% of the world cropland), while CF would decline in other countries. The United States, Western Europe and China show slight declines in CF, CalY and CalP in both scenarios, while stronger declines are identified in tropical countries. Specifically, in Brazil and Sub-Saharan African countries, CalP is projected to decline up to 30% under SSP 585. In these tropical countries, the projection suggests that without additional adaptation, future warming will result in reduced CF that will exacerbate the total CalP loss and threaten food supply.

On the basis of the estimated positive impacts of irrigation, we make a zero-order projection of how much additional irrigation expansion would be required for each country to fully offset the anticipated climate-induced decline in total CalP. We find that greater increases (increase to 5% irrigation fraction in SSP 585 scenario) are required for countries in Africa and South America than in other countries (Fig. 5). This is because a greater decline in CalP is expected in these warm countries due to climate change (Fig. 4g,h) and also because these countries are characterized by lower irrigation fractions in the historical period (Supplementary Fig. 2a). In fact, our projection suggests that irrigation fractions in countries like Congo, Angola and Zambia would have to be boosted by more than 20 times relative to their very low historical baseline levels. However, the increasing rate of irrigation fraction in Africa during the historic period is <0.2% per decade (Supplementary Fig. 2b), suggesting that the historical increasing rate in irrigation fraction has to be augmented to satisfy the future irrigation requirement. Considering the expansion of irrigation fraction in Africa is not only limited by local water resource availability but also the infrastructure investment like water pumping system construction³⁸,³⁹, there are substantial challenges to achieving these targets to fully offset anticipated food production losses in these areas. However, future precipitation changes provide opportunities for irrigation expansion. For example, countries in Eastern Africa, South Asia and Southeastern Asia are projected to have precipitation increases (Supplementary Fig. 9), which might facilitate the building of irrigation infrastructure. However, actual irrigation water withdrawal will also depend on the competition for water resources from other sectors.

## Discussion and conclusion

It has been suggested that increasing multiple-cropping practices could be an effective adaptation strategy to future warmer climate¹⁶,¹⁷,³⁷, since future warming will expand the annual window for crop production in cold areas and potentially allow for higher CF to offset loss from lower yield. However, we show that warming is expected to reduce global average CF, especially in warm areas, possibly because heat and associated drought stresses disrupt the normal planting or harvesting time window¹⁴,⁴⁰ or the depleted water supply is unable to sustain multiple-cropping²⁵. Even for northern cold areas, our predicted CF increases are smaller than previous estimations¹⁶,³⁷. This discrepancy could be explained by multiple-cropping practices not only being limited by temperature but also other environmental and socioeconomic factors. For example, converting single-cropping to double-cropping will result in higher water demand⁴¹, which might not be satisfied by water supply when local infrastructures and irrigation techniques are insufficiently developed. Multiple-cropping might also require additional adjustments in farming systems⁴², which might not always be economically practicable. As our analysis is based on FAO statistics which reflect the actual multiple-cropping practices, the estimated values of CF reported here probably account indirectly for more of these environmental and socioeconomic constraints than do the technical potential studies¹⁶,¹⁷,³⁷.

While the projected CF increase in northern cold areas is relatively small, it would be possible to achieve a substantial increase in future agricultural production by extending cultivated land into areas that are currently too cold from an agronomic perspective⁴³,⁴⁴. If cropland expansion is taken as the only approach to mitigating the projected CalP loss, based on our projection (Fig. 4), CAs would have to increase by 9.0 ± 1.7% (SSP 126 or 0.18% per year) and 12.7 ± 2.1% (SSP 585 or 0.25% per year) by 2031–2070. These calculations assume that new CAs exhibit the same average yields of existing lands, although in reality new croplands often yield lower than average¹⁵. During 1979–2018, global CAs show an increasing rate of 0.2% per year (Supplementary Table 6); therefore, if cropland expansion holds this momentum, globally, the projected CalP loss will be fully offset in SSP 126 and partially offset in SSP 585. However, climate warming might change the spatial pattern of global agricultural suitability⁴⁶.

These findings are robust to different model specifications (M1–M4), climate data (M5 and M6) (Fig. 2d), alternative time period aggregation (M7), to the use of a comprehensive management index (M8) (Fig. 2d and Supplementary Fig. 10), different polynomial functions to characterize the trends (Supplementary Figs. 11–12) and a subset of high-quality data records (Supplementary Fig. 13). We also validate FAO country-level CF with CF derived from subnational statistics and satellite data (Supplementary Fig. 14) and conduct robustness checks by running the regression model with satellite data-derived CF. This robustness check shows that satellite data-derived CF has a similar
temperature response curve and estimated warming effects (Sup-
plementary Fig. 15).

We test whether our findings could be biased by countries with extreme high or low temperatures. Removing the 10% coldest/warmest countries did not change noticeably the estimated warming effects on CF, CalY and CalP (Supplementary Figs. 16 and 17), suggesting that our estimation is robust to samples with extreme high or low temperatures. Another potential concern is the effect of crop species omission in FAO statistics on our findings. We simulate this through randomly discard- ing a portion of crop species (Methods) and find that the estimated temperature effects are robust to different levels of crop species omission (Supplementary Fig. 18).

Several caveats apply here. It is critical to find the right balance between too-simple highly biased models neglecting important inputs and too-complex models including too many inputs with overfitting. Indeed, if some important inputs are neglected, the effect of a specific input on the response can be strongly over- or under-estimated. The selection of inputs to be included in a regression model is therefore a key step. First, we use 18-month mean climate variables (or climate exposures) rather than climate variables during specific crop growing seasons, since there is no explicit crop growing season information covering all 161 crops at global scale. In addition, non-growing season climate might be also relevant to model climate effects on CF, CalY and CalP, especially for CF, as planting decisions are largely determined by preseason weather conditions⁴⁷. Therefore, we took the climatic conditions in both current year and the second half of the year pre- ceding the harvest year into account in our models. It is worth noting that our LASSO statistical model (M3) addresses this issue, as it auto- matically selects the seasonal climate variables in certain time periods which explain most of the variation of CF, CalY and CalP. Second, our statistical model does not account for the fertilization effect of rising atmospheric CO₂, which might partially offset future warming-driven crop yield losses. Since biophysical models normally account for this, such omission might cause discrepancy between statistical models and biophysical models, especially in high-warming scenarios that accompany substantially elevated CO₂ concentrations. However, we note that because CO₂ fertilization effect is further complicated by envi- ronmental conditions (soil fertility, temperature or drought stress), a better understanding of the real potential of CO₂ fertilization to offset future crop losses requires more comprehensive synthesis of field experiments⁴⁸. Third, as already identified, rainy season statistics (for example, rainy season onset, duration and cessation) might influence CF through regulating the planting season soil water status¹⁴,¹⁵, but these factors are not accounted for here. This is mainly because climate models normally show a higher uncertain projection for precipitation relative to temperature⁴⁹, which makes deriving the rainy season statis- tics more challenging. Considering the rising precipitation variability and more extreme events³⁰,⁵¹ in the future, it is important to improve the climate model's predictive capacity for precipitation and derive more accurate rainy season statistics⁵².

We show that future warming is expected to have a positive but lower than expected effect on CF in northern cold areas and will prob- ably induce a decline of CF in warm areas like Africa and South America, eventually reducing global average CF. The projected decline in global CF represents about one-third of total projected CalP losses. Increased irrigation through increased irrigated fraction or increased duration could nevertheless be an effective approach compensating for the negative impacts of heat and water stress, but its large-scale implemen- tation might be constrained by infrastructure development, economic returns and local water supply in each country⁵³,⁵⁴.

Online content

Any methods, additional references, Nature Research reporting
summaries, source data, extended data, supplementary informa-
tion, acknowledgements, peer review information; details of

author contributions and competing interests; and statements of
data and code availability are available at https://doi.org/10.1038/
s41558-022-01492-5.

**References**

1. Tilman, D., Balzer, C., Hill, J. & Befort, B. L. Global food demand and the sustainable intensification of agriculture. *Proc. Natl Acad. Sci. USA* **108**, 20260–20264 (2011).

2. Mueller, N. D. et al. Closing yield gaps through nutrient and water management. *Nature* **490**, 254–257 (2012).

3. Hong, C. et al. Global and regional drivers of land-use emissions in 1961–2017. *Nature* **589**, 554–561 (2021).

4. Laurance, W. F., Sayer, J. & Cassman, K. G. Agricultural expansion and its impacts on tropical nature. *Trends Ecol. Evol.* **29**, 107–116 (2014).

5. Cassman, K. G. & Grassini, P. A global perspective on sustainable intensification research. *Nat. Sustain.* **3**, 262–268 (2020).

6. Hodge, I., Hauck, J. & Bonn, A. The alignment of agricultural and nature conservation policies in the European Union. *Conserv. Biol.* **29**, 996–1005 (2015).

7. Heilmayr, R., Rausch, L. L., Munger, J. & Gibbs, H. K. Brazil's Amazon Soy Moratorium reduced deforestation. *Nat. Food* **1**, 801–810 (2020).

8. Diffenbaugh, N. S. et al. Quantifying the influence of global warming on unprecedented extreme climate events. *Proc. Natl Acad. Sci. USA* **114**, 4881–4886 (2017).

9. Iizumi, T. & Ramankutty, N. How do weather and climate influence cropping area and intensity? *Glob. Food Security* **4**, 46–50 (2015).

10. Davis, K. F., Downs, S. & Gephart, J. A. Towards food supply chain resilience to environmental shocks. *Nat. Food* **2**, 54–65 (2020).

11. Wang, X. et al. Emergent constraint on crop yield response to warmer temperature from field experiments. *Nat. Sustain.* **3**, 908–916 (2020).

12. Lobell, D. B., Schlenker, W. & Costa-Roberts, J. Climate trends and global crop production since 1980. *Science* **333**, 616–620 (2011).

13. Sloat, L. L. et al. Climate adaptation by crop migration. *Nat. Commun.* **11**, 1243 (2020).

14. Afifi, T., Liwenga, E. & Kwezi, L. Rainfall-induced crop failure, food insecurity and out-migration in Same-Kilimanjaro, Tanzania. *Clim. Dev.* **6**, 53–60 (2014).

15. Stigter, K. in *Applied Agrometeorology* (ed. Stigter, K.) 531–534 (Springer, 2010).

16. Seifert, C. A. & Lobell, D. B. Response of double cropping suitability to climate change in the United States. *Environ. Res. Lett.* **10**, 024002 (2015).

17. Kawasaki, K. Two harvests are better than one: double cropping as a strategy for climate change adaptation. *Am. J. Agr. Econ.* **101**, 172–192 (2019).

18. Ceglar, A., Zampieri, M., Toreti, A. & Dentener, F. Observed northward migration of agro-climate zones in Europe will further accelerate under climate change. *Earths Future* **7**, 1088–1101 (2019).

19. Cohn, A. S., VanWey, L. K., Spera, S. A. & Mustard, J. F. Cropping frequency and area response to climate variability can exceed yield response. *Nat. Clim. Change* **6**, 601–604 (2016).

20. Challinor, A. J., Simelton, E. S., Fraser, E. D. G., Hemming, D. & Collins, M. Increased crop failure due to climate change: assessing adaptation options using models and socio-economic data for wheat in China. *Environ. Res. Lett.* **5**, 034012 (2010).

21. Ray, D. K. & Foley, J. A. Increasing global crop harvest frequency: recent trends and future directions. *Environ. Res. Lett.* **8**, 044041 (2013).

22. Wu, W. et al. Global cropping intensity gaps: increasing food production without cropland expansion. *Land Use Policy* **76**, 515–525 (2018).
23. Pugh, T. A. M. et al. Climate analogues suggest limited potential for intensification of production on current croplands under climate change. *Nat. Commun.* **7**, 12608 (2016).

24. Scherer, L. A., Verburg, P. H. & Schulp, C. J. E. Opportunities for sustainable intensification in European agriculture. *Glob. Environ. Change* **48**, 43–55 (2018).

25. Qin, Y. et al. Agricultural risks from changing snowmelt. *Nat. Clim. Change* **10**, 459–465 (2020).

26. Waha, K. et al. Multiple cropping systems of the world and the potential for increasing cropping intensity. *Glob. Environ. Change* **64**, 102131 (2020).

27. Raderschall, C. A., Vico, G., Lundin, O., Taylor, A. R. & Bommarco, R. Water stress and insect herbivory interactively reduce crop yield while the insect pollination benefit is conserved. *Glob. Chang. Biol.* **27**, 71–83 (2021).

28. Ding, M. et al. Variation in cropping intensity in Northern China from 1982 to 2012 based on GIMMS-NDVI data. *Sustainability* **8**, 1123 (2016).

29. Yu, Q., Xiang, M., Sun, Z. & Wu, W. The complexity of measuring cropland use intensity: an empirical study. *Agr. Syst.* **192**, 103180 (2021).

30. Moore, F. C. & Lobell, D. B. Adaptation potential of European agriculture in response to climate change. *Nat. Clim. Change* **4**, 610–614 (2014).

31. Agnolucci, P. et al. Impacts of rising temperatures and farm management practices on global yields of 18 crops. *Nat. Food* **1**, 562–571 (2020).

32. Zhu, P. & Burney, J. Temperature-driven harvest decisions amplify US winter wheat loss under climate warming. *Glob. Change Biol.* **27**, 550–562 (2021).

33. Ortiz-Bobea, A., Knippenberg, E. & Chambers, R. G. Growing climatic sensitivity of U.S. agriculture linked to technological change and regional specialization. *Sci. Adv.* **4**, 4343 (2018).

34. Duku, C., Zwart, S. J. & Hein, L. Impacts of climate change on cropping patterns in a tropical, sub-humid watershed. *PLoS ONE* **13**, e0192642 (2018).

35. Folberth, C. et al. The global cropland-sparing potential of high-yield farming. *Nat. Sustain.* **3**, 281–289 (2020).

36. Lobell, D. B. et al. The critical role of extreme heat for maize production in the United States. *Nat. Clim. Change* **3**, 497–501 (2013).

37. Yang, X. et al. Potential benefits of climate change for crop productivity in China. *Agric. For. Meteorol.* **208**, 76–84 (2015).

38. Burney, J., Woltering, L. & Burke, M. Solar-powered drip irrigation enhances food security in the Sudano-Sahel. *Proc. Natl Acad. Sci. USA* **107**, 1848–1853 (2010).

39. You, L. et al. What is the irrigation potential for Africa? A combined biophysical and socioeconomic approach. *Food Policy* **36**, 770–782 (2011).

40. Zheng, B., Chenu, K., Fernanda Dreccer, M. & Chapman, S. C. Breeding for the future: what are the potential impacts of future frost and heat events on sowing and flowering time requirements for Australian bread wheat (*Triticum aestivum*) varieties? *Glob. Change Biol.* **18**, 2899–2914 (2012).

41. Flach, R., Fader, M., Folberth, C., Skalský, R. & Jantke, K. The effects of cropping intensity and cropland expansion of Brazilian soybean production on green water flows. *Environ. Res. Commun.* **2**, 071001 (2020).

42. Wood, S. A., Jina, A. S., Jain, M., Kristjanson, P. & DeFries, R. S. Smallholder farmer cropping decisions related to climate variability across multiple regions. *Glob. Environ. Change* **25**, 163–172 (2014).

43. Paola, A. D. et al. The expansion of wheat thermal suitability of Russia in response to climate change. *Land Use Policy* **78**, 70–77 (2018).

44. Brunelle, T. & Makowski, D. Assessing whether the best land is cultivated first: a quantile analysis. *PLoS ONE* **15**, e0242222 (2020).

45. Lark, T. J., Spawn, S. A., Bougie, M. & Gibbs, H. K. Cropland expansion in the United States produces marginal yields at high costs to wildlife. *Nat. Commun.* **11**, 4295 (2020).

46. Zabel, F., Putzenlechner, B. & Mauser, W. Global agricultural land resources—a high resolution suitability evaluation and its perspectives until 2100 under climate change conditions. *PLoS ONE* **9**, e107522 (2014).

47. Petkeviciene, B. The effects of climate factors on sugar beet early sowing timing. *Agron. Res.* **7**, 436–443 (2009).

48. Ainsworth, E. A. & Long, S. P. 30 years of free-air carbon dioxide enrichment (FACE): what have we learned about future crop productivity and its potential for adaptation? *Glob. Change Biol.* **27**, 27–49 (2021).

49. Collins, M. et al. in *Climate Change 2013: The Physical Science Basis* (eds Stocker, T. F. et al.) Ch. 11 (Cambridge Univ. Press, 2013).

50. Pendergrass, A. G., Knutti, R., Lehner, F., Deser, C. & Sanderson, B. M. Precipitation variability increases in a warmer climate. *Sci. Rep.* **7**, 17966 (2017).

51. Asadieh, B. & Krakauer, N. Y. Global trends in extreme precipitation: climate models versus observations. *Hydrol. Earth Syst.* **19**, 877–891 (2015).

52. Zhang, Y., You, L., Lee, D. & Block, P. Integrating climate prediction and regionalization into an agro-economic model to guide agricultural planning. *Clim. Change* **158**, 435–451 (2020).

53. Turner, S. W. D., Hejazi, M., Yonkofski, C., Kim, S. H. & Kyle, P. Influence of groundwater extraction costs and resource depletion limits on simulated global nonrenewable water withdrawals over the twenty-first century. *Earths Future* **7**, 123–135 (2019).

54. Zhu, W., Jia, S., Devineni, N., Lv, A. & Lall, U. Evaluating China's water security for food production: the role of rainfall and irrigation. *Geophys. Res. Lett.* **46**, 11155–11166 (2019).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2022
# Methods

## FAO agricultural statistics

Crop production, harvested area, crop yield, CA, irrigated CA and agricultural total nitrogen use at country scale spanning 1979–2018 were obtained from the FAOSTAT database$^{53}$. We restricted our analysis to 1979–2018 as FAO statistics during the earlier time period include more missing or estimated values and also because we used an hourly climate dataset starting from 1979. The 150 most populous countries are selected as they normally have more complete data records and cover 98% of global crop production. In the FAOSTAT database, harvested area and crop yield are recorded for 161 crop species, while CA, irrigated CA and agricultural total nitrogen use are the sum of all crop species and do not distinguish among crops. According to FAO’s definition for harvested area, if the crop is harvested more than once during the year as a consequence of successive cropping (for example, rice is sown and harvested more than once in the same field during the year), the area is counted as many times as harvested, which ensures multiple-cropping is accounted for in harvested area data. In terms of CA, it is the sum of areas under ‘arable land’ and ‘permanent crops’, according to FAO definitions. Arable land is the land area under temporary crops, temporary meadows and pastures and land with temporary fallow. Land under permanent crops means land cultivated with long-term crops which do not have to be replanted for several years (such as cocoa and coffee), land under trees and shrubs producing flowers and nurseries (except those for forest trees, which should be classified under ‘Forestry’). We note that FAO data were a mixture of ‘FAO estimate’ and ‘official data’. Data entries marked as ‘FAO estimate’ were often less reliable than the ‘official data’. To minimize the influence of those countries with high rates of ‘FAO estimate’ on our regression model, countries with >20% of the data on harvested area or production marked as ‘estimated data’ were excluded, as previously done$^{56}$. The excluded nations were Guinea, Kenya, Mozambique, North Korea and Zambia.

## Modelling climate effects on CalY, CF, CA and CalP

For a country producing $N$ different crops ($N$ denotes the number of crop species in a country), the total yearly crop (CalP) can be further decomposed into CalY, CF and CA:

$$ \text{CalP} = \sum_{k=1}^{N} \gamma_k \times \text{Yield}_k \times \text{HA}_k = \frac{\sum_{k=1}^{N} \gamma_k \times \text{Yield}_k \times \text{HA}_k}{\sum_{k=1}^{N} \text{HA}_k} \times \frac{\sum_{k=1}^{N} \text{HA}_k}{\text{CA}} \times \text{CA} $$

where $\gamma_k$ is the caloric conversion factor of crop $k$, and we obtained the caloric conversion factor on the basis of published dataset (http://www.fao.org/docrep/003/X9892E/X9892e05.htm#P8217_125315); HA$_{k}$ indicates harvested area of crop $k$; CA is the total CA of a country. Thus, $\frac{\sum_{k=1}^{N} \gamma_k \times \text{Yield}_k \times \text{HA}_k}{\sum_{k=1}^{N} \text{HA}_k}$ represents the country-average CalY and $\frac{\sum_{k=1}^{N} \text{HA}_k}{\text{CA}}$ represents the country-scale average CF (Supplementary Fig. 1). This country-scale average CF is the aggregate of fallow croplands (frequency = 0) due to conservation purposes or abandonment, single-cropping per year (frequency = 1) or multiple-cropping per year (CF > 1). On the basis of these definitions, we have $\text{CalP} = \text{CalY} \times \text{CF} \times \text{CA}$.

We build panel models separately for each dependent variable (CF, CalY, CA or CalP) to estimate the effects of year-to-year climate variation and farmer management practices on CF, CalY, CA or CalP (sample size $n = 5,184$), which is specified as follows:

$$ \log(Y_{i,t}) = \alpha_{1,t}t + \alpha_{2,t}t^2 + \text{country}_i + \beta W_{i,t} + \varepsilon_{i,t} $$

$Y_{i,t}$ represents CF, CalY, CA or CalP for country $i$ and year $t$. Term $\alpha_{1,t} + \alpha_{2,t}t^2$ characterizes the country-specific quadratic time trends, which capture unobserved, country-specific factors affecting CF, CalY, CA or CalP, such as technological progress. Parameter country$_i$ is the country-specific fixed-effect capturing all time-invariant,

country-specific factors that might explain variations in $Y_{i,t}$. The model component $\beta W_{i,t}$ takes into account the effects of the inputs in $W_{i,t}$ potentially affecting $Y_{i,t}$ using a set of parameters $\beta$ (common to all countries$^{12,31}$). Our baseline model for $\beta W_{i,t}$ used a quadratic function of annual mean temperature (Tmean) and annual precipitation (Prcp) to characterize the potential nonlinear effect of Tmean and Prcp on $Y_{i,t}$:

$$ \beta W_{i,t} = \beta_1 \text{Tmean}_{i,t}^2 + \beta_2 \text{Tmean}_{i,t} + \beta_3 \text{Tmean}_{i,t} \cdot \text{Irri}_{i,t} + \beta_4 \text{Prcp}_{i,t}^2 + \beta_5 \text{Prcp}_{i,t} + \beta_6 \text{Prcp}_{i,t} \cdot \text{Irri}_{i,t} + \beta_7 \text{Fert}_{i,t} $$

To account for the potential compensation effect of irrigation on temperature stress and water stress, irrigation fraction (Irri$_{i,t}$) (Supplementary Fig. 2a), which was estimated as the ratio of irrigated CA to total CA at country $i$ and year $t$, was interacted with Tmean and Prcp. The interaction terms characterize the marginal relationship between CF, CalY or CalP and irrigation, conditional on climate variables. This allows us to quantify the benefit of irrigation through evaluating the sensitivity of $Y_{i,t}$ to irrigation fraction ($\frac{\partial Y_{i,t}}{\partial \text{Irri}_{i,t}} = \beta_3 \text{Tmean}_{i,t} + \beta_6 \text{Prcp}_{i,t}$). To control the influence of fertilizer application on increasing CF or crop yield, nitrogen application rate (Fert$_{i,t}$) was included in the model as well. Parameter Fert$_{i,t}$ was estimated as the ratio of agricultural total nitrogen use to the total harvested area of all crops in a country. Our panel model was weighted by the country-level crop area from FAO to define the model output as an average over all crop areas. The weighting method was also useful to reduce heteroskedasticity and correct the influences of countries with very small CAs$^{30}$. Since presence of spatial autocorrelation could violate the assumption of normally distributed residuals, we also checked whether spatial autocorrelation occurred in each panel model residual using Moran’s $I$. The results suggest that there is no statistically significant ($P > 0.05$) spatial autocorrelation in each panel model. Further, we note that although CF might be also driven by other unobserved variables (crop rotation, cropland expansion and agricultural policies) apart from climatic factors, as long as these unobserved variables are not marginally influenced by temperature, which seems plausible, the first-order approximation of temperature effect on CF can be estimated with our linear regression model. With the established model and the historical irrigation fraction, the marginal effect of temperature on CF, CalY, CA or CalP was estimated as the difference between predicted $Y_{i,t}$ with 1 °C uniform increase in Tmean and the original $Y_{i,t}$. The global average change is then a weighted average (weighted by country crop area) of the 1 °C warming effects in each country.

Hourly reanalysis climate data ERA5 at 0.25° × 0.25° resolution spanning from 1979 to 2018 were used to characterize climate conditions during crop growth$^{57}$. Considering some crops were sown in the previous year of harvest, especially for winter crops in the Northern Hemisphere, we used both the current year and the second half of the previous year to calculate Tmean and Prcp. In addition, we also fitted a panel model with climate variables exclusively averaged over the current year (see model summary in Supplementary Table 5). When averaging climate variables to get their annual mean, we did not distinguish crop growing season and non-growing season because the 161 crop species considered here have diverse growing seasons and also crop yield or farmer planting decision might be influenced by both crop growing season and non-growing season weather$^{58}$. The annual mean temperature and precipitation at 0.25° × 0.25° resolution were then aggregated to country-level spatially weighted by CA fraction in each grid cell. We get CA fraction in each grid cell through aggregating the 1 km Global Food Security Support Analysis Data (GFSAD) crop mask into 0.25° × 0.25° resolution$^{59}$ (Supplementary Fig. 3). GFSAD dataset was created using multiple input data including remote sensing such as Landsat, Satellite Probatoire d'Observation de la Terre (SPOT) vegetation and MODIS; multiyear precipitation and temperature data; ground reference data; and country statistics data$^{59}$.
## Robustness checks

We tested other specifications of $\beta W_{i,t}$ to assess the robustness of our analysis. Similarly, nitrogen application rate and interaction between $Irri_{i,t}$ and climate variables were also considered in the following alternative specifications of $\beta W_{i,t}$ to account for fertilizer application and irrigation compensation effects, respectively.

M2: Model using natural cubic splines (NCS) function of annual mean temperature and precipitation. Studies suggest NCS function has a good capacity to capture the nonlinear response of crop productivity to climate variables$^{60}$. Here, we used five knots set at the 5th, 25th, 50th, 75th and 95th percentiles of $Tmean$ or $Prcp$:

$$ \beta W_{i,t} = \beta_1 \text{NCS}(Tmean_{i,t}) + \beta_2 Tmean_{i,t} \times Irri_{i,t} + \beta_3 \text{NCS}(Prcp_{i,t}) + \beta_4 Prcp_{i,t} \times Irri_{i,t} + \beta_5 \text{Fert}_{i,t} $$

$Tmean_{i,t,m}$ and $Prcp_{i,t,m}$ are the annual mean temperature and precipitation for country $i$ and year $t$.

M3: Model with quadratic function of 2-month mean $Tmean$ and $Prcp$. Studies have suggested that crop yield is sensitive to climate stresses in specific crop growth stages. To better characterize the influence of seasonal climate variation on CF or CalY, we used the quadratic function of 2-month average climate variables:

$$ \beta W_{i,t} = \sum_{m=1}^{9} (\beta_{1,m} Tmean_{i,t,m}^2 + \beta_{2,m} Tmean_{i,t,m} + \beta_{3,m} Tmean_{i,t,m} \times Irri_{i,t} + \beta_{4,m} Prcp_{i,t,m}^2 + \beta_{5,m} Prcp_{i,t,m} + \beta_{6,m} Prcp_{i,t,m} \times Irri_{i,t}) + \beta_{7} \text{Fert}_{i,t} $$

$Tmean_{i,t,m}$ and $Prcp_{i,t,m}$ are the $m$th 2-month mean temperature and precipitation for country $i$ and year $t$. As we used climate variables in both current year and the second half of the previous year, there are nine 2-month mean $Tmean_{i,t,m}$ and $Prcp_{i,t,m}$. Here, we did not use monthly climate variables in the model to avoid too many predictors and thus overfitting. Penalized regression (LASSO) was implemented with the R package glmnet$^{61}$ to select the most influential climate predictors. A tenfold cross-validation was performed to maximize the predictive accuracy of the model. LASSO regression allowed us to identify a subset of predictors that explain most of the variation in outcomes by shrinking the regression coefficient towards zero and discarding irrelevant predictors. This procedure thus automatically determined which time periods of climate variables were most relevant for explaining CF, CalY and CalP variation.

M4: GDD model.

$$ \beta W_{i,t} = \sum_{m=1}^{5} (\beta_{1,m} \text{GDD}_{i,t,m} + \beta_{2,m} \text{GDD}_{i,t,m} \times Irri_{i,t} + \beta_{3} \text{Prcp}_{i,t}^2 + \beta_{4} \text{Prcp}_{i,t} + \beta_{5} \text{Prcp}_{i,t} \times Irri_{i,t} + \beta_{6} \text{Fert}_{i,t} $$

In this GDD model, five levels of GDD for country $i$ and year $t$ ($GDD_{i,t,i}$, $GDD_{i,t,2}$, $GDD_{i,t,3}$, $GDD_{i,t,4}$ and $GDD_{i,t,5}$) were considered to characterize the differential response of $Y_{i,t}$ to different levels of temperature exposures with 10 °C increment: $GDD_1$ as $GDD_{0-\infty}$, $GDD_2$ as $GDD_{0^0}$, $GDD_3$ as $GDD_{20^0}$, $GDD_4$ as $GDD_{30^0}$ and $GDD_5$ as $GDD_{30^+}$. ERA5 hourly temperature at the height of 2 m was used to estimate GDD with the following equation:

$$ \text{GDD}_{T_1}^{T_2} = \sum_{h=1}^{H} \text{DD}_h / 24, \text{DD}_h = \begin{cases} 0, & T_h < T_1 \\ T_h - T_1, & T_1 \le T_h < T_2 \\ T_2 - T_1, & T_h \ge T_2 \end{cases} $$

where $T_h$ represents the temperature at hour $h$, and $H$ is the total number of hours during the growing season. $T_1$ and $T_2$ indicate the lower and upper temperature thresholds of GDD, respectively. GDD$_1$ takes the

freezing stress into account, while GDD$_5$ uses a high temperature threshold of 30 °C to characterize high-temperature stress effects on crop growth as suggested by previous studies$^{32,62,63}$. The other GDDs generally represent the mild temperature exposures. Similarly, we used LASSO regression to estimate model coefficients, as LASSO regression can automatically select GDDs most relevant for explaining CF, CalY, CA and CalP variations, thus minimizing the potential multicollinearity effect among GDDs.

M5 and M6: Model with alternative climate dataset. To test the robustness of our model to other climate datasets, we also run the baseline model with two other climate datasets: Climate Research Unit time-series datasets (CRU TS 4.0.4)$^{64}$ (M5 in Supplementary Table 5) and the University of Delaware temperature and rainfall datasets$^{65}$ (M6 in Supplementary Table 5).

M7: Model with annual climate variables aggregated over 1 year. This model uses the same specification as the baseline model M1, but the annual climate variables are aggregated over 1 year (12-month) time period.

M8: Total factor productivity input (TFPI) model (M8 in Supplementary Table 5).

$$ \beta W_{i,t} = \beta_1 Tmean_{i,t}^2 + \beta_2 Tmean_{i,t} + \beta_3 Tmean_{i,t} \times \text{TFPI}_{i,t} + \beta_4 \text{Prcp}_{i,t}^2 + \beta_5 \text{Prcp}_{i,t} + \beta_6 \text{Prcp}_{i,t} \times \text{TFPI}_{i,t} $$

In this alternative model, agricultural total factor productivity input (TFPI) replaces irrigation and fertilizer application to represent the comprehensive management practices (cultivation technology, management ability and agricultural infrastructure). Total factor productivity (TFP) measures productive efficiency, that is, the amount of agricultural output produced from the combined set of input$^{66}$. The output includes crop and livestock commodities aggregated on the basis of a common set of international prices derived by the FAO. Input (TFPI) includes agricultural land, farm labour, irrigation, capital inputs (including farm machinery) and intermediate inputs (fertilizer). TFPI thus can be used to represent the comprehensive management practices used in farm production. This dataset is obtained from the United States Department of Agriculture (USDA) Economic Research Service (ERS) International Agricultural TFP dataset (https://www.ers.usda.gov/data-products/international-agricultural-productivity/), which provides country-level TFP and TFPI index for 172 countries over the 1961–2019 period. In model M8, TFPI replaces irrigation and fertilizer application to more comprehensively account for various management practices. The interaction term between TFPI and climate variables is used to characterize the effect of agricultural management practices in offsetting the negative impact of climate change on yield and CF. To account for the statistical uncertainty of these regression models, we run each model with 1,000 bootstrap, where we sample from all 5,184 country-year observations with replacement. The 1,000 sets of estimated regression coefficients are then used to determine the CI of the model-estimated coefficients. The medians of the 1,000 sets of coefficients are used to determine the response curves in Fig. 2a–c. With the median of coefficients $\beta$, two response curves corresponding to without and with irrigation can be determined through setting irrigation fraction as zero and historical mean (3.1%), respectively (Fig. 2a–c). When LASSO regression is run using each bootstrap sample, the penalized algorithm selects new combinations of predictors for each sample. Therefore, applying bootstrap to LASSO regression model (M3 and M4) allows us to account for the uncertainties in model formulation as well.

We also checked whether the estimated coefficients will be biased by countries with extreme high or low temperatures. With the baseline model, model coefficients and marginal effects of 1 °C warming were re-estimated on the basis of samples excluding the 10% coldest or 10% warmest countries. Another potential bias is the omission of certain
crops in FAO statistics, considering there might be more than 161 crops planted in a country²¹. To simulate the influence of crop species omission on our estimations, five levels (10%, 20%, 30%, 40%, 50%) of crop species omission were considered. For example, to simulate the effects of 10% crop species omission, we re-estimated CF, CalY and CalP 1,000 times, and for each time, we randomly discarded 10% crop species. We then applied the baseline model M1 to each set of CF, CalY and CalP to estimate climate effects (1,000 estimates of climate effects) (Supplementary Fig. 18).

## Model evaluation

First, we conduct an out-of-sample model validation with 1,000 bootstraps. Each time, 80% of country-year observations are randomly sampled as the training data to build the panel model with model M1; the remaining 20% country-year observations are used as test data. The relative difference between the model predicted CF, CalY and CalP and those in test data were used to evaluate the performance of our panel model. Relative differences in each country are presented in Supplementary Fig. 19.

Second, we fit the panel model M1 using a subset of dataset flagged by the FAO as 'official data'. We consider this to be as high-quality data as possible (compared to, for example, data flagged as 'estimated data'). Similar to Fig. 2 in the main text, we estimate the temperature sensitivity of CF, CalY and CalP with model M1 and create the response curves of CF, CalY and CalP to temperature. This alternative model test has a smaller sample size but does not substantially change the response of CF, CalY and CalP to temperature (Supplementary Fig 13) and consistently suggests an overall negative association between climate warming and CF, CalY and CalP at global scale.

Third, we validate the country-level aggregated CF with subnational statistics in Brazil and satellite data-derived CF. The subnational statistics in Brazil are compiled in a recent study⁶⁷ and characterize the corn-soybean double-cropping system. The statistics contain soybean harvest area and also separate records for first-season corn and second-season corn harvest area. Since the main-season crops are soybean and first-season corn, whereas the only second crop is second-season corn, the CF can be thus calculated as:

$$CF = \frac{\text{Area}_{c1} + \text{Area}_{c2} + \text{Area}_{\text{soy}}}{\text{Area}_{c1} + \text{Area}_{\text{soy}}}$$

where Area$_{c1}$ and Area$_{c2}$ are the harvest area of first-season and second-season corn, respectively, and Area$_{soy}$ is the harvested area of soybean. The validation result can be found in Supplementary Fig. 14.

A spatially explicit estimation of CF derived from satellite data is also used. Relative to the CF derived from national statistics, it can better characterize the spatial heterogeneity of multiple-cropping practices. Similar to the CF used in the previous study¹⁹, our satellite-derived CF is also developed by detecting the number of sharp peaks followed by troughs using time series of MODIS enhanced vegetation index (EVI). This information is obtained from the established MODIS land cover dynamics product (MCD12Q)⁶⁸ at 500 m spatial resolution during 2001–2018, wherein it contains a data layer named ‘NumCycles’. Combining MODIS-derived ‘NumCycles’ which characterizes the total number of vegetation cycles per year with MODIS land cover map⁶⁹ which differentiates cropland and non-cropland, we can get the CF information for each 500 m grid cell across the global cropland. Then this satellite-derived CF information is averaged to country level and used as the dependent variable in M1, by which we can check whether our findings will hold with this satellite-based CF. These data evaluation and model test results can be found in Supplementary Figs. 14 and 15.

## Future projections

We projected future CF, CalY and CalP during 2031–2070 using the baseline model M1 with climate outputs from Coupled Model

Intercomparison Project Phase 6 (CMIP6), which is provided by Inter-Sectorial Impact Model Intercomparison Project (ISIMIP). Climate change effects are estimated as the difference between predicted CF, CalY or CalP during 2031–2070 and the reference period (1979–2018) CF, CalY or CalP. Two climate scenarios were considered to represent lower (SSP 126) and higher (SSP 585) emission scenarios. Here, SSP 585 is used to shed light on future crop production changes under the worst-case climate change scenarios, since SSP 585 is commonly used to simulate upper limits of temperature change over the next century⁷⁰. ISIMIP provides five climate model outputs (GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, MRI-ESM2-0 and UKESM1-0-LL) to account for the potential bias in future climate projections in each specific model. The five climate models were selected because they are structurally independent in terms of their ocean and atmosphere model components. Meanwhile, the five primary models are good representatives of the whole CMIP6 ensemble as they include three models with low climate sensitivity (GFDL-ESM4, MPI-ESM1-2-HR and MRI-ESM2-0) and two models with high climate sensitivity (IPSL-CM6A-LR and UKESM1-0-LL)⁷¹. Climate model outputs (daily mean temperature and precipitation) were downscaled to 0.5° resolution and bias-corrected through comparing the climate model outputs with corresponding climate observations during the training period⁷¹. The bias-corrected climate model outputs were available at https://esg.pik-potsdam.de/search/isimip/. Then a similar method as we processed the historical climate dataset was applied to the five climate model outputs to obtain country-level climate variables for 2031–2070⁷². Our future projections assume no new adaptation between now and then or holding technology equivalent to current levels, so that the estimated changes in CF, CalY and CalP are exclusively caused by changes in temperature and precipitation.

## Data availability

FAO national statistical data was obtained from http://www.fao.org/faostat/en/#data/QC. Agricultural TFP was obtained from the USDA ERS International Agricultural TFP dataset https://www.ers.usda.gov/data-products/international-agricultural-productivity/. The caloric conversion factor is based on the published dataset http://www.fao.org/docrep/003/X9892E/X9892e05.htm#P8217_125315. The bias-corrected climate model outputs are available at https://esg.pik-potsdam.de/search/isimip/.

## Code availability

The scripts used to run the regression model are available through zenodo at: https://zenodo.org/record/7038556

## References

55. FAOSTAT (Food and Agriculture Organization of the United Nations, 1997).

56. Egli, L., Schröter, M., Scherber, C., Tscharntke, T. & Seppelt, R. Crop asynchrony stabilizes food production. Nature 588, E7–E12 (2020).

57. Hersbach, H. et al. ERA5 Hourly Data on Single Levels from 1979 to Present (Copernicus Climate Change Service (C3S) Climate Data Store (CDS), accessed 1 August 2020); https://doi.org/10.24381/cds.adbb2d47 (2018).

58. Feng, P. et al. Impacts of rainfall extremes on wheat yield in semi-arid cropping systems in eastern Australia. Clim. Change 147, 555–569 (2018).

59. Teluguntla, P. et al. in Land Resources Monitoring, Modeling, and Mapping with Remote Sensing (ed. Thenkabail, P. S.) 849 (CRC Press, 2015).

60. Hawkins, E. et al. Increasing influence of heat stress on French maize yields from the 1960s to the 2030s. Glob. Change Biol. 19, 937–947 (2013).
61. Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for generalized linear models via coordinate descent. *J. Stat. Softw.* **33**, 1–22 (2010).

62. Lobell, D. B., Bänziger, M., Magorokosho, C. & Vivek, B. Nonlinear heat effects on African maize as evidenced by historical yield trials. *Nat. Clim. Change* **1**, 42–45 (2011).

63. Deryng, D., Sacks, W. J., Barford, C. C. & Ramankutty, N. Simulating the effects of climate and agricultural management practices on global crop yield. *Glob. Biogeochem. Cycles* **25**, GB2006 (2011).

64. New, M., New, M., Lister, D., Hulme, M. & Makin, I. A high-resolution data set of surface climate over global land areas. *Clim. Res.* **21**, 1–25 (2002).

65. Willmott, C. J. *Terrestrial Air Temperature and Precipitation: Monthly and Annual Time Series (1950–1996)* (Center for Climatic Research, 2000); http://climate.geog.udel.edu/~climate/html_pages/README.ghcn_ts2.html

66. Van Beveren, I. Total factor productivity estimation: a practical review. *J. Econ. Surv.* **26**, 98–128 (2012).

67. Xu, J. et al. Double cropping and cropland expansion boost grain production in Brazil. *Nat. Food* **2**, 264–273 (2021).

68. Friedl, M. & Gray, J. *MCD12Q2 MODIS/Terra+ Aqua Land Cover Dynamics Yearly L3 Global 500 m SIN Grid V006* (NASA EOSDIS, 2019).

69. Sulla-Menashe, D. & Friedl, M. A. *User Guide to Collection 6 MODIS Land Cover (MCD12Q1 and MCD12C1) Product* (USGS, 2018).

70. Schwalm, C. R., Glendon, S. & Duffy, P. B. RCP8.5 tracks cumulative CO2 emissions. *Proc. Natl Acad. Sci. USA* **117**, 19656–19657 (2020).

71. Lange, S. Trend-preserving bias adjustment and statistical downscaling with ISIMIP3BASD (v1.0). *Geosci. Model Dev.* **12**, 3055–3070 (2019).

72. Peng Zhu. Climate effects on caloric yield and cropping frequency. Zenodo https://doi.org/10.5281/zenodo.7038556 (2022).

## Acknowledgements

P.Z. and P.C. are supported by the CLAND project (grant no. 16-CONV-0003) and ISIPEDIA: The Open Inter-Sectoral Impacts Encyclopedia (grant no. ANR-17-ERA4-0006 - ISIPEDIA). D.M. is supported by the CLAND project (grant no. 16-CONV-0003) and meta-programme CLIMAE-INRAE. L.Y. is supported by Tsinghua University Initiative Scientific Research Programme (2021Z11GHX002, 20223080017). J.B. is supported by NSF/NIFA no. 1639318 INFEWS/T1. J.C. is supported by the National Key Research and Development Programme of China (2021YFE0114500). Q.X. is supported by National Key Research and Development Programme of China (grant no. 2017YFA0604300).

## Author contributions

P.Z. designed the study, performed the analysis and led the writing. J.B., J.C., Z.J., N. M., D.M. and P.C. helped the results interpretation. Q.X., J.X. and L.Y. provided additional data for comparison. All authors reviewed the manuscript and contributed to the manuscript writing.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41558-022-01492-5.

**Correspondence and requests for materials** should be addressed to Peng Zhu.

**Peer review information** *Nature Climate Change* thanks Avery Cohn, Katharina Waha and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.