Rigden et al. (2020). Combined influence of soil moisture and atmospheric evaporative demand is important for accurately predicting US maize yields

–

7

. Although water is obvi-

ously essential for crop growth and development, the rela-

tionship between moisture availability and yield remains poorly

resolved. Inclusion of precipitation has not, for example, substan-

tially improved statistical yield predictions for maize in the US

Midwest relative to models based on temperature alone

,

. It may

be that precipitation is a poor proxy for plant-available soil mois-

ture because of variable runoff, drainage and evaporation

9

, as well

as the fact that precipitation estimates are generally more uncertain

because rainfall is more heterogeneous than temperature varia-

tions

10

. Nonlinearities in the relationship between plant-available

soil moisture and yield

could also obscure the influence of

moisture availability.

A potentially important clue for the role of moisture in mediat-

ing yield damage comes from prior findings that maize yield begins

to decline by approximately 7% per degree Celsius for air tempera-

ture above approximately 30 °C on a given day

,

8

. If reductions in

maize yield come from physiological heat stress, which occurs

around 35 °C

19

, the leaf surface must be substantially warmer than

the air. Temperature gradients between a canopy and surrounding

air largely depend on available soil moisture, with higher canopy

temperatures coinciding with drier soils

20

–

22

and limited latent cool-

ing. Hence, when atmospheric temperature is below 35 °C, damag-

ing canopy temperatures are likely to be reached only when soil

moisture is low. It can be inferred that although heat stress is often

accompanied by water stress, these variables are not interchangeable

for purposes of predicting yield outcomes.

The inference of moisture-mediated exposure to damage from

temperature is consistent with previous studies. Exposure to high

temperatures has been found to cause amplified yield damage under

conditions of low water availability

13

,

18

,

23

. Conversely, well-irrigated

fields show little sensitivity to high atmospheric temperature

–

.

A further complication is that excess precipitation can cause yield

damage through waterlogging, increased pests and fungal pathogens,

and loss of nitrogen fertilizer

14

,

22

,

. Insurance data over the US sug-

gest that maize yield damage associated with excessive precipitation

is of the same magnitude as that associated with drought conditions

.

Damage from excess soil moisture implies that high tempera-

tures may sometimes benefit yields through increasing evaporative

demand and reducing soil moisture to more optimal conditions.

On the basis of the foregoing considerations, we propose a sta-

tistical model of maize yield based on interactions between the sup-

ply and demand of water. Using this model along with observations

from the Soil Moisture Active Passive (SMAP) satellite

25

, we show

that historical variations in US maize yield are better predicted as a

function of atmospheric demand and root-zone soil moisture than

using standard approaches based on temperature or precipitation

(Table

). Distinguishing the factors controlling moisture balance

is important because there is no assurance that historical covari-

ability between higher atmospheric demand and lower soil mois-

ture will hold with respect to future climate change. Coupled Model

Intercomparison Project Phase 5 (CMIP5) simulations, for example,

give order 100% increases in atmospheric vapour pressure deficit

(VPD) during the Midwestern growing season between 1980–2010

and 2070–2100 under Representative Concentration Pathway 8.5

Combined influence of soil moisture and

atmospheric evaporative demand is important

for accurately predicting US maize yields

A. J. Rigden

1

*, N. D. Mueller

2,3

, N. M. Holbrook

4

, N. Pillai

5

and P. Huybers

1

Understanding the response of agriculture to heat and moisture stress is essential to adapt food systems under climate change.

Although evidence of crop yield loss with extreme temperature is abundant, disentangling the roles of temperature and mois-

ture in determining yield has proved challenging, largely due to limited soil moisture data and the tight coupling between mois-

ture and temperature at the land surface. Here, using well-resolved observations of soil moisture from the recently launched

Soil Moisture Active Passive satellite, we quantify the contribution of imbalances between atmospheric evaporative demand

and soil moisture to maize yield damage in the US Midwest. We show that retrospective yield predictions based on the inter-

actions between atmospheric demand and soil moisture significantly outperform those using temperature and precipitation

singly or together. The importance of accounting for this water balance is highlighted by the fact that climate simulations uni-

formly predict increases in atmospheric demand during the growing season but the trend in root-zone soil moisture varies

between models, with some models indicating that yield damages associated with increased evaporative demand are moder-

ated by increased water supply. A damage estimate conditioned only on simulated changes in atmospheric demand, as opposed

to also accounting for changes in soil moisture, would erroneously indicate approximately twice the damage. This research

demonstrates that more accurate predictions of maize yield can be achieved by using soil moisture data and indicates that accu-

rate estimates of how climate change will influence crop yields require explicitly accounting for variations in water availability.

NAtURe FooD

| VOL 1 | FEBRUARy 2020 | 127–133 |

127

Articles

NATUre Food

(RCP8.5). Corresponding values of root-zone soil moisture, how-

ever, change only by order 10%, with different models variously

simulating wetter or drier conditions. Yield predictions based only

on increased VPD would give twice the damage compared with

predictions that also account for relatively stable levels of root-zone

soil moisture. Improved monitoring and modelling of soil moisture

characteristics therefore appears important for determining how

climate change will influence future yield.

A model combining soil moisture and atmospheric

evaporative demand

A simple calculation is first presented to illustrate how the response

of yield to heat stress depends on soil moisture. We use a standard

statistical approach to estimate yield sensitivities whereby the influ-

ence of atmospheric demand on yield accrues daily

8

, but apply the

model on three distinct sets of data based on seasonally averaged

soil moisture: all of the county-years, the driest 40% of the data, and

the wettest 40% of the data. Heat stress is characterized using daily-

average VPD, as opposed to temperature, because VPD is suggested

to better predict yield loss by more fully capturing water stress

26

When the results from these three regressions are compared (Fig.

),

high VPDs (

>

1.3 kPa) are found to be detrimental only in dry years,

whereas low VPDs (

<

0.9 kPa) are detrimental only in wet years. The

