Meta-analysis of yield response to climate change

The yield–temperature response functions used in this paper are derived from a database of studies estimating the climate change impact on yield compiled for the IPCC 5th Assessment Report17, also described in a meta-analysis by Challinor et al.16. Methods for sampling the literature and criteria for inclusion are described in Challinor et al.16. as a “broad and inclusive literature search” combined with quality-control procedures documented in the Supplemental Information of that paper. In this study, we focus on four major crops—maize, rice, wheat, and soybeans. The bulk of the scientific literature on yield response to temperature relates to these crops, which collectively account for about 20% of the value of global agricultural production, 65% of harvested crop area, and nearly 50% of calories directly consumed33. For the four crops, the database contains 1010 observations (344, 238, 336, and 92 for maize, rice, wheat, and soybeans, respectively) from 56 different studies (many studies report multiple yield changes for different crops, different locations, different levels of temperature change, or different assumptions about adaptation). The studies include 8 empirical studies and 48 process-based studies, published between 1997 and 2012. Supplementary Figs 14 and 15 show the geographic coverage of production areas within the database and the distribution of publication dates. Of the 1010 data points, 451 are reported as including some form of on-farm, within-crop, agronomic adaptation. The vast majority of these adaptations involve adjusting either planting date (10%) or cultivar (12%) or both (44%). Recognizing the existence of a more recent and possibly more systematic literature review for wheat yields, we perform a robustness check where we incorporate additional results identified in Wilcox and Makowski34. This substantially increases the number of observations for wheat but does not affect our estimated response curve (Supplementary Fig. 16).

We merge this database with information on baseline growing-season temperature for each data point. To do this, each data point was assigned to a country. For the 14% of studies looking at more than one country, the country assigned was the one with the highest production of the relevant crop. Average baseline growing-season temperatures were calculated using planting and harvest dates from Sacks et al35. and gridded monthly temperatures for 1979–2013 from the Climate Research Unit36. These were averaged to the country level using year 2000 crop production weights from Monfreda et al37.

The response functions are jointly estimated from the point-estimates in the database using a multi-variate:

$$\begin{array}{ccccc}\\ \Delta Y_{ijk} = {\mathrm{\beta }}_{1{{j}}}\Delta T_{ijk}{\mathrm{*}}Crop_j + {\mathrm{\beta }}_{2{{j}}}\Delta T_{ijk}^2{\mathrm{*}}Crop_j + {\mathrm{\beta }}_{3{\mathrm{j}}}\Delta T_{ijk}{\mathrm{*}}Crop_j{\mathrm{*}}\bar T_{jk} \\ + {\mathrm{\beta }}_{4{{j}}}\Delta T_{ijk}^2{\mathrm{*}}Crop_j{\mathrm{*}}\bar T_{jk} + \beta _5f_1\left( {\Delta {\rm CO}_{2ijk}} \right){\mathrm{*}}C_{3j} \\ + \beta _6f_2\left( {\Delta {\rm CO}_{2ijk}} \right){\mathrm{*}}C_{4j} + \beta _7\Delta P_{ijk} + \beta _8\Delta T_{ijk}{\mathrm{*}}Adapt_{ijk} + \beta _9Adapt_{ijk} + \varepsilon _{ijk}\end{array}$$ (1)

where \(\Delta Y_{ijk}\) is the change in yield from point-estimate i for crop j in country k (in %). \(\Delta T_{ijk},\Delta CO_{2ijk}\) and \(\Delta P_{ijk}\) are the changes in temperature (in degree C), CO 2 concentration (in parts per million (ppm)), and rainfall (in percent) for point-estimate ijk, \(\bar T_{jk}\) is the baseline growing-season temperature for crop j in country k, \(C_{3j}\) and \(C_{4j}\) are dummy variables indicating whether crop j is C 3 or C 4 , and \(Adapt_{ijk}\) is a dummy variable indicating whether the point-estimate includes any on-farm adaptation. Eq. 1 is estimated using an ordinary least squares regression.

