The impacts of future climate and carbon dioxide changes on the average and variability of US maize yields under two emission scenarios

The United States is the largest producer of maize in the world, a crop for which demand continues to rise rapidly. Past studies have projected that climate change will negatively impact mean maize yields in this region, while at the same time increasing yield variability. However, some have questioned the accuracy of these projections because they are often based on indirect measures of soil moisture, have failed to explicitly capture the potential interactions between temperature and soil moisture availability, and often omit the beneficial effects of elevated carbon dioxide (CO2) on transpiration efficiency. Here we use a new detailed dataset on field-level yields in Iowa, Indiana, and Illinois, along with fine-resolution daily weather data and moisture reconstructions, to evaluate the combined effects of moisture and heat on maize yields in the region. Projected climate change scenarios over this region from a suite of CMIP5 models are then used to assess future impacts and the differences between two contrasting emissions scenarios (RCP 4.5 and RCP 8.5). We show that (i) statistical models which explicitly account for interactions between heat and moisture, which have not been represented in previous empirical models, lead to significant model improvement and significantly higher projected yield variability under warming and drying trends than when accounting for each factor independently; (ii) inclusion of the benefits of elevated CO2 significantly reduces impacts, particularly for yield variability; and (iii) net damages from climate change and CO2 become larger for the higher emission scenario in the latter half of the 21st century, and significantly so by the end of century.