contrasting responses of yield to high and low VPD differs at the

95% confidence interval (CI). Analyses made unconditionally with

respect to soil moisture—that is, using all county-years—would

erroneously indicate that both low and high VPDs are invariably

detrimental.

The fact that the response of yield to VPD depends on soil mois-

ture indicates that a statistical model of yield should account for

both terms. To develop the model, we assume that the response of

yield to environmental variables accumulates daily over the growing

season. However, rather than modelling yield as a function of one

environmental variable, such as temperature

,

8

, we employ a two-

dimensional dependence whereby interactions of VPD and root-

zone soil moisture can be inferred (Fig.

2

). As described in detail

in the Methods, a model is fit wherein yield varies according to the

daily accrued exposure to combinations of soil moisture and VPD.

This formulation is agnostic with respect to how and whether sup-

ply and demand interact, allowing for inferences to be solely derived

from the observations.

Although strong coupling emerges between soil moisture and

VPD at the monthly to seasonal timescales over which land equili-

brates with the atmosphere

27

(Supplementary Fig. 1), these terms

are less correlated at the daily timescales resolved in our analysis

28

(Supplementary Fig. 2), particularly deeper in the column

(Supplementary Fig. 3). Exposure to combinations of soil moisture

and VPD is better resolved through the analysis of daily variability,

as opposed to using monthly or seasonal averages. Process-based

crop models

29

, of course, also resolve both soil moisture and atmo-

spheric demand but the present approach is complementary in

having few enough parameters that the response can be accurately

constrained across the US Midwest.

To examine the viability of the model, we first explore a well-

instrumented, eddy covariance site in Nebraska

30

where measure-

ments of daily root-zone soil moisture, VPD and carbon fluxes

are available. These data are unique among the FLUXNET2015

network for covering more than 3 yr of rainfed maize growth. We

verify that imbalances between soil moisture and VPD lead to

reduced rates of daily gross primary productivity (GPP, Fig.

). In

agreement with the physiology of maize

31

, a balance between root-

zone soil moisture and atmospheric demand gives optimal growing

conditions, corresponding to a high-productivity diagonal. Damage

accrues when water is in short or over supply relative to demand.

We suggest that negative productivity anomalies in the dry-soil and

high-demand quadrant of Fig.

result from stomatal limitation of

photosynthesis. Decreases in leaf water content in maize are asso-

ciated with stomatal closure

32

that prevents xylem tensions from

reaching levels that would cause cavitation

33

but also reduces pro-

ductivity by restricting CO

2

uptake and, hence, photosynthesis

34

.

The negative productivity anomalies in the wet-soil, low-

demand quadrant of Fig.

probably also reflect reduced photo-

synthetic capacity. Some species, such as tomato, respond to sodden

conditions by increasing root hydraulic resistance, but this does not

appear to be the case for maize

35

, where over-supply of water is asso-

ciated with altered photosynthetic biochemistry, including lower

chlorophyll and soluble protein levels, reduced activity of Rubisco

and increased potential for photooxidative damage

36

. Waterlogging

and associated root-zone anoxia can also lead to water and mineral

deficiencies, decreased nitrogen uptake from the soil, and increased

exposure to toxic compounds and disease organisms

. Light

limitations almost certainly play a role in reducing productivity as

well, given the likelihood for wet and cloudy conditions to coincide,

although our representation only implicitly includes light variability

in that atmospheric temperature and humidity control VPD and are

both strong functions of solar radiation.

Table 1 |

Correlation between predicted and observed yield

anomalies

Variables

r

(90% CI)

r

test

VPD, SM

r

0.83 (0.82–0.85)

0.64

T

, SM

r

0.83 (0.82–0.85)

0.59

SM

r

0.80 (0.78–0.81)

0.62

VPD,

P

0.79 (0.78–0.81)

0.58

SM

s

0.79 (0.77–0.81)

0.61

T

,

P

0.78 (0.77–0.80)

0.51

VPD

0.78 (0.76–0.79)

0.59

P

0.77 (0.75–0.79)

0.56

T

0.76 (0.74–0.78)

0.52

The correlation (

r

) between predicted and observed yield anomalies with the associated 90%

CI (from

n

=

10,000 bootstrap realizations) and the average out-of-sample

r

(

r

test

). Explanatory

variables include daily root-zone soil moisture (SM

r

), vapour pressure deficit (VPD), maximum

temperature (

T

), precipitation (

P

) and near-surface soil moisture (SM

s

). See Supplementary Table 1

for an analogous table of root mean square errors.

0.5

Daily yield sensitivity

(bushel per acre per day)

0.7

0.9

1.1

VPD (kPa)

1.3

1.5

1.7

1.9

–1.0

–0.5

0

0.5

1.0

1.5

= 0.1

= 0.3

= 0.5

All years

Dry years

Wet years

Fig. 1 |

Yield sensitivity to VPD.

The yield response to daily anomalies

in vapour pressure deficit (VPD; equation (

)) fitted to all county-years

(black,

n

=

2,552), the driest 40% of county-years (red,

n

=

1,021) and

the wettest 40% of county-years (blue,

n

=

1,021). The driest and wettest

county-years are selected on the basis of the average May–July root-zone
soil moisture. Equal-width bins from 0.4 kPa to 2 kPa are prescribed with
bin centres labelled on the

x

axis. The size of the marker symbols indicates

the fractional distribution of observations in each analysis (see scale
reference on plot), and the vertical bars indicate the 95% CI associated
with each response.

NAtURe FooD

| VOL 1 | FEBRUARy 2020 | 127–133 |

128

Articles

NATUre Food

Retrospective yield predictions across the US Midwest

Satellite observations of soil moisture allow us to extend our analy-

sis to the entirety of the US Midwest. Launched in January 2015,

SMAP provides soil moisture observations at an unprecedented

9-km and approximately 3-d resolution. Daily average soil mois-

ture is observed at a spatial resolution analogous to county-level

maize yields

38

and these outperform other remotely sensed soil