Uncertainty in the parameters is estimated through 1500 block bootstraps, with blocks defined at the study level, allowing for possible correlation between point-estimates from the same study. Error bars reported throughout the paper are based on the 2.5th and 97.5th quantiles of the bootstrapped distribution. This treatment of the errors does assume independence between studies, which may be questionable if the same model is used in multiple studies. In total, 28 models, made up of 17 process-based model families (i.e., treating CERES-maize, CERES-rice, and CERES-wheat as a single model) and 11 statistical models, are used in the 56 studies. Supplementary Fig. 17 shows response curves with standard errors based on a model block bootstrap as a robustness check. These are qualitatively similar to the error bars shows in Fig. 1, particularly for warming <3 °C that is the focus of the economic analysis, suggesting the study block bootstrap is capturing the bulk of residual covariance. All error bars reported in the paper show confidence intervals rather than prediction intervals. This is appropriate since the relevant uncertainty is in the expected response of yield to temperature change, which is given by confidence intervals.

There are a number of important things to note about this specification shown in Eq. 1. First, the impacts of temperature are modeled as crop-specific quadratics (\({\mathrm{\beta }}_{1{{j}}}\) and \({\mathrm{\beta }}_{2{{j}}}\) terms), allowing the effects of warming to vary by crop. In addition, the effects of warming are allowed to vary with baseline growing-season temperature (\({\mathrm{\beta }}_{3{{j}}}\) and \({\mathrm{\beta }}_{4{{j}}}\) terms), capturing the intuition that the impacts of a 1 °C warming should be different in a cold location than in a hot location.

Second, there is no intercept term, thereby forcing response functions without adaptation through the origin. This is consistent with the expected functional form of a climate damage function, which should have no impacts if there are no changes in climate variables. However, we include an intercept for studies that do include adaptation (\(\beta _9\)). This is prompted by the observation that, in many studies, ‘adaptation’ is represented by changing management practices that would improve yields even in the current climate, such as adoption of improved varieties or increasing fertilizer or irrigation inputs22. Failing to include an adaptation intercept in this context will lead to an overestimation of the potential of these kinds of changes to reduce the negative impacts of a warming climate. This adaptation intercept is subtracted in our estimates of the effect of climate on yield to produce an adjusted damage function that goes through the origin. (In other words, we calculate the effect of a change in temperature of X on yields to be the yield change predicted from Eq. 1 for a temperature change of X minus the yield change predicted for a temperature change of zero (i.e., \(\beta _6\,\,Adapt_{ij}\))). The true effect of adaptation is the interaction with temperature change, given by the \(\beta _8\) term in Eq. 1, which is included in all subsequent analyses. This term captures the effect of management changes that are not beneficial today but that will become beneficial under a changed climate, the standard definition of adaptation.

Finally, the functions \(f_1()\) and \(f_2()\) are concave, allowing for a declining marginal effect of CO 2 , consistent with a number of field studies20, 38. Specifically, the function takes the form \(f\left( {\Delta {\rm CO}_{2ij}} \right) = \frac{{\Delta {\rm CO}_{2ij}}}{{\Delta {\rm CO}_{2ij} + A}}\) where A is a free parameter set at 100 ppm for \(C_3\) crops and at 50 ppm for \(C_4\) crops based on a comparison of the R 2 across models using multiple possible values. The changes in CO 2 are adjusted so that all are relative to a modern baseline of 360 ppm (the most common baseline value for studies included in the analysis).

In addition to Eq. 1, our preferred specification, we investigate the effects of several alternate specifications. Specifically we first investigate whether newer studies (publication date of 2005 or later) give a different temperature response compared to the full sample; second investigate the effect of individual agronomic adaptations, specifically changing cultivar and planting date; third allow the effect of temperature to differ depending on whether the study was a process-based or empirical study; fourth add a \(\Delta T_{ijk}^3\) term in the specification; and finally perform F-tests on individual terms within the model. These findings are documented in the SI (Supplementary Figs 18–20, Supplementary Tables 6 and 7). They do not substantially alter our estimates of the yield response to climate change.

Gridded yield changes

