The 2011–16 California drought illustrates that drought-prone areas do not always experience relief once a favorable phase of El Niño–Southern Oscillation (ENSO) returns. In the twenty-first century, such an expectation is unrealistic in regions where global warming induces an increase in terrestrial aridity larger than the changes in aridity driven by ENSO variability. This premise is also flawed in areas where precipitation supply cannot offset the global warming–induced increase in evaporative demand. Here, atmosphere-only experiments are analyzed to identify land regions where aridity is currently sensitive to ENSO and where projected future changes in mean aridity exceed the range caused by ENSO variability. Insights into the drivers of these changes in aridity are obtained using simulations with the incremental addition of three different factors to the current climate: ocean warming, vegetation response to elevated CO2 levels, and intensified CO2 radiative forcing. The effect of ocean warming overwhelms the range of ENSO-driven temperature variability worldwide, increasing potential evapotranspiration (PET) in most ENSO-sensitive regions. Additionally, about 39% of the regions currently sensitive to ENSO will likely receive less precipitation in the future, independent of the ENSO phase. Consequently aridity increases in 67%–72% of the ENSO-sensitive area. When both radiative and physiological effects are considered, the area affected by arid conditions rises to 75%–79% when using PET-derived measures of aridity, but declines to 41% when an aridity indicator for total soil moisture is employed. This reduction mainly occurs because plant stomatal resistance increases under enhanced CO2 concentrations, resulting in improved plant water-use efficiency, and hence reduced evapotranspiration and soil desiccation. Imposing CO2-invariant stomatal resistance may overestimate future drying in PET-derived indices.
Observed climate change arises from the combined effects of external forcings—both anthropogenic and natural—and chaotic internal climate variability. To assess the vulnerability of different managed systems (e.g., energy, agriculture, and water supply) to future climate change, and to inform sound adaptation strategies, it is essential to examine changes in the mean state and background variability jointly, rather than individually. This is particularly important for evaluating perturbations to the hydrological cycle, where inherent variability and expected mean state changes are spatially and seasonally complex.
a. Changes in precipitation
The globally averaged tropospheric radiative cooling rate is balanced by surface sensible heat fluxes and by the latent heating released during precipitation formation. This balance provides a fundamental constraint on the changes in global mean precipitation P (Allen and Ingram 2002; Pendergrass and Hartmann 2014; Thorpe and Andrews 2014; Trenberth 2011). The P response to future climate change can be divided into two components: a “fast” response associated with the atmospheric component of CO2 radiative forcing and a “slow” response to the global surface warming independent of the forcing mechanism (Andrews et al. 2010, 2011; Myhre et al. 2017). The fast response arises because increased CO2 concentrations lead to an immediate decrease in atmospheric longwave radiative cooling (Pendergrass and Hartmann 2014) and an increase in tropospheric heating, which is largely balanced by a rapid decrease in latent heating and hence global P (Bony et al. 2013; Colman 2015; Kamae et al. 2015; Richardson et al. 2016; Samset et al. 2016; Sherwood 2015; Tian et al. 2017). The slow response occurs because global warming enhances longwave radiative cooling by the atmosphere, which is largely balanced by additional latent heating and precipitation, both of which increase by roughly 2%–3% K−1 of global warming (Allen and Ingram 2002; Andrews et al. 2010; Bala et al. 2010; DeAngelis et al. 2016).
On a regional scale, the fast P response to CO2 forcing can be complex (e.g., Richardson et al. 2016). The warming-induced changes in P are also complex, driven by multiple physical mechanisms. These include the “wet get wetter” mechanism associated with an increase in tropospheric water vapor (Held and Soden 2006) and the “warmer get wetter” mechanism linking the patterns of tropical sea surface temperature (SST) changes with rainfall changes (Xie et al. 2010). A poleward displacement of current zonal-mean wet and dry patterns is also expected as a result of shifts in atmospheric circulation (Marvel and Bonfils 2013; Polson et al. 2013; Seidel et al. 2008), while terrestrial moisture recycling introduces meridional structure in the patterns of P change over land (Greve et al. 2014). This meridional structure is largely absent in the zonal changes over oceans (Durack et al. 2012). The resulting forced pattern of P change differs markedly from patterns arising from changes in P variability alone (Power et al. 2013; Zhou et al. 2014).
Given the spatial complexity of changes in the mean state and variability of P, and in view of potential interactions between forced and unforced changes (Maher et al. 2014), it is important to quantify the relative importance of different contributory factors to future P changes. Such an analysis was performed by Bonfils et al. (2015), who decomposed the twenty-first-century P response into the sum of three components: 1) the long-term change in mean P, 2) the historical twentieth-century P response to ENSO variability, and 3) an intensification of ENSO-driven P responses in the twenty-first century. Understanding where and when the projected changes in P will rise above the historical range of ENSO-driven variability can be of great benefit for assessing the potential vulnerability of agriculture, water resources, and certain energy sectors (such as hydropower) to climate change.
b. Changes in evapotranspiration
In drought research, understanding changes in the mean and variability of P only informs the “supply side.” It is equally important to assess the behavior of actual evapotranspiration1 (ET). On time scales longer than a few months, global mean ET increases at 2%–3% K−1 of global warming to be in balance with changes in P. Over oceans, where there is an infinite moisture source, the fast response to rising CO2 concentrations is an increase in low-level atmospheric stability and near-surface humidity, which lead to a suppression of ET (Cao et al. 2012; Kamae and Watanabe 2012). In contrast, the warming of the land surface decreases the vertical stability of the lower atmosphere, increases the export of moisture flux from the boundary layer, and enhances convection. These changes tend to increase ET and P (Cao et al. 2012), particularly in tropical regions (Richardson et al. 2016).
Over land, however, an important competing mechanism exists, involving transpiration of water by vegetation (Schlesinger and Jasechko 2014). Plant stomata regulate transpiration; the degree of regulation depends on the availability of energy (in the form of radiation), water (in the soil and as atmospheric moisture), and CO2 (Farquhar and Sharkey 1982).
c. Impact of vegetation on changes in P and ET
Under elevated atmospheric CO2 levels, plants optimize photosynthesis and water loss per unit leaf area by reducing the opening of stomata in their leaves. Elevated CO2 induces an increase in water-use efficiency (WUE), which is the amount of CO2 absorbed by photosynthesis per amount of water lost during transpiration (Goudriaan and Unsworth 1990).2 The reduced transpiration and surface latent heat flux to the atmosphere lead to a significant drying and warming of the boundary layer over land, and a decrease in P (Andrews et al. 2011; Cao et al. 2012; Dong et al. 2009; Doutriaux-Boucher et al. 2009). These physiological effects can also modulate the relative strengths of latent and sensible heat flux responses to increased CO2, thereby influencing the net atmospheric heating itself, and hence the so-called fast response of P to CO2 (DeAngelis et al. 2016). While substantial progress has been made in understanding these mechanisms, their contributions to recent observed changes in ET remain uncertain. For example, the recent multiyear decline in global terrestrial ET has been ascribed to both a persistent reduction in soil water availability (Jung et al. 2010) and to a transition to El Niño conditions (Miralles et al. 2014), with a very recent recovery in ET attributed to vegetation greening (Zhang et al. 2015).
d. PET and changes in aridity
The message from the research described above is that the mechanisms driving future changes in P and ET are complex, and are not easily partitioned into contributions from background variability, mean state changes, and changes in variability. In contrast, most climate projections show a common signal of twenty-first-century increases in atmospheric moisture demand and aridity in response to CO2-induced warming. These projections, however, often rely on aridity metrics based on potential evaporation (PE) (Burke and Brown 2008; Burke et al. 2006; Dai 2011a,b, 2013; Sheffield et al. 2012; Zhao and Dai 2015) or on potential evapotranspiration (PET) (Berg et al. 2016; Cook et al. 2014; Cook et al. 2015; Feng and Fu 2013; Fu and Feng 2014; McCabe and Wolock 2015; Scheff and Frierson 2014; Scheff and Frierson 2015; Swann et al. 2016) rather than on ET. PET represents the atmospheric evaporative demand from a well-watered reference grass. In contrast, PE is the maximum amount of evaporation that would occur from an open-water surface.
In most projection studies, PE and PET estimates assume a space- and time-invariant value for the stomatal resistance rs (Allen et al. 1998; Sheffield et al. 2012). While this is standard practice because it avoids problems associated with calibration, vegetation mixture, and morphology, this assumption does not account for the sensitivity of PET to plant adaptation in a CO2-enriched atmosphere. Recent scientific literature suggests that in simulations investigating the plant physiological response to elevated CO2, the choice of aridity index can have considerable influence on the estimated future drying, Drying can be enhanced by the use of aridity indices based primarily on atmospheric variables [such as surface relative humidity, the Palmer drought severity index (PDSI) or P/PET; Swann et al. 2016; Berg et al. 2016]. Alternatively, use of aridity indices based mainly on hydrological variables or vegetation properties (such as P minus ET, or runoff; Swann et al. 2016; Roderick et al. 2015) can reduce the estimated future drying.
Independently, Milly and Dunne (2016) show that PE and PET estimators based on time-invariant rs may severely overpredict the changes in aridity in non-water-stressed regions. The assumed space and time invariance of rs is therefore questionable. Whether this overprediction in aridity is systematic is unclear. Mixed results are found in studies comparing future drying trends over twenty-first-century North America obtained with either offline PET-derived PDSI or model soil moisture estimates (which implicitly account for the variations in rs). While the drying trends predicted with PDSI and near-surface soil moistures are similar (Cook et al. 2015), disparities across regions and models exist when PDSI and full-column soil moisture trends are compared (Feng et al. 2017; Smerdon et al. 2015). These mixed results highlight important uncertainties across different aridity predictors. The likely cause of such disparities is that we have relatively poor observational constraints on the impacts of CO2-induced structural and physiological changes in vegetation on soil moisture, runoff, and water resources (Frank et al. 2015; Ukkola et al. 2016).
e. The recent California drought
The 2011–16 drought in California can be used to study the relative importance of interannual and long-term changes in P, ET, and PET. The rainfall deficit over California has at least three proposed causes: 1) a purely internal mechanism initiated by equatorial Pacific SST fluctuations (Seager and Henderson 2016); 2) a human-induced, nonuniform change in geopotential height (Swain et al. 2016); and 3) a response to Arctic sea ice loss (Sewall 2005; Sewall and Sloan 2004). These explanations are not mutually exclusive. Other California drought studies based on PDSI indices require a contribution from human-induced increases in evaporative demand (Diffenbaugh et al. 2015; Williams et al. 2015). Although progressive forest canopy water loss has been observed over 2011–15 (Asner et al. 2016), the role of vegetation in explaining the California drought is presently unclear.
Despite these scientific uncertainties, the 2015/16 “super El Niño”3 (as it was referred to in media headlines) was expected to provide much-needed drought relief in California.4 Instead, the state received only near-average P during the event, which resulted in a worsening of the drought (as of December 2016, roughly 73% of California was still experiencing drought conditions). In contrast, the weak La Niña of 2016/17 has unexpectedly been associated with near-record rainfall and snowfall, bringing drought relief throughout California by April 2017. These two unusual years merit more scientific scrutiny. As they clearly illustrate, the expectation that a particular ENSO phase will almost always be a regional source of relief is significantly challenged by the limited predictability of the “far field” responses to ENSO.
In this study, we use different sets of atmospheric model simulations to identify regions where projected changes in aridity are large enough to overwhelm the effect of local drying/moistening associated with ENSO variability. To provide a comprehensive assessment of relative contributions to future drought, we consider a variety of hydroclimatic variables, including changes in P, PET, and ET, and several aridity indices. The latter are based on soil moisture, runoff, or on combinations of P and PET or P and ET. For each variable, the combined effects arising from mean-state changes and ENSO variability are examined for three components of greenhouse warming: 1) the so-called slow ocean warming (+SST), 2) the ocean warming and the response of plant physiology to rising atmospheric CO2 levels (+SST+Veg), and 3) ocean warming and plant physiology effects, together with the radiative forcing from enhanced CO2 (+SST+Veg+Rad). These sets of simulations are used to estimates the individual contributions of different factors and processes to future droughts.
a. Description of experiments
We use results from six different sets of numerical experiments (Table 1) available in the CMIP5 database (Taylor et al. 2012). Our objectives are as follows: 1) to identify regions where a projected change in aridity exists; 2) to determine whether this change exceeds the effect of local drying/moistening associated with ENSO variability; 3) to examine the robustness of results to model uncertainties; 4) to assess the sensitivity of results to the choice of aridity indices; and 5) to determine the relative contributions of ocean warming, plant physiology, and CO2-induced radiative forcing to future changes in aridity. We focus on ENSO and the ENSO season (January–March) for two reasons: ENSO is the primary contributor to internal variability in P, and ENSO-induced P effects have significant global socioeconomic impacts.
In estimating these relative contributions, we consider the incremental inclusion of 1) SST warming of 4 K (+SST), 2) the fast carbon cycle response associated with the impact on plant physiology of an imposed global increase of CO2 concentration (+Veg), and 3) the fast radiative responses to rising CO2 (+Rad). Atmosphere-only model experiments amip and amipFuture from CMIP5 are used to address objectives 1–4. Three other types of atmosphere-only simulations from CMIP5 (amip4xCO2, sstClim, and sstClim4xCO2) are analyzed to address objective 5. Although these five experiments are nominally a part of the Cloud Feedback Model Intercomparison Project (CFMIP), they have also been very useful for studying changes in the hydrological cycle (Bony et al. 2013; Chadwick 2016; Chadwick et al. 2014; DeAngelis et al. 2016).
Brief details of these simulations are as follows:
The amip experiment (herein Hist) is forced by observed time-varying SST and sea ice changes over 1979–2008, and the CO2-equivalent concentration is set at 345 ppmv (Gates et al. 1999).
The amipFuture experiment (herein Hist+SST) is almost identical to the amip simulation, with the exception that a composite SST warming pattern ( appendix A) scaled to a global mean of 4 K is added to the time-varying SST observations over the 1979–2008 period. Atmospheric constituents are held fixed at their amip values.
The amip4xCO2 experiment (herein Hist+Rad) is almost identical to the amip simulation; the sole difference is that the atmospheric CO2 concentration “seen” by the radiative scheme is quadrupled. The vegetation, however, only experiences the 1 × CO2 value.
The sstClim experiment (herein PI) is driven by prescribed climatological preindustrial SSTs, which vary over the seasonal cycle, but not interannually. Sea ice and atmospheric constituents are taken from each model’s own preindustrial control run. The simulation is 30-yr long; the radiation and vegetation schemes use the preindustrial value of equivalent CO2.
The sstClim4xCO2 experiment (herein PI+Veg+Rad) is almost identical to the sstClim simulation, with the exception that atmospheric CO2 concentration is quadrupled. Both the radiative and the vegetation schemes “see” the increase in CO2.
In all simulations, the dynamic vegetation schemes are turned off (i.e., the spatial distribution of vegetation is fixed). The vegetation, however, can respond to climate variations by altering leaf area, modifying stomatal resistance, and modifying productivity. We estimate the slow ocean-warming response by computing the amipFuture − amip anomalies:
We estimate the fast CO2-induced plant responses by computing the (sstClim4xCO2 − sstClim) − (amip4xCO2 − amip) anomalies:
The first part of this operation (sstClim4xCO2 − sstClim) allows the combined effect of vegetation and radiative responses to elevated CO2 to be inferred:
A total of four models (HadGEM2-ES, IPSL-CM5A-LR, MPI-ESM-LR, and MPI-ESM-MR) in the CMIP5 archive (Table 2) have contributed simulations and output variables suitable for addressing our objectives (with one single realization per model). Prior to analysis, all model data were regridded to a common 72 × 144 latitude–longitude grid using an area-conserving scheme. The significant advantages of using atmosphere-only simulations instead of coupled ocean–atmosphere simulations are described in appendix A.
b. Analysis of the amip and amipFuture experiments
Our analysis begins with the amip and amipFuture experiments, which isolate the effects of the change in mean climate and the variability induced by the slow ocean warming alone (+SST). We use model composites to identify the spatial patterns of land surface climate response to the two different phases of ENSO. Model composite response patterns are calculated both before to and after applying the SST warming. For the historical (Hist) climate, and for each variable of interest, we first identify the region R1 that is sensitive to ENSO. The R1 region is defined with two simple criteria. It encompasses continental grid cells that 1) show a sign change in anomalies between the two ENSO phases and 2) have a minimum temporal correlation of |0.15| with the model Niño-3.4 SST time series. The region R1 is expressed as an area-weighted fraction of the total land area between 60°S and 75°N (this excludes the large glaciated landmasses of Greenland and Antarctica). We then estimate how this fraction will be altered by a future change in the mean state. For further details see appendix B.
c. Analysis of the amip4xCO2, sstClim, and sstClim4xCO2 experiments
We add the +Veg and +Veg+Rad anomalies to the Hist+SST data to test the sensitivity of our results to the incremental inclusion of the fast physiological and radiative responses to imposed quadrupling of CO2. This adjustment involves adding the local monthly climatological anomaly of +Veg or +Veg+Rad to each month on each grid cell of the Hist+SST output. The +Veg anomalies are added to estimate the additional responses imposed by changes in CO2 due to plant physiology, yielding Hist+SST+Veg. The +Veg+Rad anomalies are added to provide information on the role of both radiative and physiological responses, yielding Hist+SST+Veg+Rad.
d. Aridity indices
We measure terrestrial aridity using six indices: the climate moisture index (CMI; appendix C), P/(P + PET), PDSI, surface soil moisture (10-cm depth), total soil moisture (summed over all soil layers, with a full depth ranging between 3 and 7 m across models; Table 2), and total runoff. CMI (Willmott and Feddema 1992) is a dimensionless drought indicator of the potential water availability, often used by the World Bank to assess the risks of climate change. The index P/(P + PET) is a modified version of P/PET, the aridity index adopted by the United Nations Environment Programme (UNEP) (Middleton and Thomas 1992), but it is better suited for analyzing future aridity in a warming world (Scheff and Frierson 2015). PDSI (Palmer 1965) is a widely used operational meteorological drought index based on both water supply and demand. The near-surface and full-column soil moisture variables are indicators of hydrological and agricultural droughts that incorporate different root zone depths. The two soil domains interact differently with the atmosphere (e.g., deeper soil moisture pool interacts mainly through plant transpiration and has a longer memory than the shallow soil layer, which is exposed to bare soil evaporation). Runoff, a primary concern for water managers, is directly related to terrestrial hydrological conditions; it can be affected by both CO2-induced structural (“greening”) and physiological change in vegetation (e.g., Ukkola et al. 2016).
Here, soil moisture and runoff are direct outputs from the AMIP-style simulations and depend on actual ET calculated in the model. In contrast, CMI, P/(P + PET) and PDSI are indices derived from both P and PET, where P denotes the water supply (a direct output of the simulations) and PET represents the atmospheric demand. PET is not supplied as model output and needs to be estimated. For this purpose, we used the Penman–Monteith equation ( appendix C), as recommended by the United Nations Food and Agriculture Organization (Allen et al. 1998). As in many drought projection studies (see section 1), our calculation of PET employs the CO2-invariant rs of a reference grass. The implications of increasing rs and reducing PET in a CO2-enriched atmosphere (by an increase in the denominator of the Penman–Monteith equation) will be addressed further in section 3f.
a. Changes in T in response to ocean warming
The left and right panels of Fig. 1a show the present-day land surface temperature T responses (or teleconnections) to El Niño and La Niña events (ΔHh),5 respectively, derived from the multimodel composite analysis of the Hist simulations. Only regions where the temperature anomaly changes sign between the two ENSO phases are displayed (region R1). This region represents 60% of the continental surface between 60°S and 75°N. Similar teleconnections patterns are found in the future (Hist+SST) in response to the ocean warming [Fig. 1b; ΔFf (see again footnote 5)]. Note that in both Figs. 1a and 1b, only region R1 from the Hist simulations is displayed. Figure 1c shows the corresponding temperature responses to El Niño and La Niña events in the Hist+SST simulations expressed relative to the present-day mean state (ΔFh). By the end of the twenty-first century, warming dwarfs temperature changes associated with ENSO. The phase of ENSO modulates the amplitude of warming everywhere but does not reverse the sign of the mean changes in temperature. In the southwestern United States, for example, the slight cooling during El Niño events does not fully compensate for the warming induced by the large increase in SSTs.
b. Changes in P and PET in response to ocean warming
From the Hist simulations, we find that only 41% of the continental surface between 60°S and 75°N exhibits P responses that change sign between the two ENSO phases (Fig. 2a and Table 3, first column). Meteorological droughts (P deficits) associated with El Niño events occur during boreal winter in western Australia, India, Indonesia, Brazil, southeastern Africa, Central America, and the northwestern United States. In contrast, La Niña events are characterized by drier than normal P conditions along the west coast of tropical South America, and at subtropical latitudes of North America and South America. Regions with positive T anomalies (Fig. 1a) are often collocated with regions experiencing negative P anomalies (Fig. 2a). Model expectations regarding the future behavior of precipitation are spatially more complex than with temperature, which generally shows a strong latitudinal dependence in all models. Estimates of the ENSO-driven P responses in the future relative to the present-day mean state (ΔFh; Fig. 2c) suggest that most regions north of 30°N will become wetter, independent of the phase of ENSO. In contrast, many subtropical regions will become drier, exceeding the range of ENSO-induced variability in P. In the future climate, only 20% of the original region R1 continues to show P anomalies that change sign with the phase of ENSO.
To highlight the results from ENSO response maps for all variables and experiments,6 we display them in a more compact graphical and tabular form in Figs. 3 and 4 and in Table 3. Note that Figs. 3 and 4 and Table 3 focus solely on the R1 regions, which are currently sensitive to ENSO.
The two left-hand bars of Fig. 3a summarize the precipitation results from Figs. 2a,c (the two right-hand bars are discussed in section 3d). The R1 region for P comprises the 41% of the land area with P responses changing sign during the two ENSO phases in the Hist experiments. During historical El Niño events, 56% of the R1 region is drier (yellow) and 44% is wetter (light green), as indicated by the bar labeled Historical. The +SST bar in Fig. 3a indicates that only one-fifth (20%) of the R1 region continues to be sensitive to ENSO in response to imposed future SST warming (combined light green and yellow portions of the bar; Table 3, sixth column). Nearly 41% (33%–44% across models) of the R1 region becomes wetter (dark green) independent of the phase of ENSO. Drying (dark brown) exceeds the range of ENSO-driven P variability in 39% (30%–44% across models) of the R1 region.
Figure 3b repeats the analysis in Fig. 3a, but shows results for PET instead of P. In the Hist simulations, the sign of PET anomalies is affected by ENSO variability over 53.5% of the global land surface. Roughly 54% of the R1 area has higher PET during El Niño events. The PET anomalies (see Fig. S1a in the supplemental material) correspond closely to the ENSO-driven T anomalies (Fig. 1a): anomalously warm (cool) regions yield anomalously high (low) PET.7 In the western United States, for example, warmer than normal conditions during La Niña events induce anomalously high PET and low P; both contribute to drier conditions over the region. In response to +SST, virtually all of the R1 land (99%) experiences increased PET, exceeding by far the PET changes induced by ENSO variability. The increase in PET is mainly due to the simulated increase in water vapor deficit, as is documented in a number of different studies (e.g., Roderick et al. 2015; Fu and Feng 2014; Scheff and Frierson 2014).
c. Changes in aridity indices in response to ocean warming
We now repeat the above-described T and P analyses using five of the six aridity indices: CMI, P/(P + PET), PDSI, surface soil moisture, and total soil moisture (see Table 3 and Figs. 3c–g, two left-hand bars). Across these indices, the percentage of the global land surface area with values for aridity sensitive to the phase of ENSO ranges from 39% to 50%. From 57% to 60% of these regions sensitive to ENSO are drier during El Niño events (Hist). This partitioning is consistent with results from Miralles et al. (2014), who found global-scale El Niño–induced limitations in terrestrial moisture supply. Ocean-forced warming (+SST) yields a pronounced decrease in the percentages of region R1 that remain influenced by ENSO variability. Increases in aridity overwhelm ENSO-driven variability in 67%–72% of the R1 regions (48%–72% when all models are considered individually), while only 9%–17% of the R1 region experience a decrease in aridity (9%–29% when all models are considered individually). Results are relatively insensitive to the choice of aridity index or to model uncertainties (Figs. 3c–h): preliminary conclusions indicate that both model surface and total soil moisture responses are largely consistent with the three offline, PET-derived aridity responses. Total runoff is the least affected aridity indicator, with 41% (27%) of the “ENSO sensitive” land area experiencing less (more) runoff in response to +SST.
d. Sensitivity of results to inclusion of physiological and radiative forcings
The above-described changes in aridity in the amipFuture simulations (+SST) occur in response to an imposed ocean warming and in the absence of increasing atmospheric CO2. In reality, however, regional hydroclimate responds not only to CO2-induced SST increases, but also to the influence of enhanced CO2 levels on the vertical profile of radiative heating and to altered surface energy fluxes arising from plant responses to increased CO2. Including the responses to both the radiative and physiological forcings can therefore lead to substantial changes in projected aridity: CO2 radiative heating enhances land surface temperature and increases evaporation and precipitation, whereas CO2 fertilization causes plant stomata to close, increasing rs and hence reducing terrestrial evapotranspiration (Andrews et al. 2011).
The simulations analyzed here enable us to assess the sensitivity of our results to the individual and combined radiative and physiological effects of increased CO2 (section 2c). Results for P and PET are synthesized in Figs. 3a and 3b, respectively. In the columns marked +SST+Veg, the effect of enhanced water-use efficiency of plants is added to the ocean warming effect. The combined effects of SST warming, changes in water-use efficiency, and atmospheric radiative forcing are shown in the columns marked +SST+Veg+Rad.
Consider the P results first. Adding the physiological effect of increased CO2 (+SST+Veg) to ocean warming only (+SST) has minimal impact on the total land area that changes sign during ENSO phases (Fig. 3a, yellow and light green; also compare Figs. 2d,e to Fig. 2c): there is a very small increase in the percentage of land area that receives consistently less precipitation (from 39% to 41%), but this change is smaller in amplitude than the inherent intermodel uncertainties (denoted by the error bars; all models simulate an increase in R1 except MPI-ESM-MR). This small reduction in precipitation caused by the incorporation of CO2-induced plant physiological effects (Fig. 2f) corroborates previous findings (Andrews et al. 2011; Berg et al. 2016; Swann et al. 2016; Cao et al. 2012). Globally, the P reduction caused by the physiological effects is more than offset by the P increase associated with the radiatively induced fast surface warming and increase in atmospheric water vapor (Fig. 3a; also cf. Figs. 2f and 2g). In the case of PET (Fig. 3b; see also Figs. S1c–e), results are insensitive to the inclusion of the physiological and radiative forcings. In all three sets of simulations (+SST, +SST+Veg, and +SST+Veg+Rad), most regions experience an increase in future PET that exceeds ENSO-driven variability.
In terms of PET-derived CMI, P/(P + PET), and PDSI, our results show a very small increase in permanently arid area in response to the inclusion of plant physiological forcing. For CMI, P/(P + PET), and PDSI, the permanently arid area is 69%, 70%, and 67% in +SST, respectively, increasing to 70%, 71%, and 68% in +SST+Veg, respectively (see Figs. 3c–e and 5; see also Figs. S2c,d and S3c,d in the supplemental material). This increase in aridity is in accord with the P/PET and PDSI results of Berg et al. (2016) and Swann et al. (2016), who attributed it primarily to increasing temperature and PET. PET increases even further with the inclusion of radiative forcing alone (e.g., Swann et al. 2016), leading to an additional increase of the proportion of land becoming arid beyond the range of responses induced by ENSO variability [to 78%, 79%, and 75% for the CMI, P/(P + PET), and PDSI indices, respectively, in +SST+Veg+Rad; see Fig. 5 and Figs. S2e and S3e]. With the exception of MPI-ESM-LR, all models simulate the increase in permanently arid area for all three indicators between the +SST and +SST+Veg cases. This increase in arid land in response to the inclusion of the radiative forcing (+SST+Veg+Rad) is robust across all four models considered and the three aridity indicators (Table S1 in the supplemental material). In summary, these three model-averaged PET-derived indices imply that a state of permanent aridity will occur in roughly 75%–79% of the regions where aridity is currently sensitive to ENSO variability.
e. Sensitivity of surface and total soil moisture changes to inclusion of physiological and radiative responses
In contrast to the case of the three PET-derived aridity indices, the inclusion of plant physiological effects has a pronounced impact on total soil moisture (Figs. 3g and 6). The percentage of ENSO-sensitive land area that is permanently arid in response to the imposed 4-K warming decreases from 72% in +SST to only 41% in +SST+Veg, a result robust in sign across models.8 In contrast, the soil moisture in the top 10 cm [Fig. 3f; see also Fig. S4 in the supplemental material and Dirmeyer et al. (2016)], which is primarily controlled by soil evaporation rather than by plant transpiration, is much less sensitive to the changes in stomatal closure than to the evaporative demand of the atmosphere. This explains why the shallow soil moisture responses are very similar to those obtained with the PET-based aridity indices (cf. Fig. 3f with Figs. 3c–e). In contrast, total soil moisture (Fig. 3g) reflects the deeper root zone soil moisture (the rooting depth for many ecosystems is estimated between 1.5 and 3 m; Hagemann and Stacke 2015), which is directly controlled by ET. The latter decreases during nondrought conditions, sequestering water in the soil, allowing plants to bridge longer dry spells. Owing to this additional soil moisture, runoff increases in many areas (Fig. 3h and Fig. S6 in the supplemental material).
To understand these results in more detail, we analyzed additional terrestrial variables: actual ET and its three components: evaporation from soils (EVs), evaporation from water intercepted by plants (EVv), and transpiration (TR). We also examined plant carbon uptake [i.e., gross primary productivity (GPP)] and leaf area index (LAI).
Under present-day conditions, water-limited regions with negative P anomalies during El Niño experience negative anomalies in ET, GPP, and LAI due to deficits in water availability for plant growth and evapotranspiration. No ENSO-related effects on GPP, LAI, and TR occur in regions where growth is limited by temperature, such as the boreal forest (see Fig. 1a of Nemani et al. 2003). In simulations of future warming, inclusion of the plant response to enhanced CO2 levels (from +SST to +SST+Veg) generally yields 1) larger land area with more carbon uptake (Fig. 4b, red shading, and Figs. S7c,d in the supplemental material) and 2) less transpiration (Fig. 4f, dark brown shading, and Figs. S8c,d in the supplemental material) and less ET per unit leaf area associated with increased stomatal resistance. The enhanced carbon uptake also favors more greening (Fig. 4c, red) and more canopy layers that can intercept (EVv; Fig. 4e, dark brown) and transpire water. Plant transpiration and evaporation (Figs. 4e,f) react more strongly than soil evaporation (Fig. 4d) to inclusion of physiological effects (+Veg). The relative magnitudes of changes in soil evaporation (Fig. 4d) and surface soil moisture (Fig. 3f) are in good accord; the same holds for the relative sizes of changes in transpiration (Fig. 4f) and total soil moisture (Fig. 3g). These results support the above-described discussion of the primary drivers of soil moisture changes. The changes in actual evapotranspiration (Fig. 4a) are influenced by changes in both soil evaporation and plant transpiration.
f. Implications of fixing stomatal resistance in PET-derived aridity estimates
In most studies predicting an increase in aridity under climate change (including this study), the PET-derived aridity metrics assume a stomatal conductance rs that is invariant to changes in CO2. The lack of rs output in the AMIP experimental database precludes the use of CO2-dependent values of rs in the PET calculations. As an alternative, we analyze two aridity indices based on actual ET instead of PET: CMI ET and P/(P + ET). The reason for choosing these ET-based aridity indices is that ET calculations account for variations in rs in response to solar radiation, temperature, and humidity, as well as in the response to CO2 levels imposed by the experimental design. Although not ideal, we use these alternative aridity indices in an attempt to bound the problem by assuming that the “true” projection of terrestrial aridity should lie somewhere between the PET-based indices (the upper bound) and the ET-based indices (the lower bound).9
In the +SST case (Figs. 4g,h), the ET-defined indices imply that more than half of the R1 region will become more arid in the future, beyond the expected range of ENSO variability [55%–57% for CMI ET and P/(P + ET)], with little change when both the physiological forcing and the fast response to radiative forcing are included [54%–57% for CMI ET and P/(P + ET)]. The regions predicted to face more arid climate in the future remain widespread, although these “lower bounds,” as defined by the ET indices, show smaller changes in aridity than those inferred from the “upper bounds” defined by the PET indices [78%–79% for CMI and P/(P + PET)]. Similar results between PET- and ET-derived indices are found when models are analyzed one by one.
4. Limitations of the study and future directions
The results presented above have several key caveats. First, our analyses do not consider the impact of possible changes in ENSO (in terms of amplitude, frequency, or ocean patterns in the tropical Pacific) or how such changes might affect the quantitative partitioning of the different drought influences that we have examined. At our present stage of scientific understanding, there are still significant uncertainties in model projections of future changes in ENSO magnitude and structure (e.g., Bonfils et al. 2015). In view of these uncertainties, our underlying assumption that ENSO characteristics will remain unchanged in the future is not unreasonable.
Another underlying assumption is that the regions defined as ENSO sensitive under present-day climate do not change with global warming. There is, however, evidence of an eastward shift of the teleconnection patterns between ENSO and hydroclimate variables in response to ocean warming (Zhou et al. 2014). This suggests that the locations of regions categorized as ENSO sensitive are not identical in the present-day and future climate simulations analyzed here (e.g., Figs. 2a,b show opposite El Niño–driven P signals in Mexico and in the southern United States). However, we do not expect this small locational shift to substantially alter our results. Although we do not consider the eastward shift in teleconnections identified by Zhou et al. (2014), we do explore an issue they did not examine: the impact of physiological and radiative forcings, which can also impact the strength or position of ENSO teleconnections.
A drawback of the simulations analyzed here is that the slow +SST response did not interact with the fast +Rad and +Veg responses. Additionally, the simulations did not consider the effects of dynamic vegetation or land use or land cover changes, and changes in the nitrogen cycle.
5. Discussion and conclusions
A number of studies have investigated changes in the mean state of hydroclimatic variables in response to global warming (Allen and Ingram 2002; Held and Soden 2006; Manabe et al. 1981; Marvel and Bonfils 2013; Wetherald 2010; Wetherald and Manabe 2002). The size and physical significance of projected changes in hydrological cycle components should be viewed in the context of an understanding of their intrinsic variability (Wetherald 2010). Such context provides insights into the potential risks of droughts in a warmer climate.
Understanding the likely sign and magnitude of projected changes in regional aridity indices is a challenging scientific problem. In this study, we consider 1) the relative contributions from changes in P versus evaporative demand, 2) the combined effects of mean-state changes versus variability associated with ENSO, and 3) the sensitivity of results to the “slow” ocean warming, the feedbacks from physiological adaptation of plants to enhanced CO2, and the “fast” radiative forcing changes caused by elevated CO2 concentrations.
Using simulations of present-day climate, we identified land surface areas where aridity indices currently change sign between the El Niño and La Niña phases of ENSO. This “ENSO sensitive” region is calculated for five different meteorological and agricultural aridity indices (three PET-dependent aridity indices and two soil moisture indicators). Depending on the index examined, 57%–60% of the ENSO-sensitive area in the Hist simulation experiences drying during the El Niño phase, whereas the other 40%–43% of the area experiences wetting during the same phase.
Both PET and P anomalies can drive changes in terrestrial aridity. In response to the slow global-scale ocean warming alone, we find that changes in meteorological and agricultural aridity are dominated by the increase in atmospheric demand, and that approximately 67%–72% of the ENSO-sensitive region becomes permanently arid. This increase in permanently arid area is not sensitive to the choice of the aridity index employed, and is robust across models. In these regions, even the large ENSO-induced intrinsic variability will no longer completely offset the expected regional-scale warming and drying signals associated with the ocean warming. In contrast, the simulated future changes in total runoff indicate that the projection of increasing hydrological droughts (with 41% of the ENSO-sensitive region) is less dire than that of increasing meteorological and agricultural droughts.
CO2 forcing alters the surface energy balance, reduces plant ET, increases land surface warming, enhances evaporative demand, and reduces the increase in P that is expected to result from global warming. According to the recent literature, the net impact of this physiological effect depends on the type of aridity index considered. Inclusion of plant physiology induces more drying in atmosphere-based aridity indices (such as relative and specific humidity, P/PET, or PDSI; Andrews et al. 2011; Berg et al. 2016; Milly and Dunne 2016; Swann et al. 2016), but less drying when aridity indices are considered that incorporate information directly relevant to plant growth (e.g., soil moisture, runoff, or P minus ET, Swann et al. 2016; see footnote 5). Over North America, aridity trends over the twenty-first century estimated with PET-derived PDSI and near-surface soil moisture are in good agreement, but disparities exist when compared to trends in deeper soil moisture (Cook et al. 2016; Cook et al. 2015; Smerdon et al. 2015). The uncertainties in the treatment of land processes can contribute to such disparities. For instance, transpiration may decline with a CO2-induced stomatal closure, but increase with greening and a lengthening of growing season (Ukkola et al. 2016). The impact of such factors on future aridity is not well understood.
Our results are in agreement with this recent literature: when using the three PET-dependent aridity indices dominated by the behavior of meteorological variables [CMI, P/(P + PET), and PDSI] and 10-cm soil moisture (layer that interacts closely with the atmosphere), we obtained a relatively small increase in permanently arid land area in response to the inclusion of plant physiology feedbacks. The permanently arid area in the +SST simulations (67%–70%) increased by only several percent in +SST+Veg (to 68%–71%). The area affected rises to 75%–79% in the +SST+Veg+Rad case, a result consistent across all models and the three aridity indices. Regions where the sign of local response associated with ENSO variability continue to change despite the change in mean state typically include the northwestern and southeastern United States and the region at the southeastern border of the Amazon rain forest (the boundaries of these regions are, however, spatially sensitive to the choice of aridity index).
We obtained a large reduction in future arid land when total soil moisture (a more plant-based index) is considered (from 72% in +SST to 41% in +SST+Veg). This large reduction mainly arises because the inclusion of CO2–physiology feedback preserves soil moisture in the root zone. While the differences in future projections of aridity between PET-based indices and total soil moisture have a plausible physical explanation, neglecting the increase in stomatal resistance in PET estimators may have led to an overprediction of the future drying estimated with the PET-derived indices. Our results highlight the importance of simulating the biological components of the Earth system in climate change studies, and clearly illustrate that the response of terrestrial aridity to greenhouse warming can be strongly dependent (in both sign and magnitude) on the choice of aridity index.
Model-specific rs information is not currently used in calculating PET-based aridity indices. There is therefore an inconsistency in the projections of aridity between the model’s actual ET-based estimate and the PET-based estimate that is approximated from the available model standard output (using the FAO’s spatially and temporally invariant value of rs). To reduce biases in terrestrial projections of aridity estimated with PET-driven aridity indices, we recommend that modeling groups provide calculated PET as a diagnostic output of the climate model (e.g., using model values of rs, wind speed, and surface net radiation). This new standard output would represent the simulated atmospheric water demand with greater fidelity than is presently possible in CMIP multimodel studies.
Sound decision-making requires climate change projections that are based on an understanding of both intrinsic variability and the likely changes in the climatic mean state associated with global warming. This is particularly important in developing climate change adaptation strategies, in assessing agricultural and infrastructure risks, and in managing water resources and land use. In some regions, ENSO variability will continue to be the major determinant of year-to-year changes in aridity. Other regions are likely to experience twenty-first-century changes in mean state that overwhelm the effects of intrinsic variability. The analysis performed here provides a sound scientific basis for identifying regions where societies may be more vulnerable to future climate change, and are likely to face drought conditions well outside of the historical range.
This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. GA and IC were supported by the DOE-BER Early Career Research Program award; CB, BDS, TJP, KT, MZ, and PD were supported by the Regional and Global Climate Modeling Program of the Office of Science of the U.S. Department of Energy (DOE); KM by DE-FOA-0001036 funding; and BI Cook by the NASA Modeling, Analysis, and Prediction program. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups for producing and making available their model output. For CMIP the U.S. Department of Energy’s PCMDI provides coordinating support and leads development of software infrastructure in partnership with the Global Organization for Earth System Science Portals. California information was obtained from https://www.drought.gov/drought/dews/california.
Advantages of Using Atmospheric Simulations
The use of atmospheric simulations presents certain diagnostic advantages. First, while most coupled climate models consistently project an increase in the frequency of extreme ENSO events in the future climate (Cai et al. 2014; Capotondi 2015), there is considerable model disagreement in projected twenty-first-century changes in the spatial structure and amplitude of ENSO (Bonfils et al. 2015; Coelho and Goddard 2009; Collins et al. 2010; Kao and Yu 2009; Vecchi and Wittenberg 2010). By prescribing the same sequence of observed SSTs in both the amip and amipFuture simulations, we assume identical ENSO evolution, and thus avoid introducing coupled model errors in the pattern and amplitude of ENSO. Second, coupled model biases in zonal P and SST features can obscure robust physical changes in hydroclimate (Marvel and Bonfils 2013). The AMIP-type simulations alleviate this problem by imposing the same SST changes in all models. This eliminates P errors arising from model-specific biases in SST magnitude and pattern. The use of the ensemble average of these simulations focuses attention on climate change features that are robust across models, and serves to reduce the noise associated with atmospheric internal variability.
The SST warming pattern (Fig. 1b of Zhou et al. 2014) is derived from a normalized composite of 13 individual CMIP3 model SST response patterns, obtained from simulations with a 1% yr−1 increase in atmospheric CO2. The response pattern was estimated at the time of CO2 quadrupling. Each model’s SST response was divided by its global mean SST change and multiplied by 4 prior to computing the multimodel ensemble mean, ensuring a global-mean SST forcing of 4 K (http://cfmip.metoffice.com/CMIP5.html).
There may be one minor inconsistency in the experimental design: A quadrupling of CO2 concentrations from the AMIP historical value (345 ppm) corresponds to a radiative forcing pathway leading to 8.5 W m−2 and 1380 ppm CO2 by 2100 (van Vuuren et al. 2011). Based on Fig. 12.5 of the IPCC Fifth Assessment Report, this radiative forcing is expected to induce a global annual-mean surface warming of about 4 K by about 2090. It is not clear whether a SST perturbation of 4 K by 2100 in the amipFuture simulation is slightly larger than expected. This possible inconsistency should not, however, alter the main conclusions of our study presented.
The main analysis consists of the following five steps:
For each model and each amip (Hist) simulation, we compute local January–March (JFM) anomalies relative to the climatological seasonal mean of that run (ΔHh). A similar calculation is performed for each amipFuture (Hist+SST) simulation (ΔFf).
We also calculate each model’s future JFM anomaly relative to its respective historical amip climatological value (ΔFh). (Note: for PDSI, the amipFuture PDSI values are already calculated relative to the amip baseline in the PDSI code.)
From the ΔHh, ΔFf, and ΔFh anomalies, we calculate the multimodel El Niño and La Niña composites, using nine El Niño events (1982/83, 1986/87, 1987/88, 1991/92, 1994/95, 1997/98, 2002/03, 2004/05, and 2006/07) and eight La Niña events (1983/84, 1984/85, 1988/89, 1998/99, 1999/2000, 2000/01, 2005/06, and 2007/08).
We then identify all continental grid cells where the variable of interest is influenced by ENSO variability (the R1 region). Inclusion in the R1 region requires that two criteria are fulfilled: first, that the correlation between the hydroclimate variable at this location and the model Niño-3.4 SST time series exceeds a minimum threshold of |0.15|, and second, that the local sign of ΔHh in the El Niño and La Niña composites changes between these two phases of ENSO. For example, a grid cell would qualify as “sensitive to ENSO” if it satisfies the correlation criteria with the Niño-3.4 index and also has land surface temperatures that are anomalously warm during El Niño and anomalously cold during La Niña. These criteria allow us to investigate whether the projected changes in mean state are large enough to overwhelm the sign of local response to ENSO variability. These criteria also provide a uniform basis for comparison across all models, simulations, and variables.
Finally, we identify all continental grid cells where the composited ΔFh anomalies continue to change sign between the two phases of ENSO (R2 region). Interpretation of these results is relatively straightforward: while both the R1 and R2 regions are primarily dominated by ENSO variability, the climate in the nonoverlapping portion of R1 and R2 (i.e., the land area present in R1 but excluded from R2) exhibits a change in mean state that overwhelms the range of ENSO-driven variability.
We note that steps 1, 3, and 4 closely follow the procedure of Zhou et al. (2014). By exploring the differences between maps of ΔHh and ΔFf, they investigated how large-scale warming may affect the variance in P associated with El Niño–driven teleconnections. In contrast, our goal is to focus on the El Niño and La Niña ΔFh composite maps, which provide information on the contributions of mean-state change and ENSO variability.
CMI and Penman–Monteith Potential Evapotranspiration
The CMI ranges from −1 to +1, with wet climates indicated by positive values:
The PET is not supplied as a model output and thus needs to be calculated. We use the Penman–Monteith equation (Allen et al. 1998) for a well-watered reference crop of 0.12-m height:
Here Rn is the net radiative energy available at the surface (MJ m−2 day−1) and G is the heat flux into the soil, which is close to zero on daily or longer time scales. The term es − ea is the vapor pressure deficit of the air (kPa), where es is the saturation vapor pressure (kPa) at air temperature Ta. Also, Δ = des/dT at Ta is the local slope of the saturation vapor pressure temperature relationship (kPa K−1), ea = (RHes)/100 is the actual vapor pressure (kPa), and RH (%) is relative humidity. The variable u2 is the wind speed at 2 m above the canopy (m s−1) fixed at 1 m s−1. The psychrometric constant is γ = 0.000 665Ps (kPa °C−1), with Ps the atmospheric surface pressure (kPa). Finally, Rn, Ta, RH, and Ps are all model outputs.
The factor u2/208 = 1/ra, where ra is associated with the friction of air flowing over vegetative surfaces, and rs is the bulk stomatal resistance that regulates vapor flow through the stomata openings of a reference crop. This reference vegetation is a 0.12-m-high green grass, with fixed rs = 70 s m−1, and an albedo of 0.23. For assessing future changes in aridity in a warming world, the Penman–Monteith equation is more reliable than other empirically based PET calculations, such as the Thornthwaite equation (Dai 2011b; Sheffield et al. 2012; Wehner et al. 2011), which is excessively sensitive to surface warming. Although using a constant value of rs in the Penman–Monteith equation is unrealistic under increasing CO2 levels, this is standard practice in many drought projection studies.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-17-0005.s1.
The actual amount of water removed from the surface through soil land canopy evaporation and transpiration.
Increasing CO2 can lead to increased photosynthesis but to similar transpiration in regions where plants are limited by water availability. Where plants are limited by energy availability, increased CO2 may cause little or no increase in photosynthesis but a decrease in transpiration (Goudriaan and Unsworth 1990). In both cases, there is an increase in WUE.
The oceanic Niño index (ONI; http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml) computed over the Niño-3.4 region indicates that super El Niño events have occurred in 1982/83 (ONI = 2.1), 1997/98 (ONI = 2.3), and 2015/16 (ONI = 2.2).
Most seasonal forecast models predicted a wetter-than-average California winter in 2015/16.
The notation ΔHh should be interpreted as follows: Δ indicates that this is a change in some variable (e.g., temperature and precipitation), H indicates that the quantity considered is from the historical period, and h indicates that the historical mean value of that quantity is removed to obtain an anomaly. We will present other combinations, such as ΔFf and ΔFh, where F and f refer to the simulations of idealized future conditions, respectively.
This close T–PET coupling is well known in drought research and is attributed to the use of PET formulations that are overly sensitive to temperature (Sheffield et al. 2012), such as the Thornthwaite equation. In our analysis, we find close T–PET coupling even when the Penman–Monteith equation is used for estimating PET (cf. Roderick et al. 2015).
However, this change is much smaller in amplitude in IPSL-CM5A-LR than in other models. Figure S5 in the supplemental material shows that IPSL-CM5A-LR, as expected, conserves soil moisture by stomatal closure. In this particular model, the effect is smaller than in other models and is overwhelmed by the increase in aridity induced by ocean warming (+SST case). The smaller preservation of soil moisture in IPSL-CM5-LR (relative to the other three models) is probably attributable to differences in land models, soil structure, root zone depth, vegetation mixture, or morphology.
Because ET is limited by both P and available soil moisture, using ET in aridity metrics is not ideal. For example, P/(P + ET) or P − ET may increase (i.e., indicating “wetter” conditions) because ET declines faster than P. If the decline in ET reflects a reduction in atmospheric demand, the conditions are indeed wetter. However, if the strong reduction in ET is caused by a sharp decline in soil moisture (indicating drier conditions), soil moisture and P/(P + ET) or P − ET results are inconsistent.