moisture data products when compared with in situ soil mois-

ture measurements

39

,

40

. SMAP observes soil moisture to a depth of

approximately 5 cm, which we use as the surface boundary forcing

of a simple hydrologic model

41

to estimate root-zone soil moisture.

Comparisons of SMAP-derived estimates of root-zone soil moisture

against available in situ observations indicate a high degree of asso-

ciation (Supplementary Fig. 4). A more complex land surface model

could be used to infer root-zone soil moisture but would entail

major assumptions regarding soil characteristics and make connec-

tion with the SMAP data less apparent.

Anomalies in Midwestern maize yields

38

are analysed over the

SMAP interval from 2015 to 2018. To determine yield anomalies

over such a short interval, they are estimated relative to a linear

trend fit to historical data since 1998.

Applying our model to end-of-season yields across 2,552 total

county-years between 2015 and 2018 gives a sensitivity pattern

(Fig.

) similar to that of productivity (Fig.

2

), consistent with a

tight coupling between photosynthesis and yield. Yield sensitivities

indicate that 1 d of sodden conditions induces 0.27 bu per acre of

yield loss, whereas losses associated with high VPDs and dry condi-

tions account for 0.47 bu per acre of yield loss. The yield sensitivity

of sodden conditions is estimated as the average of the eight boxes

under the green diagonal in Fig.

and the sensitivity of dry con-

ditions is estimated as the average of the six boxes over the green

diagonal. Integrating damage between May and July over 2015 to

2018 indicates that approximately half of total yield losses are due to

over-saturated conditions (47%) and half to dry conditions (53%).

These results may partly reflect SMAP’s sampling of unusually wet

conditions during 2015 (Supplementary Fig. 5) but are consistent

with other recent studies covering a longer time interval

.

We examine the predictive skill of our yield model using a leave-

one-out approach (Table

and Supplementary Table 1). If all years

are used to parameterize the model using root-zone soil moisture

and VPD as explanatory variables, predicted yield anomalies are

well correlated with actual yield anomalies (

r

=

0.83 with a 90%

CI: 0.82–0.85; as listed in Table

). Training on each combination of

3 yr of data at the county level and testing predictions in the fourth

leads to an average correlation across the four resulting out-of-

sample tests of 0.64. In contrast, if an analogous model is fitted using

maximum daily temperature and precipitation, the

r

value decreases

from 0.78 when trained on 4 yr of data to 0.51 when tested out-of-

sample, indicating that VPD and soil moisture provide a stronger

basis for predicting the effects of environmental change on yield

outcomes.

Implications of soil moisture for future yield

Yield response to increases in VPD will ultimately depend on the

availability of root-zone soil moisture, but how the balance between

supply and atmospheric demand of moisture will change remains

highly uncertain. Relative to a 1980–2010 climatology, the CMIP5

ensemble predicts increases in vapour pressure in 2070–2100 with

a range of 0.12 to 1.2 kPa in the US Midwest under the RCP8.5,

and changes in root-zone soil moisture that range from

−

0.04 to

0.01 cm

3

cm

−

3

(Fig.

). The increase in VPD reflects an increase

in temperature and saturation vapour pressure, but negligible

changes in actual vapour pressure, due to stable relative humidity

42

.

Increasing VPD will tend to dry near-surface soils, and changes

in VPD are strongly anti-correlated with changes in near-surface

soil moisture across the 24 CMIP5 models (

r

=

−

0.71). Root-zone

soil moisture is, however, less correlated with VPD across models

(

r

=

−

0.52). This decoupling between the deep soil and atmosphere

has been attributed to variations in winter and spring precipitation

persisting into the growing season

43

Increasing VPD by 1.2 kPa for each daily observation reduces

the predicted median yield across all county-years by 64 bu per

acre per season, or approximately a third of the median yield of

176 bu per acre (Fig.

). An increase in root-zone soil moisture of

0.01 cm

3

cm

−

3

, which is the maximum predicted increase across the

CMIP5 models, would decrease losses to 54 bu per acre per season.

Conversely, a decrease in root-zone soil moisture of 0.04 cm

3

cm

−

3

,

which is the maximum predicted decrease across CMIP5 mod-

els, would increase losses to 85 bu per acre per season. If median

yield loss were instead estimated on the basis of a one-dimensional

dependence on VPD (Supplementary Fig. 6), estimated damage

would more than double to 131 bu per acre, signifying the impor-

tance of including the interactions between soil moisture and evap-

orative demand.

Low demand

High demand

Low supply

High supply

Y

′

< 0

Y

′

> 0

Y

′

> 0

Y

′

< 0

Hypothesized yield response

0.25

0.30

0.35

SM

r

(cm

3

cm

–3

)

SM

r

(cm

3

cm

–3

)

GPP

max

′

(

µ

mol m

–2

s

–1

)

0.5

1.0

VPD (kPa)

1.5

2.0

2.5

36

70

41

32

40

9

60

41

6

–2

0

2

Observed productivity response

a

b

Fig. 2 |

Moisture supply and demand drive yield.

a

, The observed daily response of maize productivity to daily vapour pressure deficit (VPD) and root-

zone soil moisture (SM

r

) measured at 25 cm.

b

, Hypothesized daily yield sensitivity (

Y

′

) to daily VPD and SM

r

conditions. In

a

, the colour indicates the

average anomaly in daily GPP

max

(GPP

max

′

) and the number indicates the number of days in each bin.

NAtURe FooD

| VOL 1 | FEBRUARy 2020 | 127–133 |

129

Articles

NATUre Food

Conclusions

Given its importance in determining yield, both with respect to

drought

44

and sodden

conditions, continued improvement in the

representation of root-zone soil moisture in general circulation

models appears important for gauging the implication of climate

change for food production. This study also implies that adapta-

tion strategies that optimize the root-zone water balance, including

plants with deeper roots

45

or crops that better conserve root-zone

soil moisture under conditions of high atmospheric demand

46

, are