After estimating Eq. 1, we developed global gridded yield change scenarios for the four major crops (maize, wheat, rice, and soybeans). Although IAM damage functions are typically based on global temperature changes, it is important to account for the fact that local warming may differ significantly from global warming in estimating impacts. Local yield impacts will depend on local temperature changes, which scale in a predictable way with global temperature change. We estimate this scaling using the CMIP5 multi-model ensemble mean for the high emissions scenario RCP 8.539. For each grid cell, we take the change in temperature between a future (2035–2065) and baseline (1861–1900) period and divide by the mean global warming over this time period, giving the pattern scaling relationship between global and local temperature change for each grid cell (Supplementary Fig. 21). For a given increase in global mean temperature, warming is larger over land than over the ocean and at high latitudes compared to the tropics.

These gridded temperature changes are combined with the yield–temperature response function estimated using Eq. 1 and baseline growing-season temperature to give yield changes at different levels of global warming. We calculate yield changes for warming of 1, 2, and 3 degree Celsius including the estimated effect of on-farm adaptation. Any predicted yield losses >100% are set to losses of 99%. The CO 2 fertilization effect is included for all crops. CO 2 concentrations for a given level of global temperature change are determined based on a fitted quadratic relationship between global temperature change and CO 2 concentrations from the RCP 8.5 CMIP5 multi-model ensemble mean (adjusted R 2 > 0.999, 98 degrees of freedom).

AgMIP GGCMI ensemble

The AgMIP GGCMI is the one other source of global, multi-crop, multi-model yield responses and so we compare the results of our meta-analysis against these results. This ensemble of gridded crop model outputs includes up to seven process-based crop models, run using five General Circulation Models (GCMs)18. Yield changes are calculated relative to the 1981–2000 average. In order to determine yield changes for specific levels of temperature change (1–3 °C), we find the year in which warming passes each specific level for each GCM for the RCP 8.5 emissions scenario (taking the average of multiple ensemble members, if available) and take the 11-year yield average around that year40. We determine irrigated areas using crop-specific irrigation areas from Monfreda et al37. and use irrigated results for cells where irrigated crop area exceeds non-irrigated crop area. We use runs including CO 2 fertilization for all analyses.

The results reported in the main text use a preferred AgMIP ensemble that only includes models that explicitly represent nitrogen stress (EPIC, GEPIC, PEGASUS, and pDSSAT). We believe these results are preferred given crop response to changing temperature and, in particular, CO 2 conditions is known to depend on nutrient availability41, 42 and crops in many areas of the world are currently under-fertilized43. Moreover, the distinction between models based on representation of nitrogen stress has been identified as significant in understanding ensemble results by the AgMIP team18. Results in the main text should therefore be interpreted as impacts assuming continuation of current nutrient management practices. In the SI, we report results using the full AgMIP ensemble, which differ substantially from those of the restricted ensemble (Supplementary Fig. 11 and Supplementary Table 5). For both ensembles, the mean is calculated as a simple mean of yield change for each level of warming using all crop model×GCM combinations.

Welfare consequences of yield changes

To estimate the economic implications of warming-induced yield shocks, we use the Global Trade Analysis Project (GTAP) general equilibrium model and its accompanying database23, 24. GTAP is a widely used, comparative static general equilibrium model that exhaustively tracks bilateral trade flows between all countries in the world and explicitly models the consumption and production for all commodities of each national economy. Producers are assumed to maximize profits, while consumers maximize utility. Factor market clearing requires that supply equal demand for agricultural and non-agricultural skilled and unskilled labor and capital, natural resources, and agricultural land, and adjustments in each of these markets in response to the climate change shocks determines the resulting wage and rental rate impacts. The model has been validated with respect to its performance in predicting the price impacts of exogenous supply side shocks, such as those that might result from global climate change44. Additional information on the structure of GTAP is given in Supplementary Fig. 22.

GTAP captures a number of dimensions important for determining the welfare implications of climate change impacts on agriculture. These include the shifting of land area between crops, potential intensification of production, shifting of consumption between commodities and sources of goods, and the adjustment of global trade patterns (Supplementary Table 2). For the purposes of this study, GTAP is run with 140 regions and 14 commodities—with the latter designed to place an emphasis on the agricultural sector. Productivity changes are introduced to GTAP as a Hicks-neutral shift in the production function such that farmers employing the same combination of inputs would experience X% lower output in the presence of a X% climate-driven yield shock.