Introduction
A growing body of research has emphasized the critical importance of high temperature days, often measured in terms of degree days above 29°or 30°C and called Extreme Degree Days (EDD), as a predictor of yields in rainfed systems (Schlenker and Roberts 2009). The main reason behind the importance of EDD appears to be that it is a good measure of cumulative evaporative demand during the growing season, and more direct measures such as average vapor pressure deficit (VPD) perform equally well (Lobell et al 2013. A warming climate in which high-EDD growing seasons occur more frequently would therefore be expected to decrease mean yields, and has also been projected to increase variability (Schlenker 2006, Urban et al 2012. Despite these findings, disentangling the effects of heat and water stress remains a challenge for understanding agricultural impacts of climatic change. Some researchers have countered that attributing severe yield losses to evaporative demand reflects omitted variables bias (Ortiz-Bobea 2013, Basso and Ritchie 2014), in which the effect of high demand in fact largely owes instead to low moisture. They argue that the predicted yield response to a measure of evaporative demand (such as VPD, EDD, or maximum daily temperature) can become less severe when also properly controlling for some measure of soil water supply (such as precipitation or model-derived soil moisture).
While it is true that omission of accurate moisture measures could potentially bias results toward overstating the importance of temperature, it is also true that failing to account for moisture could understate the impacts of warming in some settings. For example, during the vulnerable stages of crop growth that occur in roughly the third month after planting, low-supply conditions can significantly exacerbate the demand response. This is physiologically expected, because low soil moisture will impede a plant's ability to access water needed to sustain high transpiration rates. Indeed, this interaction is accounted for in the APSIM Maize model, which calculates the ratio of moisture supply to demand when determining crop growth rates (Holzworth et al 2014). Despite strong evidence for the coupling of hot days and low soil moisture (Mueller and Seneviratne 2012), as well as crop simulations highlighting the role of moisture stress in yield variability (Andresen et al 2001), empirical models have largely ignored this interaction, with the possible consequence of underestimating impacts in future climate scenarios with simultaneous warming and drying.
The interaction between temperature and moisture is further influenced by atmospheric carbon dioxide (CO 2 ) levels. In C4 crops such as maize, it is well established that intercellular CO 2 concentration is already saturated, and thus maize exhibits little to no direct stimulation of photosynthesis at higher ambient CO 2 under well-watered conditions (Samarakoon andGifford 1996, Leakey et al 2009). However, an indirect photosynthetic improvement arises through decreased stomatal conductance in the presence of water stress (Ghannoum et al 2000, Leakey et al 2006. This improvement in transpiration efficiency with elevated CO 2 leads to positive effects under dry conditions. The ability to account for this effect is a major need in empirical approaches, particularly in cases such as rainfed maize where the main mechanism behind empirical temperature effects is related to water stress. In this study, we first use a new dataset on maize yields in the central US to investigate the effect of simultaneous high evaporative demand (hereafter referred to as 'demand') and low soil moisture supply ('supply') during the period roughly two to three months after planting, a period during which the crop is especially sensitive to this damaging combination (Çakir 2004). We then examine the implications of potential changes in temperature and moisture for maize yields under alternate maize model specifications, including ones that ignore measures of supply and/or interactions between supply and demand to gauge potential errors in previous approaches. We also introduce a new approach for accounting for CO 2 effects based on reducing effective demand. Finally, we estimate how differing levels of CO 2 and climate changes according to different RCPs could influence maize yields in this region, both in terms of mean and variability of yields. Our work advances the current literature by empirically demonstrating robust evidence for an interaction between moisture and temperature in a large US maize dataset, quantifying its implications for yield variability, and the potential of elevated CO 2 to mitigate this effect through improved water use efficiency.

Methods
We use yield data from the Risk Management Agency containing 100 randomly selected fields per countyyear for the years 1995-2012 in Iowa, Illinois, and Indiana (Lobell et al 2014). We fit regressions to both these and county-averaged values of the same data. Having found that both fieldand county-level regressions give similar responses, we used the county-level panels for computational tractability of spline fitting and bootstrapping. For historical weather, we use 30 day average PRISM-derived precipitation, VPD, and temperature values, downscaled to 800 m. Each 30 day interval is taken relative to each field's reported sow date (or average sow date in the case of the countylevel data). We consider the interaction of supply and demand specifically in the 61-90 day period following sowing. This period,which roughly corresponds to July and the main window of flowering in this region, is the most sensitive growth stage to the combination of soil moisture supply and evaporative demand. We also consider VIC-derived 1/8°soil moisture data from the NLDAS-2 dataset (Xia et al 2012) as a supply measure, but find that precipitation (P) gives slightly better model skill, explaining about 5% more of the variance in residuals from a model without weather. For this reason, as well as the fact that precipitation data are generally more available than soil moisture data, we use P as our supply variable in the models and results that follow.
In order to identify the interaction between supply and demand in these data, we potentially require a more flexible fitting procedure than linear interaction terms within an ordinary least squares regression. For example, the yield sensitivity to demand might depend on supply, but only at low supply levels. Such a relationship would not be captured by a linear interaction model. To search for this interaction's functional form, we used generalized additive models (Wood 2006a(Wood , 2006b, which include tensor splinebased smoothing terms, together with linearly additive terms for other (non-interacting) covariates. While overfitting can be a liability of such a flexible model, we found that under sufficient smoothness penalization, the response surfaces were stable under rigorous cross validation and bootstrapping robustness checks. Furthermore, the most skillful spline-based fit could be well approximated by an OLS model in which both linear and quadratic demand terms interact linearly with supply (see results).
In order to ensure our results did not depend on an overly specific choice of interaction form, we considered a broad suite of regression models that are summarized in table 1. In all models, if year or supply does not appear in an interaction, it is still included as a linear control term, as are sowing date, 91-120 day post-sowdate maximum temperature, and a categorical county identifier. These controls serve to significantly (p < .001) increase variance explained, but importantly to our purposes, they ensure that the supply-demand effects are not confounded by late sowing or high heat during grain filling, both of which could be correlated with supply and demand in our target window. Our preferred demand measure is average maximum VPD, as it has been shown to be a more effective predictor than temperature and has a more direct physiological connection to the water demand experienced by the crop (Lobell et al 2013. We consider both spline-based GAM models and their OLS analogs that best approximate the spline-based smoothing terms, and within each of these two classes we consider increasing degrees of interaction, starting with models containing no interactions and ending in those in which year, supply, and demand all interact jointly. Using curvilinear relationships between yield and weather variables such as T and P in a regression model has a long history in the empirical maize modeling literature (Thompson 1969). Such studies have often used monthly or seasonal averages of T and P, but in more recent years as more extensive daily weather datasets have become available, daily measures such as growing degree and EDD have been used to detect nonlinear relationships between these variables and yields (Roberts et al 2013). A common criticism of empirical models is their lack of accounting for the timing of weather events, but recent work has leveraged new datasets to compute weather measures for sowdate-specific intervals, with significant model improvement over monthly data (Lobell et al 2014).
Our study belongs to this latest category, but with the additional step of quantifying the impact of interacting heat and water stress during such an interval. While previous models explain an increasingly large amount of variance and are useful for certain projections, they are likely to underestimate yield loss under the particular combination of extremely high demand and low supply.
To check the robustness of the observed interaction effects, we apply a bootstrapping procedure by resampling only years, thereby treating the entire spatial extent as a single block. This is a conservative approach in that the number of clusters in our data without correlated errors is likely much higher than the number of years (18). In effect we are controlling for spatial autocorrelation over the entire Corn Belt region, when in reality any weather-related autocorrelation is likely to be smaller than the entire region. We cross validate both by random subsets (10fold) and by year, in which we withhold all observations belonging to the same year as the test set, and train on the remaining 17 years. The by-year cross validation serves as a check against models that overfit based on high-leverage years, a particular concern in our data which includes the hot, dry, and low-yielding year of 2012.
To estimate the plausible range of supply-demand interaction effects under future climates, we use the output of 29 CMIP5 models run under low and high emission scenarios of RCPs 4.5 and 8.5, respectively. We compute the supply and demand variables' relative changes between 2010 and 2050, 2070, and 2090 for each GCM and emission scenario, and apply those changes to the historical data. As a benchmark, we also consider the hypothetical scenario of supply decreasing uniformly by one standard deviation of its historical distribution, and demand increasing by one standard deviation.
Because future changes in the supply-demand relationship could also be modulated by the improved water use efficiency of elevated CO 2 , we scale the differences between current and future VPD by the percent increase of future CO 2 from a 2010 baseline of 389 ppm. The scaling factor used by the APSIM maize model to modify transpiration efficiency and effective water demand increases linearly at 10.6%/100 ppm beginning at 350 ppm (Harrison et al 2014). Here we use it to artificially reduce the VPD seen by the regression model. For example, if projected VPD increases Table 1. Each row represents a different level of interaction between year, supply (s), and demand (d), and each row contains two model specifications-a GAM model (top) with tensor spline smooths indicated by the te() terms, and the OLS model (bottom) that most closely approximates the GAM model. In all cases, 'controls' refers to linear (i.e. non-interacting) additions of planting date, maximum daily temperature in the 91-120 day post-planting period, and a categorical county identifier. In the OLS models, the '*' symbol represents an interaction including all lower order terms.

Results and discussion
We find a stark interaction between supply and demand in these data, as illustrated in figure 1. The leftmost panel shows the joint response of yield to both supply and demand when fitting a highly flexible tensor spline-based model (abbreviated as te(s,d) in table 1). High evaporative demand is more damaging at low levels of rainfall, consistent with physiological expectations that crops are more sensitive to heat when soil moisture is inadequate to sustain high rates of transpiration. When constraining the fit to further penalize unevenness in the surface (center panel), we achieve a more stable result that still exhibits a steeper yield decline with VPD at low P than with average P. We can approximate this fit very well by replacing the spline-based smoothers with a P*(VPD + VPD 2 ) interaction in an OLS model (right panel, abbreviated as s*d2 in table 1). We find this interaction to be significant under stringent robustness checks. ANOVA comparisons of models with and without interactions show the interaction model to be highly significant (p < .001), and this significance holds up under a much more conservative test of block-bootstrapping the data (figure 2). We conclude that this interaction is a robust feature in these data. We also find some evidence for an increasing yield sensitivity to demand over time, as captured models in which supply and demand also interact with year, but this is considerably less robust than the interaction between supply and demand only, and we do not consider further in this study.
The GAM models, while explaining more variance and performing better under random-sampling cross validation, become too sensitive to extremely high demand values. Cross validating by year similarly reveals that some of these higher-flexibility models are too sensitive to extreme years like 2012 (figure 3). Despite this, our range of regression models shows good agreement on the magnitude of yield effects in response to demand increases and supply decreases. The models that included an interaction between supply and demand, but not year, struck the best balance between high variance explained and low test error. These are the models in the second row of table 1, and for the sake of simplicity, we hereafter present the The flexibility of the spline fits with their marginal bases lightly penalized (left) identify a strong interaction in the low-P-high-VPD corner, but are prone to overfitting. By increasing the smoothness penalty on the marginal fit of both P and VPD (center), the interaction remains and the tendency to overfit is greatly reduced. A traditional least-squares model with an interaction between P and (VPD + VPD 2 ) shows a very similar fit. Figure 2. The interaction of precipitation and VPD during the 61-90 day period following planting has a large impact on final yield. While few data and hence large uncertainty exist for high supply high demand conditions, the yield response to demand at very low supply (2nd percentile precipitation, blue curve) is significantly lower than the response at median precipitation levels (red curve). Error bars represent the 5th and 95th bootstrap confidence quantiles.
results of the OLS model 's*d2' as described in table 1. We proceed with this choice for the easier interpretability of its functional form over that of a tensor spline smooth, but emphasize that the shape of the GAM model's interaction motivated this OLS formulation, and that both the s*d2 and te(s,d) models give nearly identical results. The full OLS model specification is then: where supply and demand are P and VPD in the 61-90 day post-sowing period, T max is the average maximum daily temperature in the 91-120 day post-sowing period, and the interaction expression contains all lower order terms (e.g. supply*demand). Our preferred supply variable is P rather than VICderived soil moisture, though both are imperfect proxies for true soil water supply. The VIC model does not account for certain factors affecting soil moisture such as year to year phenology changes due to crop rotation, or management decisions such as tiling and tilling. It does, however, represent spatial differences in soil texture (albeit crudely for an agricultural setting), and accounts for soil moisture memory, both of which are important and neither of which is accounted for by precipitation, and has been shown to reproduce observed inter-annual and intra-seasonal variability and persistence of soil moisture across the US (Xia et al 2014). However, the shortcomings of the VIC model appear to outweigh its benefits in this particular application, as using P as our supply measure results in a modest but significant improvement in model skill (R2 of .67 versus 65). Nonetheless, soil properties, especially as they relate to retention, storage capacity, matric potential and runoff, are important to this topic, and improved datasets with large spatial and temporal coverage could provide significant advantages in future empirical research of crop water and yield relationships.
Having established a robust interaction between P and VPD, we turn to the question of whether this effect has meaningful bearing on yields. As seen in figure 2, the values of VPD for which the yield response significantly differs between high and low P represents a small regime of the data. It could be that the interaction, while statistically significant, has little impact, owing to the relatively small mass of data it represents. We find, however, that this is not the case. Figure 4 shows the county-level effects of a hypothetical increase of VPD and decrease in P by one standard deviation from their 1995-2012 levels. We see that the response to this level of warming and drying differs by as much as 6% for mean losses and 25% for interannual standard deviation increases between a model in which P and VPD explicitly interact (right panels) and one in which they do not (left panels). We attribute the comparatively larger effect on variability than mean to the form of the interaction identified in figure 1. Increasing VPD and decreasing P pushes more cases toward the low-P-high-VPD corner of the interaction surface. The majority of observations still lie in the flat region of the surface (such that the mean is only modestly affected), but the addition of cases to steeper portions of the surface creates the potential for some county-years to experience dramatic yield losses, with consequent increase in interannual variability. It is therefore likely that many previous studies, even those including moisture controls, have likely underestimated yield variability in response to high temperatures.
Finally, we consider the yield effects in response to the changes between the historical baseline and future climate projections of supply and demand measures. Of particular interest here is whether the stark supplydemand interaction we find could in fact be substantially mitigated by higher CO 2 , which effectively scales down demand through improved water use efficiency. It is possible that higher emissions scenarios lead to hotter and drier conditions, but that the effective demand owing to CO 2 mitigation results in effective demand similar to lower-emission scenarios, with consequently similar yield effects. We address this question by considering the output of 29 CMIP5 ensemble members, the projected changes of which from a 2010 baseline are seen in figure 5. We apply these changes to the 1995-2012 historical values, once when factoring in the CO 2 scaling of demand, and once without. In the former case, we reduce the projected VPD changes by the percent improvement in transpiration efficiency per CO 2 increase as prescribed by the APSIM maize model (see methods), while in the latter we apply no changes to the projected VPD increases. In both emission scenarios, climate trends toward higher evaporative demand push yields toward  table 1 for a more complete description of each model). In the 'random' group, we perform 10-fold cross validation on a random shuffling of the observations, while for 'year', each test set consists of one year. The models with supply interactions in general perform slightly better than those without with regard to both 10-fold and by-year test error, but the spline fits including interactions with year tend to overfit (two rightmost bars).  greater mean losses and greater interannual variability. These trends are more severe in scenario 8.5, but its higher CO 2 levels mean comparatively greater mitigation of these yield impacts than scenario 4.5. A key question, therefore, is whether the greater climaterelated stress of the higher emission scenario outweighs the alleviation potential of higher CO 2 , or vice versa. We illustrate the results in figure 6, which shows the distributions of yield mean and variability changes aggregated over all counties in response to rising emissions over the 21st century under the RCP 4.5 and 8.5 scenarios. By 2050, the large CO 2 benefits of RCP 8.5 counterbalance its higher demand, to the point that its yield effects are not significantly different than RCP 4.5. By 2070, however, the demand increases of the high emissions scenario outpace the CO 2 gains for most models, and by 2090, 90 to 95% of models show 8.5 as the less favorable scenario.
To further examine the dependence of impacts on emission scenario, figure 7 illustrates the distribution of differences between model projections for the scenarios 8.5 and 4.5. The inter-scenario differences increase over time, with 8.5 showing greater median changes in yield mean and CV by 2070 than 4.5. By 2090, these differences are much greater, and significant at the 5% level for yield mean changes and at 10% for yield CV changes. This is true both when demand is scaled down by CO 2 (blue bars) and when it is not (red bars). This indicates that while CO 2 has potential to mitigate both mean and variability changes in yields, the higher emissions scenario nonetheless becomes less favorable by mid-century, and significantly so by end of century.
We repeat the above analysis using VIC-derived soil moisture as the supply measure, and estimate future changes with CMIP5 top-layer soil moisture projections, obtaining close agreement to the Figure 6. The yield impacts under lower (RCP 4.5) and higher (RCP 8.5) emissions scenarios are small in the mid 21st century, but become substantial by the end of the century with respect to both mean and variability changes. While smaller yield effects result from the reduced effective demand of higher CO 2 (light-colored bars), higher emissions result in demand-related damage that outweigh the CO 2 benefits, and the distinction between emissions scenarios remains significant. Horizontal dark bars indicate the medians across 29 climate models while crosshairs indicate means. Box widths are defined by the interquartile range, while lower and upper whiskers indicate the 5th and 9th percentiles. Figure 7. The differences between scenarios 8.5 and 4.5 with respect to yield mean and variability increase over time. The dashed line at 0% indicates no difference between scenarios. For mean changes, values below the dashed line indicate models for which mean yield losses in 8.5 are greater than 4.5, while for CV, values above the line indicate models for which yield variability increases in 8.5 are greater than in 4.5. Horizontal dark bars indicate the medians across 29 climate models. Box widths are defined by the interquartile range, while lower and upper whiskers indicate the 5th and 9th percentiles.
precipitation-based results. We also examine the robustness of our results to extrapolation beyond the regime of the training data. In most scenarios, between 92 and 99% of the projected supply and demand distributions lie within the range of the historical distribution. In the most extreme case of RCP 8.5 in 2090, however, nearly 24% of the increased demand lies above the historical maximum when not considering CO 2 effects, and 23% lies below the historical minimum when scaling according to CO 2 increases. However, we find little difference in the model average mean and CV changes when removing these observations from the analysis, and our conclusions regarding the inter-scenario differences remain unaffected.
While noting that management factors can also affect yield sensitivity to moisture stress (Basso and Ritchie 2012), FACE studies show that yield improvements under elevated CO 2 are below what might be expected from the increased transpiration efficiency (Leakey et al 2009). Thus, our model's yield benefit via reduced demand could overstate the true CO 2 yield benefit, with the consequence of understating the lower emission scenario's advantage.
As with any statistical approach, we caution that our coefficients could be affected by measurement error in the predictors, and that our results are largely restricted to the domain of our data. Though suggestive, they do not immediately extend to other regions, crops, or supply and demand regimes outside those used to fit our model. Future statistical (as well as simulation-based) studies that expand upon the topics presented here would benefit from a deeper physiological understanding of the nonlinear responses arising from interactions between evaporative demand, soil moisture supply, and CO 2 .

Conclusions
The interaction between soil water supply and evaporative demand during the most sensitive period of crop development is a robust phenomenon that has not previously been adequately represented in empirical models. We find that this interaction leads to larger mean yield losses, and larger still yield variability increases, when conditions become both hotter and drier. The CMIP5 ensemble suggests a higher frequency of years in which crops will face both high evaporative demand and water-limited conditions, but we find that under higher emissions scenarios, the demand-reducing effect of elevated CO 2 can substantially reduce these changes. However, by end of century, a lower emissions scenario with lower associated crop water demand becomes distinctly favorable.