especially important for ensuring food security in future climates.

Methods

Data sources. To characterize soil moisture in each county, we utilize the

SMAP Level-3 9-km radiometer soil moisture product

25

, which has an exact

repeat frequency of 8 d and a nominal return frequency of 3 d. To minimize

retrieval errors, soil moisture estimates are omitted if the water fraction exceeds

5% of the pixel or if the vegetation density exceeds 5 kg m

−

2

. It is worth noting

that ideally soil moisture values estimated over the 9-km footprint would be

homogeneous and representative exclusively of croplands. In reality, noise could

be introduced by vegetation or other surfaces surrounding croplands that have

different dry-down dynamics. Daily average VPD, which is approximated as the

average of the maximum and minimum daily VPD, maximum temperature (

T

),

and precipitation accumulation (

P

) are acquired from the Parameter-elevation

Relationships on Independent Slopes Model (PRISM) dataset

47

. The PRISM data

are originally derived from weather station measurements. The data were then

spatially interpolated to a 4-km grid using an algorithm that accounts for elevation

and physiographic factors (for example, coastal proximity and topographic facet

orientation)

47

. VPD and temperature measurements used in the PRISM dataset are

from instrument height, which is typically 1–2 m above the land surface. Although

instrument height measurements are commonly used to model vegetation

and yield, we would use leaf-level humidity and temperature values if such

measurements were available, as they are more physiologically relevant. If we use

data from weather stations instead of PRISM, we find similar relationships between

VPD, soil moisture and yield (Supplementary Fig. 7).

County-level maize yield data are from the US Department of Agriculture/

National Agricultural Statistics Service

38

. If counties reported non-irrigated maize

yields in addition to total yields, the non-irrigated yields are used. Non-irrigated

yields are preferentially chosen because we are interested in explaining yield

variability in rainfed areas, and the non-irrigated soil moisture conditions are

probably more representative of the 9-km SMAP soil moisture observation.

Data from PRISM and SMAP are aggregated to the county scale using county

boundary maps obtained from the 2016 US Census Bureau. Specifically, we average

over all grid box centres within the county boundary. Our results do not change

if we averaged only grid boxes categorized as dominantly ‘croplands’ or ‘cropland

mosaics’ (Supplementary Fig. 7).

Projections of monthly relative humidity (variable name ‘hurs’ in CMIP5),

temperature (‘tas’) and multi-level soil moisture (‘mrlsl’) are obtained from CMIP5

RCP8.5. Models were chosen on the basis of the availability of monthly multi-layer

soil moisture and surface-level humidity and temperature. Monthly average VPD is

estimated from near-surface relative humidity and temperature. To characterize the

plant-available water, root-zone soil moisture is estimated between 0.2 and 1.0 m by

linearly interpolating soil moisture between reported levels. Surface soil moisture is

defined between 0 and 0.1 m. Soil depth and the number of levels in the soil profile

are model dependent (Supplementary Table 2). Averages are calculated over two

different intervals, 1980–2010 and 2070–2100 (Fig.

).

Estimating root-zone soil moisture. The soil moisture estimates from SMAP

represent near-surface soil moisture to a depth of approximately 5 cm (ref.

25

). To

more accurately characterize plant-available moisture, which has been shown to

extend below 5 cm for maize

48

, we estimate a root-zone soil moisture at time

t

n

(SM

r,

n

) using an auto-regressive exponential framework where innovations are

from the surface

41

SM

r

;

n

¼

SM

r

;

n

�

1

þ

K

n

SM

s

;

n

�

SM

r

;

n

�

1

ð

1

Þ

In equation (

), SM

r,

n

−

1

is the predicted root-zone soil moisture estimate at

t

n

−

1

,

SM

s,

n

is the surface soil moisture retrieved from SMAP at time

t

n

and

K

n

is the gain

at time

t

n

estimated as

K

n

¼

K

n

�

1

K

n

�

1

þ

e

tn

�

tn

�

1

N

ð

2

Þ

In equation (

2

),

N

represents the timescale of soil moisture variation in units

of days. Although the optimal value of

N

varies spatially, it has been demonstrated

that using an average

N

does not significantly decrease model accuracy

41

49

. Thus,

we assume a spatially uniform

N

=

8 d, which is based on calibrations from the Soil

Moisture and Ocean Salinity satellite over Nebraska and Oklahoma using near-

surface (5–10 cm) and root-zone (25–60 cm) soil moisture measurements

49

. Results

and conclusions are not sensitive to the choice of parameter

N

(Supplementary

Fig. 7). We initialize the exponential filter with SM

r,1

=

SM

s,1

and

K

1

=

1. The

recursive exponential filter smooths the near-surface soil moisture observations,

as demonstrated for one SMAP 9-km grid box in Supplementary Fig. 8. We also

evaluate this methodology using soil moisture observations at four depths from

24 US Climate Reference Network sites

50

51

, finding a high degree of association

(Supplementary Fig. 4). This model of root-zone soil moisture is simplified

compared to those in land surface and climate models, which represent

hydrologic processes in more detail and, thus, have longer soil moisture memory.

To 2.5 kPa

2.2

2.0

1.8

1.6

1.4

1.2

1.0

0.8

0.6

0.4

0.2

0.15

0.20

0.25

0.30

SM

r

(cm

3

cm

–3

)

VPD (kP

a)

0.35

0.40

Daily yield sensitivity (b

u

d

–1

)

0.2

–0.2

–0.4

–0.6

–0.8

–1.0

0

2015–2018

observations

0

–50

–100

–150

0

0.2

f

(SM

r

, VPD)

f

(VPD)

0.4

0.6

Increase in VPD (kPa)

0.8

1.0

1.2

Yi

eld response

(b

ushel per acre per season)

a

b

Fig. 3 |

Maize yield response to moisture supply and demand in present

and future climates.

a

, Using observations, we estimate the yield response

to vapour pressure deficit (VPD) and root-zone soil moisture (SM

r

). Each