Wheat and rice are modeled as individual sectors within each region. Maize is part of the coarse grains sector and soybeans is part of the oilseeds sector. Impacts in these sectors are scaled downwards based on the relative importance of maize and soybeans for sectoral production in each region. Yields of crops not covered in the meta-analysis (coarse grains nec., oilseeds nec., sugarcane, cotton, and fruits and vegetables) are not altered. Absent normalization, this will lead to an underestimate of potential climate impacts, since these other sectors are also likely to be affected by climate change. Therefore, in the results that follow, welfare changes are normalized by the value of production of the crops covered in the meta-analysis. Global and regional welfare changes are measured in terms of equivalent variation and are decomposed into the three components shown in Fig. 2 following Hertel and Randhir25.

In order to explore uncertainty in the economic modeling, we perform a systematic sensitivity analysis of GTAP output to perturbations in four sets of key parameters governing the supply and demand behavior in this model. On the supply side, these pertain to the parameters determining the intensive (substitution of other inputs for land) and extensive (land supply elasticities) margins of crop supply response to commodity price. On the demand side, these are the parameters that govern the price elasticity of demand for food and the price elasticity of demand for imports—which in turn govern the price responsiveness of export demands. Parameters vary by commodity/sector. We develop symmetric, triangular distributions for each parameter value, based on estimates in the literature (Supplementary Table 8) and sample from these distributions using the Gaussian Quadrature approach implemented by Arndt45. This approach has been shown to perform nearly as well as a complete Monte Carlo analysis in the context of CGE modeling, but it is much more efficient, requiring far fewer model solutions46. Due to the computational burden of conducting a complete, systematic sensitivity analysis in the 140 region model, we collapse those regions down to the 16 FUND regions for purposes of this robustness check. The resulting mean and standard deviations for regional welfare are reported in Supplementary Fig. 13. Because on the dominance of direct effects (i.e., the impact of climate change in yields) in many of the regions’ total welfare, variation of the economic parameters has a modest impact on the underlying uncertainty.

Calculating the social cost of carbon

Results of the economic modeling are used to create damage functions that relate changes in economic welfare (measured as percentage of the value of agricultural sector output) with temperature change. GTAP results are aggregated from the country level to the 16 FUND regions. Damage functions are based on a linear interpolation between the point-estimates of welfare changes at 1, 2, and 3 °C of warming and then a linear extrapolation beyond 3 °C (results reported in main text) or on a quadratic fitted through the point estimates (Supplementary Table 5).

These agricultural damage functions are then incorporated into a sectorally and regionally disaggregated SCC damage module based on the FUND model, keeping the rest of the impact sectors unchanged28, 47. Damage functions in the module use the central parameter estimates of FUND. The full FUND model includes probability distributions over many parameters and is designed to be run in a Monte Carlo mode47. This uncertainty is not dealt with in this paper, meaning uncertainty reported in the SCC reflects only the uncertainty in the yield response derived from the meta-analysis. The damage module is driven by a standardized socio-economic and emissions pathway and climate model28. We use a business-as-usual emissions scenario (Scenario 2 in ref. 3), paired with the DICE climate module48. This produces a warming of 4 °C of warming above preindustrial by 2100. The SCC is calculated by adding a 1 Gt pulse of CO 2 emissions to this reference emissions path in 2020 and comparing the time path of damages along the perturbed pathway to the reference case. Then these incremental damages (or benefits) are discounted back to 2020 at a 3% discount rate and normalized by the CO 2 pulse volume to give the SCC. Results using alternative discount rates are given in Supplementary Table 5. As the SCC is additive, it can be decomposed by sector and region, allowing a detailed comparison of the regional impacts in agriculture between FUND and the revised regional damage functions.

Code availability

Code for GTAP is open source and available for download at http://www.gtap.agecon.purdue.edu/resources/res_display.asp?RecordID=2458. The FUND model is open source and available at http://www.fund-model.org/source-code. All other code is available from the authors upon request.

Data availability

Data for the meta-analysis is available at ag-impacts.org. Gridded yield changes based on the meta-analysis are available at 10.6084/m9.figshare.5417548. Welfare changes for GTAP regions based on yield shocks from both the meta-analysis and AgMIP are available at 10.6084/m9.figshare.5417557 and 10.6084/m9.figshare.5417560. Other data are available from the authors upon request.