Crop-climate feedbacks boost US maize and soy yields

US maize and soy production have increased rapidly since the mid-20th century. While global warming has raised temperatures in most regions over this time period, trends in extreme heat have been smaller over US croplands, reducing crop-damaging high temperatures and benefiting maize and soy yields. Here we show that agricultural intensification has created a crop-climate feedback in which increased crop production cools local climate, further raising crop yields. We find that maize and soy production trends have driven cooling effects approximately as large as greenhouse gas induced warming trends in extreme heat over the central US and substantially reduced them over the southern US, benefiting crops in all regions. This reduced warming has boosted maize and soy yields by 3.3 (2.7–3.9; 13.7%–20.0%) and 0.6 (0.4–0.7; 7.5%–13.7%) bu/ac/decade, respectively, between 1981 and 2019. Our results suggest that if maize and soy production growth were to stagnate, the ability of the crop-climate feedback to mask warming would fade, exposing US crops to more harmful heat extremes.

Despite global warming, damaging heat has not increased over the northern Midwest, home to the most productive US maize and soy croplands [6,[22][23][24]. KDD trends inversely covary with yield trends over maize-and soy-growing regions in the US, with areas of strong yield growth having less warming (or cooling) of hot extremes [22,23] (cf figures 1(a)-(f); SI figure 15). Reduced warming trends over croplands have been shown to benefit yields and to be most pronounced on the hottest days of the growing season [6,24,25]. At the same time, GDDs have increased across the entire US (figures 1(c) and (d) and SI figure 2), lengthening the growing season, benefiting maize and soy. Thus, the effects of climate change have been unexpectedly beneficial in recent decades: increased warmth (measured by GDDs) without an increase in damaging heat [6] (measured by KDDs). Crops affect local climate by modifying soil moisture, albedo, and turbulent fluxes, which shape the surface energy budget [26]. Increased crop growthwhether due to rising yields or expanding cropped area-increases evapotranspiration (ET) given a constant water use efficiency, a stable growing season, and sufficient moisture availability. This ET growth increases latent heating and reduces sensible heating, cooling local air temperatures, a process which has been demonstrated in modeling studies [27][28][29][30][31].
However, the extent to which US crops are influenced by climate trends arising from agricultural intensification itself rather than other climate drivers like aerosols, greenhouse gasses, or oceanatmosphere variability [32], is uncertain. Here we show that maize and soy production trends can explain nearly all of the local cooling across the US Midwest. We demonstrate the magnitude and spatial variability of a crop-climate feedback, a positive (negative) feedback whereby increased (decreased) crop production enhances (reduces) latent heating, reduces (increases) sensible heating, and reduces (increases) damaging heat, further increasing (reducing) crop production. This feedback has occurred as the result of rising yields driven by technological change over predominantly rainfed croplands, despite little growth in harvested area (see SI figure 1(b)). Our quantification of this feedback allows us to assess how agricultural intensification has interacted with climate variability and change to shape yields over the past four decades in the US. Most importantly, our methodology provides an understanding of whether recent crop-favorable climate trends should be expected to continue, or whether potential future yield stagnation [1,33] or moisture limitation could result in diminished local cooling, exposing the background warming trend and resulting in production losses begetting further production losses.

Data and methods
We calculate growing season mean precipitation, sensible heat flux, latent heat flux, net radiation, and wind speed by extracting monthly mean values from the ERA5-Land [34] reanalysis for months within the growing season as defined by a location-specific crop calendar for the primary maize growing season [35] (growing seasons for soy are very similar to those of maize, see SI figure 12). Multiple cropping is not common in the US Midwest [36]. We average monthly values across the local growing season. We also extract the same data from the Global Land Data Assimilation System (GLDAS)-V2-Noah and GLDAS-V2-Vic products [37], average them together, and repeat the analysis as a robustness check of our results (SI figure 4). We calculate GDDs and KDDs using daily maximum and daily minimum 2 m temperature data from ERA5-Land. While ERA5-Land does not represent land cover change, it does include the effects of agricultural intensification by assimilating temperature, humidity, and other observations that are influenced by crops. In our robustness check using GLDAS, we calculate GDDs and KDDs using daily temperature data from the Climate Prediction Center (CPC) [38]. GDDs and KDDs are calculated using the following formulas: Other definitions of GDDs and KDDs are used in agronomic literature, but the definitions we use here have been used in recent crop-climate studies [6]. GDDs and KDDs are not intrinsic drivers of crop growth, but instead capture sunlight, instantaneous temperature exposure, and moisture availability, and are predictive of crop yields [6]. We obtain countyscale, annual maize and soy yields, along with their harvested areas, for 1981-2019 from the USDA NASS QuickStats [39] tool. The gridded climate data is upscaled to county resolution by spatially averaging grid cells falling within each county. We calculate total maize and soy production using the following formula: where Y is county yield and HA is county harvested area. We calculate linear trends for climate and crop quantities using an ordinary least-squares fit. We calculate the fraction of each county's maize and soy that are irrigated using an irrigated area map provided by the Food and Agriculture Organization [40]: Irr frac = Irr area HA maize + HA soy .
For each county, we linearly detrend production and climate data and subtract each variable's mean value.
We also repeat the analysis using singular spectrum analysis to detrend the production and climate data and find no change in the results: estimates of the crop-climate feedback's effect on maize and soy yields are not significantly different (K-S test, p = 0.99) as compared to results using linear detrending. Between 1981 and 2019, maize and soy production has grown linearly at a rate of 113 M bu yr −1 (see SI figure 11).
We estimate how latent heating varies with maize and soy production using the following regression in each county: where maize and soy production are detrended anomalies, and LH, Pr, R net , and W are detrended anomalies of growing season mean latent heating, precipitation, net radiation, and surface wind speed, respectively (see SI figures 3(a)-(d) for regression coefficient magnitudes). This equation considers the primary factors controlling ET and thus latent heating [41]. Across counties, the median R 2 is 0.78, indicating that these predictors account for most of the interannual variance in LH. Because precipitation and KDDs are negatively correlated (r = −0.6), the precipitation term in the regression partially captures the response of LH to high temperatures. We exclude counties which fail to meet any of the following criteria: Across the US, 398 counties meet these criteria, most in the central US where maize and soy cultivation is most intense, and 375 of these counties have significant (p < 0.05) coefficients on production. We then exclude counties with non-significant coefficients on production. Included counties account for 59% of national maize and soy production, and the median interannual Pearson correlation coefficient between maize and soy yields is 0.69. Production is minimally colinear with climate variables: the absolute value of Pearson correlation coefficients between production and precipitation, net radiation, and wind are all less than 0.2. In counties with low harvested fractions, β Production becomes relatively large (>1 W m −2 per million bu), potentially because of the influence of non-crop vegetation on latent heating. However, maize and soy production is highly concentrated, with 80% of production occurring in counties with more than 50% of their area harvested, meaning that the counties that contribute most to crop production have the least non-crop vegetation. We find that crop-driven LH declines with harvested fraction, showing that our regressions are unlikely to be strongly influenced by non-crop vegetation (SI figure 5).
Using the fitted coefficients from equation (1), we estimate the time series of growing season latent heating anomalies produced by recorded production anomalies, denoted LH ag-int . For each county, we insert the time series of maize and soy production anomalies (relative to the intercept after fitting a linear trend to county production) into the regression. Unlike the production data used to fit the regression, this production time series is not detrended to capture agricultural intensification over time. We control for variability in precipitation, net radiation, and wind by setting those anomalies to zero in the regression. Because production has an upward trend and β Production is positive in most counties, LH ag-int also has a positive trend in most counties ranging from 0.5 to 3 W m −2 per decade.
To relate latent heating variability to sensible heating (SH) variability, we regress detrended growing season mean SH anomalies against detrended growing season mean LH and R net anomalies between 1981 and 2019 (see SI figures 3(e)-(f) for regression coefficient magnitudes): LH and R net almost entirely predict SH: the median R 2 value across counties is 0.99. Partitioning R net into SH and LH using equation (2) is equivalent to calculating the Bowen ratio ( SH LH ) and using this to partition R net into the turbulent fluxes. This partitioning is done for each county separately to take into account variations in energy partitioning due to albedo, cropped area, natural vegetation type, and other factors. We use this regression to estimate the sensible heating anomaly from maize and soy production, SH ag-int , by inserting the LH ag-int time series into the regression and setting R net to zero. In all counties, the estimated SH ag-int time series is nearly equal to LH ag-int , but with an opposite sign: trends range from −0.5 to −3 W m −2 .
Next, we fit regressions to estimate how detrended growing season mean GDD and KDD anomalies vary with detrended growing season mean SH anomalies for each county (see SI figures 3(g) and (h) for regression coefficient magnitudes): Coefficient magnitudes are positive and significant in all counties for both KDDs and GDDs, indicating that growing seasons with higher mean SH have more KDDs and GDDs. Because KDDs and GDDs are daily-scale temperature metrics, however, they are only partially correlated with growing season scale SH: we find median R 2 values across counties of 0.56 and 0.28 for KDDs and GDDs, respectively. We use regression to estimate SH effects on GDDs and KDDs, rather than the exact GDD and KDD functions, because the seasonal-scale effects of production on SH cannot be disaggregated to a daily scale. Then, we estimate the time series of KDD and GDD anomalies resulting from the time series of production-related SH anomalies by inserting SH ag-int into equations (3) and (4). We denote these time series as KDD ag-int and GDD ag-int , and we find that they have negative trends in most counties, indicating that agricultural intensification has cooled local temperatures. We estimate the effects of climate trends on maize and soy yields using a statistical yield model fit over 1981-2019 detrended data for growing season total GDD, KDD, and mean precipitation for each county: We fit this model independently for maize and soy (see SI figures 3(i)-(n) for regression coefficient magnitudes). For both maize and soy, β GDD is positive, β KDD is negative, and β Pr is near zero in most counties. GDD, KDD, and precipitation coefficients are significant (p < 0.05) in most counties. The near-zero precipitation coefficients are consistent with other published results and likely reflect the strong correlation between precipitation and temperature and potentially also the differing impact of rainfall intensity on yields [42]. The regression R 2 and the coefficient on KDDs do not substantially change when using a quadratic precipitation term. The median R 2 value across counties is 0.49 for maize and 0.37 for soy, which are similar to other studies [43] and indicate that GDDs, KDDs, and precipitation substantially predict interannual yield variability. We assess uncertainty by bootstrapping each regression estimate by sampling the input data with replacement 100 times. We also sample each coefficient estimate 100 times from a student's tdistribution defined by the mean coefficient estimate and its standard error, giving 10 000 estimates of the feedback. We propagate uncertainty via the following steps: (a) We bootstrap each regression estimate 100 times by sampling the input data with replacement. (b) For each regression, we fit a student's tdistribution with mean equal to the coefficient estimate and standard deviation equal to the coefficient standard error. (c) We select 100 random values from this coefficient distribution. (d) Using each randomly selected value, we estimate the production-driven LH (or SH, GDD, KDD, or yield for later regressions). (e) For the following regression, we randomly select a coefficient for the current step, and randomly select one of the coefficients from the previous step and use that estimate of LH, SH, GDD, or KDD.
This method produces 10 000 estimates of the crop-climate feedback (100 from bootstrapping times 100 from random coefficient selection). It both samples uncertainty in the regression estimates and propagates uncertainty through the regression chain. All results are presented as the mean value and the 5th-95th percentile across the 10 000 estimates.
We use equation (5) to estimate the effects of observed GDD (GDD obs ), KDD (KDD obs ), and precipitation trends on maize and soy yields between 1981 and 2019. These climate-driven yield trends are denoted as Yield obs . Next, we estimate GDD and KDD trends in a counterfactual world with no agricultural intensification: We then insert KDD no-ag-int , GDD no-ag-int , and observed precipitation trends into equation (5) for maize and soy to estimate the impacts of climate trends on maize and soy yields in the counterfactual world with no agricultural intensification. These counterfactual yield impacts are denoted as Yield no-ag-int . Using these results, the yield impact from climate altered by agricultural intensification is calculated as: The agricultural-intensification component of the yield trend is shown as a fraction of the total yield trend for each county: where ∆t is the 1981-2019 period. The additional production associated with agricultural intensification is calculated as: Production ag-int = Yield maize-ag-int × HA maize + Yield soy-ag-int × HA soy where HA is the mean county harvested area. We use bushels per acre (bu/ac) as our unit of yield and production measurement. A bu/ac of maize is equivalent to 63 kg ha −1 and a bu/ac of soy is equivalent to 67 kg ha −1 . While a regression framework such as this can only show association and not causation, each regression step is grounded in physical mechanisms: our methodology and conclusions draw on both correlation and physical inference.

Results and discussion
Analysis of observations for 375 US counties shows that rising maize and soy production has increased regional ET, which has modified the surface energy budget by increasing latent heating (LH) and reducing sensible heating (SH) (see figures 2(a)-(d)). Reduced SH from crop production growth has led to fewer KDDs, which has benefitted maize and soy yields (figures 2(e)-(k)). These results come from a set of empirical estimates at each stage of the crop-climate feedback (figure 2; see SI figure 3 for coefficient magnitudes). We estimate the effect of interannual maize and soy production variability on the latent heat flux at county scale by regressing total maize and soy production on growing season mean latent heating, precipitation, net surface radiation, and surface wind speed (see section 2, SI figures 1-3). While weather and disease drive short term yield variation and technological and climate change drive long term yield trends, higher crop production produces more ET, allowing us to estimate production-LH relationships from interannual variability. We present results using the ERA5-Land [34] reanalysis, although our findings are robust when alternative datasets are used (see SI figure 4). Results are also robust to data detrending methodologies (see section 2). We focus on maize and soy because they are the predominant crops grown in the central US, have similar growing seasons, and have yields that are highly correlated in space and time (see section 2). While rates of ET for maize and soy are different for biophysical reasons [44][45][46], we analyze the sum of maize and soy ET together in our regression framework because it is this total latent heat flux that influences the surface energy budget. More discussion of this analytical choice is in extended data table 2.
Maize and soy production has a positive (mean coefficient of 0.77 W m −2 million −1 bu; standardized coefficient shown in figure 2(a)) and significant (p < 0.05) effect on LH in 90% of counties that Coefficients are estimated for each included maize-and soy-growing county separately. Coefficients that are significant (p < 0.05) in at least 75% of included counties (see section 2) are denoted with three asterisks. Error bars show one standard deviation in coefficient magnitude across all counties. Green bars denote quantities shown to be modified by production variability [27,31] (figures (b)-(k)). Changes in LH, SH, KDD, GDD, and yield, respectively, associated with agricultural intensification as estimated by equations (1)- (5). LH, SH, KDD, GDD were evaluated over the growing season, and all counties that met inclusion criteria are shown (see section 2). meet our inclusion criteria (see section 2; figure 2, SI figure 4(a)), indicating that years with higher levels of production have more growing season mean latent heating. We restrict our analysis to counties with significant production coefficients. Counties with a positive and significant coefficient on production account for 60% of national maize and soy production. Counties that irrigate at least 10% of their harvested maize and soy (14% of all included counties) have a mean production coefficient of 0.89 W m −2 million −1 bu, higher than the mean across all included counties. The estimated coefficients on precipitation and net radiation are also positive and significant in nearly all counties, as both moisture and energy availability strongly contribute to the latent heat flux. The estimated LH from production increases with county harvested fraction (SI figure 5(a)).
Using county-specific regressions, we estimate the LH trend resulting from each county's production trend from 1981 to 2019 (figure 2(b); see section 2). Rising production has increased LH in almost all included counties by a mean of 2.2 (5th-95th percentile: 0.8-3.6) W m −2 decade −1 , with the largest increases occurring in counties with strong production growth. This production-driven LH trend is large compared to total LH trends of −1 to 2 W m −2 decade −1 across the central US (SI figure 2(c)), representing the substantial influence that crops have on the surface energy balance in intensely cultivated regions.
An increasing share of incoming radiation partitioned into LH by growing crops should reduce the remaining energy available for SH. We estimate how this production-driven LH trend has modified SH (figures 2(c) and (d)) by regressing SH against LH and net radiation at the county scale (figure 2(c)). We find mean production-driven reductions in SH across counties of −2.2 (5th-95th percentile: −2.3 to −2.0) W m −2 decade −1 , though many counties still have overall upward SH trends driven by rising net radiation associated with global warming (see SI figure 2(d)). Reduced SH is associated with reduced temperatures and therefore decreased GDDs and KDDs (see SI figure 6), though we find that KDDs vary ∼5 times more strongly with SH than do GDDs (figure 2(e), standardized coefficients of 1.5 for KDDs, 0.3 for GDDs). Propagating these production-based changes in energy partitioning into heat accumulation indices leads to mean reductions in KDD and GDD trends of −15.5 (5th-95th percentile: −28.8% to −2.7%) KDDs/season/decade (−25.0% to −2.3%/decade) and −23.2 (5th-95th percentile: −58.6-11.1) GDDs/season/decade (−3% to 0.5%/decade) (figures 2(f) and (g)), respectively, counteracting greenhouse gas warming. We find that these effects of agricultural intensification on the local surface energy balance have produced a yield bonus across much of the US. Slower increases (or even decreases) in KDDs benefit maize and soy yields, while slower increases in GDDs harm them (figures 2(h)-(j)). However, the KDD trends dominate due to the high sensitivity of yields to extreme heat, and we estimate that production-driven changes in KDD and GDD trends have boosted maize and soy yield trends in nearly all counties by a national mean of 3.3 (2.7-3.9) and 0.6 (0.4-0.7) bu/ac/decade (figures 2(i)-(k)), respectively. These results do not substantially differ for only rainfed counties, indicating that the crop-climate feedback is not a consequence of irrigation [30,47] (see SI figure 8).
Without the cooling driven by maize and soy production growth, KDDs would have increased everywhere, including in the Upper Midwest, and GDDs would have risen by a larger amount than was observed (see SI figure 7). Using our estimate of the crop-climate feedback, we show that if agricultural intensification had not occurred since 1981, climate trends would have caused a decline in US maize and soy yields. We compare the yield trends caused by observed climate trends (figures 3(b)-(e)); our estimates of recent climate-driven maize yield trends are similar to other studies [6,25] since 1981 with an estimated counterfactual world in which climate was not altered by agricultural intensification (figures 3(c)-(f)). Observed climate trends changed national mean maize and soy yields by 0.2 (5th-95th percentile: −0.1 to 0.4) and 0.4 (5th-95th percentile: 0.3-0.4) bu/ac/decade, respectively, since 1981 (figures 3(b)-(e)). Without the climate benefits of agricultural intensification, we estimate that additional climate warming would have instead changed national mean maize and soy yields by The crop-climate feedback has buffered the harmful effects of greenhouse gas warming and enhanced national mean maize and soy yields since 1981 (figures 4(a) and (b)). We estimate that observed climate trends have had a slightly positive effect on US maize (mean change of 0.4%, relative to mean 1981-2019 yields) and soy (mean change of 3.1%, relative to mean 1981-2019 yields) yields. However, without agricultural intensification, we estimate that the national mean effect of climate trends on maize and soy yields would have been −7.3% and −0.7%, respectively. Aggregated over the 1981-2019 period, the yield benefit from agricultural intensification has boosted total US maize production in our included counties by approximately 514 M bu, comparable to Canada's annual maize production. Similarly, soy production has been boosted by approximately 76 M bu, more than South Africa's annual production.
The positive climate effects of agricultural intensification on yields have occurred in nearly every maize and soy-growing county, with the size of the effect dependent on the magnitude of production-driven KDD change (figures 4(c) and (d)). This benefit has occurred even in counties with rising overall KDD trends: without the reduction in KDDs from agricultural intensification, the KDD increase would have been greater. The counties that have had the most yield damage from recent climate trends have had the greatest relative benefit from the crop-climate feedback (figures 4(c) and (d)); red dashed lines have slopes less than one, indicating that a larger fraction of recent production trends can be attributed to the crop-climate feedback in counties with more negative climate-driven yield change. In 41% of maizegrowing counties and 29% of soy-growing counties, the crop-climate feedback reversed the sign of the climate-driven yield effect from negative (reduced yields) to positive (increased yields). We find that the crop-climate feedback contributed a national mean 16.8% (5th-95th percentile: 13.7%-20.1%) of total county-level maize yield trends between 1981 and 2019 ( figure 4(e)). For soy, the crop-climate feedback contributed a national mean of 10.5% (5th-95th percentile: 7.5%-13.7%) of total county-level yield trends ( figure 4(f)). The contribution of the cropclimate feedback to yields is dependent on non climate-driven yield trends, the differential sensitivity of crops to GDDs and KDDs, maize and soy cultivars, fertilizer use, and other farming practices [5].
Our results suggest that the crop-climate feedback provides the most benefit to regions with the most concentrated crop production (SI figure 5(b)). We show that the crop-climate feedback has contributed to production growth over the past four decades in the central US, but it has also likely boosted production in other major agricultural zones across the world [23]. While substantial global land area is dedicated to agriculture, few places have yields as high as in the central US. Accordingly, crop production-and production-driven LH-is less spatially concentrated in most global regions, meaning that the crop-climate feedback will have a smaller beneficial impact on yields. This implies that while extremely productive US croplands have been somewhat shielded from detrimental global warming over the past 40 years, many low-yielding agricultural regions have been harmed by it, a result highlighted in recent work [25,48]. Thus, the crop-climate feedback may have increased the gap in yields between high-and low-yielding agricultural regions, contributing to global warming's effect on economic inequality [49].
The crop-climate feedback likely acts to increase interannual maize and soy yield variability: in years with poor yields or a smaller harvested area, production-driven LH, cooling, and thus yield benefits are smaller, depressing yields. This effect could play a role in recent increases in drought sensitivity across US croplands [18,50] by effectively making drought years hotter and drier. In contrast, for years with excellent yields or a larger harvested area, LH and cooling are more robust, increasing yields further.
Several limitations of our study may affect our estimation of the strength of the crop-climate feedback. While the counties included in our analysis are intensely cultivated, some land area is covered with natural vegetation, crops other than maize or soy, or urban development. This non-harvested land may affect the estimation of latent heating attributable to maize and soy. This correlation could lead to either an overestimation or an underestimation of production-driven LH, as ET from other vegetation is captured in our coefficient on production (see equation (1) in section 2). However, maize and soy production is heavily concentrated in counties with the highest harvested fractions, meaning that the counties with the most production-driven latent heating also have the least non-crop vegetation. Indeed, we find that the crop-climate feedback's yield effect increases with harvested fraction, implying that the influence of non-crop vegetation should be limited in our results (see SI figure 5(b)). We also do not consider the effect of crop sensitivity to GDD and KDD timing on yields. While such timing is important, it does not change the sign of the yield impact from GDDs and KDDs [5,6]. Furthermore, the cooling effect from production-driven LH affects KDDs much more strongly than GDDs, and recent work has found that KDD timing changes have had minimal influence on maize yield trends [6].
While we show that the crop-climate feedback has played an important role in historical maize and soy production trends, its future effect is uncertain.
Global warming is increasing net radiation, pushing temperatures and, by extension, KDDs and GDDs upwards. Unless continued increases in crop production intensify the crop-climate feedback, this warming trend will further raise KDDs in the southern US and will likely reverse the cooling trend that has occurred over the Upper Midwest, harming yields. The potential for continued maize and soy production growth in the central US is uncertain. In many counties, most or all productive land is already cultivated. Additionally, yield growth may slow due to biophysical limitations, insufficient water, or slower technological change [1,3,33,51]. Together, these challenges suggest that in the coming decades, the ability of the crop-climate feedback to mask some of the global warming trend may fade, exposing crops to more harmful climate conditions.

Data availability statements
The data that support the findings of this study are openly available at the following URL/DOI: https://github.com/ecoffel/2020-ag-climate.
Crop yield and harvested area data are freely available from the USDA NASS QuickStats tool [39]. ERA5-Land climate data are freely available from the European Center for Medium Range Forecasting [34]. GLDAS climate data are freely available from NASA [37]. CPC temperature data are freely available from the CPC [38]. Analysis code available at https://github.com/ecoffel/2020-ag-climate.