square represents inferred daily yield sensitivities (that is,

β

h

coefficients)

that are estimated via observational data and the constraints in equation (

).

The size of the square represents the relative number of days within each
bin, such that the largest box in the centre represents 12.5% of the data
and the smallest box in the upper right corner represents less than 1% of
the data (see Supplementary Fig. 13 for the distributions of days in each
bin for each year). A black dot within the square indicates that the fit
coefficient is significant (

P

<

0.05). A polynomial surface fit to these daily

yield sensitivities is plotted behind the squares. Residuals of the polynomial
fit are shown in Supplementary Fig. 10. Each orange arrow represents the
change in SM

r

and VPD from 1980–2010 to 2070–2100 from one CMIP5

model (see Supplementary Table 2 for a list of models). For reference,
the average of summertime observations (2015–2018) is plotted with
an open circle and labelled.

b

, We predict the yield response to elevated

VPD by increasing the observed VPD and estimating yields using a fitted
polynomial function for VPD (magenta) and a fitted polynomial surface for
VPD and SM

r

(teal). The shaded CIs represent the 10th to 90th percentile

from a distribution of possible median yields estimated by bootstrapping
data according to county (

n

=

10,000 realizations). The polynomial function

for VPD is shown in Supplementary Fig. 6, and the polynomial surface for
VPD and SM

r

is shown above.

NAtURe FooD

| VOL 1 | FEBRUARy 2020 | 127–133 |

130

Articles

NATUre Food

A more complex land surface model could be used to infer root-zone soil moisture

in this study, but would entail major assumptions regarding soil characteristics, and

could make connection with the SMAP data less apparent.

Estimation of yield anomalies. Maize yield anomalies are estimated for 2015 to

2018 in each county as follows: fitting a linear regression to the yield in each county

from 1998 to 2014 to determine the long-term trend external to climate variability;

and subtracting the predicted yields from the observed yields, respectively. In the

first step, we include only counties with over 14 yr of data excluding 2012, which

had extremely low yields due to severe drought. We obtain similar results if yield

anomalies are estimated relative to the 2015–2018 average (Supplementary Fig. 7).

First, we resolve the univariate yield response to individual climate variables,

including daily near-surface soil moisture, root-zone soil moisture, precipitation,

VPD and maximum temperature. To facilitate comparisons between the climate

variables, we analyse only days with available soil moisture observations (as there

are missing data only in soil moisture), and counties that have a minimum of 21 d

of available soil moisture observations from May to July. Note that we do not fill

missing soil moisture data because preserving the covariance between climate

variables is critical to accurately resolve how they interact to affect yield. In total,

there are 2,552 county-years that meet these criteria, aggregating across all counties

east of 100° W, west of 80° W and north of 35° N and years 2015–2018.

To model the yield response for each climate variable, we assume the effects are

cumulative over time during the months of May, June and July, and thus additively

substitutable over time

. We model yield as a function of exposure time to climate

conditions, rather than fitting a parametric yield model to the climatic variable, to

avoid prescribing a functional relationship. For each climate variable, we fit a step

function to the county-level yield anomalies

Y

0

c

¼

X

n

h

¼

1

ρ

h

;

c

β

h

þ

f

c

þ

ϵ

c

ð

3

Þ

In equation (

),

ρ

h

represents the exposure time in units of days to each of

n

different environmental states measured for each county, c;

β

h

represents sensitivity

to different exposures in units of bushels per acre per day and is held fixed across

counties and years;

f

c

represents county-level fixed effects and

ϵ

c

I

represents the

error. To account for the fact that soil moisture observations are not necessarily

available every day and counties may have more or less missing data, we scale

ρ

h

assuming that the available data are representative of the whole growing season.

Specifically, the value of

ρ

h

is estimated for each county-year as

N

b

/

N

o

×

N

g

, where

N

b

is the number of days that fall in the bin,

N

o

is the total number of days in

which soil moisture is observed, and

N

g

is the total number of days in the growing

season (

N

g

=

92 d). The

β

h

and

f

c

values are fitted using a linear mixed-effects model

without a global intercept. We estimate 95% CIs for

β

h

coefficients (Fig.

), and

P

values for a hypothesis test of whether each parameter is equal to zero (Fig.

). To

estimate CIs on the correlation coefficients (Table

), we utilize non-parametric

bootstrapping techniques wherein data from individual counties are resampled in

10,000 iterations. These CIs have a slight asymmetry due to the fixed effects

52

To illustrate the dependence of heat stress on water availability (Fig.

), we fit

equation (

) to VPD using three sets of data: all of the county-years, the driest 40%

of county-years based on seasonally averaged soil moisture, and the wettest 40% of

county-years. For visual clarity, we fit the yield response using eight equal-width

bins ranging from 0.4 to 2 kPa, although this response function looks similar to the

one fitted with five bins (Supplementary Fig. 6).

To compare the yield response to different climate variables, May–July daily

observations are binned into five equal-width bins. We use five bins to allow for

comparisons between models with interactions, which increase the number of

bins by its square (from 5 to 25), as described in the next paragraph. When fitting

equation (

), bin edges are defined on the basis of the 1% and 99% quantiles, which

we estimate using all days in May to July across all 2,552 county-years for each

explanatory variable. Bins range from 0.1 to 0.5 cm

3

cm

−

3

for surface soil moisture,

0.14 to 0.44 cm

3

cm

−

3

for root-zone soil moisture, 0.14 to 2.1 kPa for VPD, 10 to

36 °C for temperature and 0 to 40 mm for precipitation. Data outside bin edges

are counted in the edge bin. Five bins and their respective ranges are selected to

ensure sufficient data resides within each to allow for sensitivity estimates, but the

conclusions of the study are insensitive to the specification of the bin edges and

number of bins, provided that the bins adequately span the observational range of

data (Supplementary Fig. 7).

To incorporate moisture supply and demand interactions, data are binned

conditioned on both supply (root-zone soil moisture or precipitation) and demand

(VPD or temperature), resulting in 25 bins in total (for example, 5 soil moisture

bins (

h

) by 5 VPD bins (

i

))

Y

0

c

¼

X

n

h

¼

1

X

m

i

¼

1

ρ

h

;

i

;

c

β

h

;

i

þ

f

c

þ

ϵ

c

ð

4

Þ

In this two-dimensional framework,

ρ

h,i

represents an exposure matrix and

β

h

represents a sensitivity matrix whose product is summed across rows and columns.

The range of environmental data is maintained as previously stated. As more

soil moisture data become available from SMAP, additional covariates should be

explored, such as interactions with solar radiation

13

. Residuals from equation (

)

fitted to root-zone soil moisture and VPD are shown in Supplementary Fig. 9.

Underlying our statistical framework is the assumption that moisture

availability influences yield variability similarly from May to July. Field-level

studies have shown that certain growth stages (for example, tasselling and ear

formation) are more sensitive to water stress

23

, although the general conclusions

of this study do not change on the basis of the summertime months selected

(Supplementary Table 3).

Extrapolating yield model to future climate conditions. To extrapolate the yield

response beyond environmental conditions observed in 2015–2018 (Fig.

), we

fit polynomial functions to the derived daily yield response (

β

h

) from equations

(

) and (

). Specifically, we fit a quadratic curve to the five VPD coefficients from

equation (

), as

β

h

VPD

ð

Þ ¼

W

h

q

0

þ

q

1

´

VPD

h

þ

q

2

´

VPD

2

h

ð

5

Þ

In equation (

5

), we incorporate unequal weights,

W

h

, that reflect the fraction

of total number of days in each bin summed across the 2,552 county-years,

N

h

,

to the total number of observations summed across all 2,552 county-years,

Σ

N

h

,

as

W

h

=

N

h

/

Σ

N

h

. This quadratic fit is shown in Supplementary Fig. 6. To account

for the interactions between VPD and soil moisture, we fit a two-dimensional

polynomial surface of order two for VPD and root-zone soil moisture to the 25

coefficients in equation (

)

β

h

VPD

;

SM

ð

Þ ¼

W

h

p

0

þ

p

1

´

SM

h

þ

p

2

´

VPD

h

þ

p

3

´

SM

2

h

þ

p

4

´

VPD

h

´

SM

h

þ

p

5

´

VPD

2

h

ð

6

Þ

Weights are similarly based on the number of observations in each of the 25

bins. The weights are important so that coefficients based on a small number of

observations are not over-represented in the fit, and explain why, for example, the

coefficient in the low soil moisture and low VPD combination differs from the

polynomial fit. The polynomial fit with interactions is shown in the background

in Fig.

and the residuals are shown in Supplementary Fig. 10. These polynomial

fits allow us to estimate yield under elevated VPD conditions (Fig.

) that extend

beyond current observational ranges. To estimate uncertainty in the polynomial

fits, we fit polynomials to 10,000 resampled versions of the

β

h

coefficients in

equations (

) and (

) obtained from bootstrapping yield through resampling

according to county. Yield response is predicted for an increase in VPD for each

county and the median yield is recorded as a function of VPD, with the shaded

CIs in Fig.

showing the 10th to 90th percentile of the distribution of possible

median yields given uncertainties in regression coefficients.

When extrapolating the yield response, we focus on VPD and root-zone soil

moisture because these are physiologically relevant variables

26

that are the most

skillful out-of-sample predictors. When applied in the context of climate change,

we assume that the joint distribution between VPD and soil moisture is unchanged

except for a shift in the mean and that sensitivity to these variables remains fixed.

There is the potential, however, for altered sensitivities from increased carbon

dioxide concentrations and associated improved water relations

53

that should be

explored in future work.

Out-of-sample testing. To test the predictive ability of the yield model in equations

(

) and (

), we train the model on 3 yr of data and test it on the remaining year (for

example, train using 2015–2017 data, and test on 2018) and average the statistics

across all four sets of out-of-sample tests (Table

and Supplementary Table 1). The

fits for temperature or VPD and root-zone soil moisture or precipitation are shown

in Supplementary Figs. 11 and 12.

Eddy covariance analysis. We explore the relationship between daily VPD, soil

moisture and GPP using data from the Mead rainfed eddy covariance site in

Nebraska

30

(site abbreviation: US-Ne3). Isolating the years with maize (2001,

2003, 2005, 2007, 2009 and 2011), we estimate anomalies in daily maximum GPP

(GPP

max

) during the growing season. We define the growing season by smoothing

GPP

max

with a two-week window and isolating days when the smoothed seasonal

cycle of GPP

max

exceeds 40 umol m

−

2

s

−

1

. To quantify the influence of soil moisture,

measured at 25 cm, and VPD on productivity (Fig.

https://nsidc.org/data/smap

), we average GPP

max

anomalies on the basis of the daily daytime VPD and soil moisture status, with

daytime hours identified by when incoming solar radiation exceeds 50 W m

−

2

. In

the conditional averaging scheme, we use 3

×

3 bins due to the smaller number of

observations at the eddy covariance site relative to county-year yield data.

Reporting Summary. Further information on research design is available in the

Nature Research Reporting Summary linked to this article.

Data availability

SMAP science data are available at the NASA National Snow and Ice Data Center

Distributed Active Archive Center (

). Meteorological

data are available at the PRISM Climate Group at Oregon State University

NAtURe FooD

| VOL 1 | FEBRUARy 2020 | 127–133 |

http://prism.oregonstate.edu

131

Articles

NATUre Food

(

). Maize yield data are available at the US Department

of Agriculture/National Agriculture Statistics service (

https://www.nass.usda.gov

).

Observations from the US-Ne3 eddy covariance site are available at the AmeriFlux

(soil moisture;

http://ameriflux.lbl.gov

) and FLUXNET2015 (VPD, GPP and solar

radiation;

http://fluxnet.fluxdata.org/

) data repositories.

Received: 23 July 2019; Accepted: 6 January 2020;

Published online: 18 February 2020

References

1. Butler, E. E. & Huybers, P. Variations in the sensitivity of US maize yield to

extreme temperatures by region and growth phase.

Environ. Res. Lett.

10,

034009 (2015).

2. Handmer, J. et al. in

Managing the Risks of Extreme Events and Disasters to

Advance Climate Change Adaptation

(eds Field, C. B. et al.) 231–290 (IPCC,

Cambridge Univ. Press, 2012).

3. 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).

4. 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).

5. Porter, J. R. et al. in

Climate Change 2014: Impacts, Adaptation, and

Vulnerability.

(eds Field, C. B. et al.) 485–533 (IPCC, Cambridge Univ.

Press, 2014).

6. Schlenker, W. & Roberts, M. J. Nonlinear temperature effects indicate severe

damages to U.S. crop yields under climate change.

Proc. Natl Acad. Sci. USA

106, 15594–15598 (2009).

7. Zhao, C. et al. Temperature increase reduces global yields of major crops in

four independent estimates.

Proc. Natl Acad. Sci. USA

114, 9326–9331 (2017).

8. Butler, E. E. & Huybers, P. Adaptation of US maize to temperature variations.

Nat. Clim. Change

3, 68–72 (2013).

9. Basso, B. & Ritchie, J. T. Temperature and drought effects on maize yield.

Nat. Clim. Change

4, 233 (2014).

10. Fezzi, C. & Bateman, I. The impact of climate change on agriculture:

nonlinear effects and aggregation bias in Ricardian models of farmland

values.

J. Assoc. Environ. Resour. Econ.

2, 57–92 (2015).

11. Anderson, C. J., Babcock, B. A., Peng, Y., Gassman, P. W. & Campbell, T. D.

Placing bounds on extreme temperature response of maize.

Environ. Res. Lett.

10, 124001 (2015).

12. Carter, E. K., Melkonian, J., Riha, S. J. & Shaw, S. B. Separating heat stress

from moisture stress: analyzing yield response to high temperature in

irrigated maize.

Environ. Res. Lett.

11, 094012 (2016).

13. Carter, E. K., Melkonian, J., Steinschneider, S. & Riha, S. J. Rainfed maize

yield response to management and climate covariability at large spatial scales.

Agric. For. Meteorol.

256–257, 242–252 (2018).

14. Rosenzweig, C., Tubiello, F. N., Goldberg, R., Mills, E. & Bloomfield, J.

Increased crop damage in the US from excess precipitation under climate

change.

Global Environ. Change

12, 197–202 (2002).

15. Schauberger, B. et al. Consistent negative response of US crops to high

temperatures in observations and crop models.

Nat. Commun.

8,

1–9 (2017).

16. Shaw, S. B., Mehta, D. & Riha, S. J. Using simple data experiments to explore

the influence of non-temperature controls on maize yields in the mid-West

and Great Plains.

Clim. Change

122, 747–755 (2014).

17. Troy, T. J., Kipgen, C. & Pal, I. The impact of climate extremes and irrigation

on US crop yields.

Environ. Res. Lett.

10, 054013 (2015).

18. Urban, D. W., Sheffield, J. & Lobell, D. B. The impacts of future climate and

carbon dioxide changes on the average and variability of US maize yields

under two emission scenarios.

Environ. Res. Lett.

10, 045003 (2015).

19. Siebert, S., Webber, H., Zhao, G. & Ewert, F. Heat stress is overestimated in

climate impact studies for irrigated agriculture.

Environ. Res. Lett.

12,

054023 (2017).

20. Siebert, S., Ewert, F., Eyshi Rezaei, E., Kage, H. & Graß, R. Impact of heat

stress on crop yield—on the importance of considering canopy temperature.

Environ. Res. Lett.

9, 044012 (2014).

21. Webber, H. et al. Simulating canopy temperature for modelling heat stress in

cereals.

Environ. Model. Softw.

77, 143–155 (2016).

22. Kaur, G., Zurweller, B. A., Nelson, K. A., Motavalli, P. P. & Dudenhoeffer, C. J.

Soil waterlogging and nitrogen fertilizer management effects on corn and

soybean yields.

Agron. J.

109, 97–106 (2017).

23. Ortiz-Bobea, A., Wang, H., Carrillo, C. M. & Ault, T. R. Unpacking the

climatic drivers of US agricultural yields.

Environ. Res. Lett.

14,

064003–064013 (2019).

24. Li, Y., Guan, K., Schnitkey, G. D., DeLucia, E. & Peng, B. Excessive rainfall

leads to maize yield loss of a comparable magnitude to extreme drought in

the United States.

Global Change Biol.

5, 143–113 (2019).

25. Entekhabi, D. et al. The Soil Moisture Active Passive (SMAP) mission.

Proc. IEEE

98, 704–716 (2010).

26. Hsiao, J., Swann, A. L. S. & Kim, S.-H. Maize yield under a changing

climate: the hidden role of vapor pressure deficit.

Agric. For. Meteorol.

279,

107692 (2019).

27. Seneviratne, S. I. et al. Investigating soil moisture–climate interactions in a

changing climate: a review.

Earth Sci. Rev.

99, 125–161 (2010).

28. Novick, K. A. et al. The increasing importance of atmospheric demand for

ecosystem water and carbon fluxes.

Nat. Clim. Change

6, 1023–1027 (2016).

29. Tao, F., Yokozawa, M. & Zhang, Z. Modelling the impacts of weather and

climate variability on crop productivity over a large area: a new process-based

model development, optimization, and uncertainties analysis.

Agric. For.

Meteorol.

149, 831–850 (2009).

30. Suyker, A.

AmeriFlux US-Ne3 Mead – Rainfed Maize-Soybean Rotation Site

from 2001–Present

(AmeriFlux, 2016);

http://ameriflux.lbl.gov/sites/siteinfo/

US-Ne3

31. Kramer, P. J. & Boyer, J. S.

Water Relations of Plants and Soils

(Elsevier, 1995).

32. Bennett, J. M., Sinclair, T. R., Muchow, R. C. & Costello, S. R. Dependence

of stomatal conductance on leaf water potential, turgor potential, and

relative water content in field-grown soybean and maize.

Crop Sci.

27,

984–990 (1987).

33. Cochard, H. Xylem embolism and drought-induced stomatal closure in

maize.

Planta

215, 466–471 (2002).

34. Farquhar, G. D. & Sharkey, T. D. Stomatal conductance and photosynthesis.

Annu. Rev. Plant Physiol.

33, 317–345 (1982).

35. Wenkert, W., Fausey, N. R. & Watters, H. D. Flooding responses in

Zea mays

L.

Plant Soil

62, 351–366 (1981).

36. Yordanova, R. Y. & Popova, L. P. Flooding-induced changes in photosynthesis

and oxidative status in maize plants.

Acta Physiol. Plant.

29, 535–541 (2007).

37. Voesenek, L. A. C. J. & Bailey-Serres, J. Flood adaptive traits and processes:

an overview.

New Phytol.

206, 57–73 (2015).

38.

National Agricultural Statistics Service

(US Department of Agriculture, 2019);

https://www.nass.usda.gov

39. Cui, C. et al. Soil moisture mapping from satellites: an intercomparison of

SMAP, SMOS, FY3B, AMSR2, and ESA CCI over two dense network regions

at different spatial scales.

Remote Sens.

10, 33–19 (2018).

40. Ma, H. et al. Satellite surface soil moisture from SMAP, SMOS, AMSR2 and

ESA CCI: a comprehensive assessment using global ground-based

observations.

Remote Sens. Environ.

231, 111215 (2019).

41. Albergel, C. et al. From near-surface to root-zone soil moisture using an

exponential filter: an assessment of the method based on in-situ observations

and model simulations.

Hydrol. Earth Syst. Sci.

12, 1323–1337 (2008).

42. Ficklin, D. L. & Novick, K. A. Historic and projected changes in vapor

pressure deficit suggest a continental-scale drying of the United States

atmosphere.

J. Geophys. Res. Atmos.

122, 2061–2079 (2017).

43. Berg, A., Sheffield, J. & Milly, P. C. D. Divergent surface and total soil

moisture projections under global warming.

Geophys. Res. Lett.

44,

236–244 (2017).

44. Ukkola, A. M. et al. Evaluating CMIP5 model agreement for multiple drought

metrics.

J. Hydrometeorol.

19, 969–988 (2018).

45. Kell, D. B. Breeding crop plants with deep roots: their role in sustainable

carbon, nutrient and water sequestration.

Ann. Bot.

108, 407–418 (2011).

46. Messina, C. D. et al. Limited-transpiration trait may increase maize drought

tolerance in the US corn belt.

Agron. J.

107, 1978 (2015).

47. Daly, C. et al. Physiographically sensitive mapping of climatological

temperature and precipitation across the conterminous United States.

Int. J. Climatol.

28, 2031–2064 (2008).

48. Nicoullaud, B., King, D. & Tardieu, F. Vertical distribution of maize roots in

relation to permanent soil characteristics.

Plant Soil

159, 245–254 (1994).

49. Ford, T. W., Harris, E. & Quiring, S. M. Estimating root zone soil moisture

using near-surface observations from SMOS.

Hydrol. Earth Syst. Sci.

18,

139–154 (2014).

50. Bell, J. E. et al. U.S. climate reference network soil moisture and temperature

observations.

J. Hydrometeorol.

14, 977–988 (2013).

51. Diamond, H. J. et al. U.S. climate reference network after one decade

of operations: status and assessment.

Bull. Am. Meteorol. Soc.

94,

485–498 (2013).

52. Slaets, J. I. F., Piepho, H.-P., Schmitter, P., Hilger, T. & Cadisch, G.

Quantifying uncertainty on sediment loads using bootstrap confidence

intervals.

Hydrol. Earth Syst. Sci.

21, 571–588 (2017).

53. Leakey, A. D. B. et al. Elevated CO

2

effects on plant carbon, nitrogen,

and water relations: six important lessons from FACE.

J. Exp. Bot.

60,

2859–2876 (2009).

Acknowledgements

We acknowledge the World Climate Research Programme’s Working Group on Coupled

Modelling, which is responsible for CMIP, and we thank the climate modelling groups

(listed in Supplementary Table 2) for producing and making available their model

output. For CMIP, the US Department of Energy’s Program for Climate Model Diagnosis

and Intercomparison provides coordinating support and led development of software

infrastructure in partnership with the Global Organization for Earth System Science

NAtURe FooD

| VOL 1 | FEBRUARy 2020 | 127–133 |

132

Articles

NATUre Food

Portals. We also acknowledge the US Department of Energy’s Office of Science for

funding the AmeriFlux data resources, as well as the FLUXNET community, whose

data processing and harmonization was carried out by the European Fluxes Database

Cluster, the AmeriFlux Management Project and the Fluxdata project of FLUXNET,

with the support of the CDIAC and ICOS Ecosystem Thematic Center, and the OzFlux,

ChinaFlux and AsiaFlux offices. N.P. acknowledges partial financial support from the

Office of Naval Research. A.J.R. acknowledges financial support from the Rockefeller

Foundation Planetary Health Fellows programme at Harvard University.

Author contributions

A.J.R. and P.H. conceived and drafted the manuscript, A.J.R. performed the analysis, and

all authors contributed to the writing and interpretation.

Competing interests

The authors declare no competing interests.

Additional information

Supplementary information is available for this paper at

https://doi.org/10.1038/

s43016-020-0028-7

Correspondence and requests for materials should be addressed to A.J.R.
Reprints and permissions information is available at

www.nature.com/reprints

.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in

NAtURe FooD

| VOL 1 | FEBRUARy 2020 | 127–133 |