Potential evapotranspiration (PET) is a supply-independent measure of the evaporative demand of a terrestrial climate—of basic importance in climatology, hydrology, and agriculture. Future increases in PET from greenhouse warming are often cited as key drivers of global trends toward drought and aridity. The present work computes recent and “business as usual” future Penman–Monteith PET fields at 3-hourly resolution in 13 modern global climate models. The percentage change in local annual-mean PET over the upcoming century is almost always positive, modally low double-digit in magnitude, usually increasing with latitude, yet quite divergent between models.
These patterns are understood as follows. In every model, the global field of PET percentage change is found to be dominated by the direct, positive effects of constant-relative-humidity warming (via increasing vapor deficit and increasing Clausius–Clapeyron slope). This direct-warming term accurately scales as the PET-weighted (warm-season daytime) local warming, times 5%–6% °C−1 (related to the Clausius–Clapeyron equation), times an analytic factor ranging from about 0.25 in warm climates to 0.75 in cold climates, plus a small correction. With warming of several degrees, this product is of low double-digit magnitude, and the strong temperature dependence gives the latitude dependence. Similarly, the intermodel spread in the amount of warming gives most of the spread in this term. Additional spread in the total change comes from strong disagreement on radiation, relative humidity, and wind speed changes, which make smaller yet substantial contributions to the full PET percentage change fields.
a. Why potential evapotranspiration?
Potential evapotranspiration (PET), a basic land climate variable (e.g., Hartmann 1994), is the rate at which a given climate is trying to evaporate water from the soil–vegetation system. In other words, for given atmospheric and radiative conditions, PET is the surface evapotranspiration (ET) rate that would hold if the soil and vegetation were well watered. Synonymous and near-synonymous concepts include reference evapotranspiration, potential evaporation, evaporative demand, and pan evaporation. Critically, PET may be thought of as the water required to maintain a garden or irrigated crop, or the water “price” a plant must pay to maintain open stomata. A higher-PET climate is thus a more arid, evaporative climate. Therefore, in this study we attempt to understand how local PET will scale with global greenhouse warming, using global climate models (GCMs) as well as basic physical principles.
PET is also of interest because it is a key factor explaining other hydrologic and climatic quantities. Several prominent conceptual models of land hydrology, including the Palmer Drought Severity Index (PDSI) (Palmer 1965) and the Budyko and Miller (1974) ecohydrologic theory, take precipitation (water supply) and PET (water demand) as climate-supplied forcings and give soil moisture, actual ET (latent heat) flux, runoff, and/or drought index as land-generated responses. In these sorts of frameworks, understanding precipitation and PET changes is necessary for understanding other land hydroclimatic changes. In particular, recent studies using the PDSI to warn of widespread drought increases with future greenhouse warming (e.g., Dai 2013; Burke et al. 2006) cite systematic global PET increases as the main driver of their alarming results. Understanding the nature, magnitude, and pattern of these projected increases is the motive of the present work.
Additionally, PET is a more natural choice than actual ET for the evaporative component of land “aridity” metrics because changes in actual ET often just reflect supply (precipitation) changes. For example, the well-known study of Seager et al. (2007) uses precipitation minus actual ET (P − E) to quantify modeled aridification due to greenhouse warming in a subtropical terrestrial region where the model precipitation declines a great deal. The model ET (the Seager et al. E) in this area also significantly declines, not surprisingly. However, the analysis, by its nature, interprets the ET decline as if it is some other factor helping to offset or mitigate the precipitation decline. In fact, the model climate is probably becoming more evaporative, not less, owing to warming and (presumably) cloud-cover and relative-humidity reduction, and this should not mitigate but aggravate the local ecological effect of the precipitation reduction, even though the actual evaporative flux necessarily decreases owing to the supply decrease (e.g., Brutsaert and Parlange 1998). To avoid this type of pitfall, the aridity of a climate is usually quantified using the ratio P/PET of annual water supply to annual water demand, or similar (e.g., Budyko and Miller 1974; Middleton and Thomas 1997; Mortimore 2009), which has the additional advantage of being dimensionless. Then P/PET < 0.05 is defined as hyperarid, 0.05 < P/PET < 0.2 as arid, 0.2 < P/PET < 0.5 as semiarid, and so forth. Feng and Fu (2013) show that global climate models project systematic future decreases in P/PET (i.e., aridification) over most of the earth’s land, again owing to the (projected) systematic PET increases that we attempt to understand in this work.
b. Quantifying PET
Except where an evaporation pan, lysimeter, or other direct method is available, PET cannot be measured in the field, so it is usually estimated from its meteorological and/or radiative causes. Several estimation methods are in wide use. All of the above studies of greenhouse-driven future drought or aridity expansion (Dai 2013; Burke et al. 2006; Feng and Fu 2013) use the Penman–Monteith equation, a fundamental physics-based method (Penman 1948; Monteith 1981). Given some near-surface air temperature Ta, water vapor pressure ea, wind speed |u|, and net downward broadband radiation Rn, this equation simply gives the latent heat flux (LH) (equivalent to ET) that solves the system:
for the three unknowns SH (sensible heat flux), LH, and Ts (skin temperature that would hold under well-watered conditions). Here, (1) is the bulk formula for SH, (2) is the bulk formula for LH under well-watered conditions, (3) is the surface energy budget, e*(T) is the saturation vapor pressure at a given temperature T, ra = 1/(CH|u|) is the aerodynamic resistance between the canopy surface and the level where Ta and ea are measured, CH is a scalar transfer coefficient, rs is the bulk stomatal resistance under well-watered conditions, G is the heat flux into the ground or soil (usually parameterized or ignored), ρa is the air density, cp is the air specific heat, γ = (cpps)/(εLυ) is the collection of constants from having written (2) in a manner analogous to (1), ps is the surface air pressure, ε ≈ 0.622 is the ratio of molar masses of water vapor and dry air, and Lυ is the latent heat of vaporization of water.
The solution to this system proceeds by noting that, if Ts − Ta is not too large, then e*(Ts) ≈ e*(Ta) + de*/dT(Ta)(Ts − Ta), which allows Ts to be cleanly eliminated between (1) and (2), giving [with the help of (3)]
which is the surface latent heat flux that would hold under well-watered conditions with the given meteorology and radiation. Here Δ:= de*/dT(Ta) is the standard shorthand for the local slope of the Clausius–Clapeyron curve, which will be used from now on. By definition, this flux (divided by Lυ) is the potential evapotranspiration. The resulting Eq. (5) is the Penman–Monteith equation:
The first term in the numerator of (5) is known as the radiative term, and the second is called the aerodynamic term. Note that in the latter we have rewritten [e*(Ta) − ea] the vapor pressure deficit appearing in (4), as e*(Ta)(1 − RH), where RH is the near-surface relative humidity. This allows changes in the vapor pressure deficit to be separated into constant-RH changes in e* (from Ta changes), and constant-Ta changes in RH. [Henceforth we are dropping the (Ta) and simply writing e* for e*(Ta), since Ts has been eliminated.]
Many of the input variables in (5) will change with significant greenhouse warming. Most immediately, the surface net radiation Rn will tend to increase (absent any cloud feedbacks) because of the extra longwave emitters in the atmosphere, sending more longwave energy back at the surface. This alone would tend to increase PET (5). However, the warming itself will also directly change PET through e* and Δ, which both increase with Ta by the Clausius–Clapeyron law. Constant-RH increases in e* will increase PET by widening the vapor pressure deficit, especially where and when RH is low. [The discussion in the review paper of Roderick et al. (2009) omitted this mechanism.] Increases in Δ may increase or decrease PET depending on the magnitudes of various terms in (5). It is not clear a priori whether the radiation changes or these direct-warming changes will dominate.
In addition, RH might change in either direction, through a common theoretical expectation for RH is that it should remain roughly constant (e.g., Held and Soden 2000), as generally observed thus far (e.g., Held and Soden 2006). This is the main motivation for considering constant-RH e* changes separately from changes in RH.
Finally, raw observations indicate that |u| decreased in most land areas over the past several decades (McVicar et al. 2012), in sufficient magnitude to overcome the concurrent Ta increases in (5) and explain the widespread observations of declining pan evaporation, that is, declining PET (McVicar et al. 2012; Roderick et al. 2009; Wang et al. 2012). However, it is still unclear whether this terrestrial wind “stilling” is a measurement artifact, as it does not appear in reanalyses (e.g., Pryor et al. 2009; McVicar et al. 2008) or marine observations (McVicar et al. 2012), and some of the pan-evaporation declines themselves are also raw and unadjusted for observing-system changes. Even if real, it is highly unclear whether the stilling is due to global warming (McVicar et al. 2012), and it may have reversed course after about 1998 (Wang et al. 2012). Therefore, in this study we take the future model output of |u| at face value, which contains no such systematic declines. However, if any real global stilling trend of the proposed magnitude were to continue unabated into future decades, PET would presumably continue declining and the conclusions of our study (as well as those mentioned above) would not apply.
Other, non-Penman methods of estimating PET are also in use, as mentioned at the beginning of this subsection. The Thornthwaite (1948) method and other temperature-proxy methods empirically relate PET to Ta alone, for a given location and time of year. This simplicity has encouraged their frequent use for variability in the current climate (e.g., Palmer 1965), which has led some studies to use them, or models containing them, to assess future climate change (e.g., Wehner et al. 2011; Price and Rind 1994); see also references in Lofgren et al. (2011). However, within a given climate (especially during warm, high-PET parts of the year), anomalous warmth is associated with anomalous sunshine (higher Rn), and often also with anomalous low RH, significantly enlarging the positive response of (5). By contrast, future climate change should warm Ta without the sunshine and RH changes that might accompany a similar warm anomaly in year-to-year variability. Thus one would expect from (5) that an empirically determined dependence of PET on Ta from year-to-year variation would overestimate the greenhouse climate change response. Indeed, several studies (e.g., McKenney and Rosenberg 1993; Hobbins et al. 2008) have found that the same long-term climatic changes can imply large increases in the Thornthwaite PET but much smaller increases, or even slight decreases in Penman–Monteith PET. Similarly, negative PDSI responses to future global climate model output are 2–3 times stronger using the default Thornthwaite PET than using Penman–Monteith PET (A. Dai 2012, personal communication) Thus, studies that use a simple temperature-proxy method to assess future PET changes may be severely flawed.
Other studies of future climate change (Lofgren et al. 2011; Arora 2002) simply estimate PET as Rn/Lυ, which we will call the energy-only method. While this works reasonably well for spatial differences in the present climate (Budyko and Miller 1974), one would presume that it underestimates future PET increases because it does not include the independent physical effects of Ta through the Clausius–Clapeyron law, discussed above.
Still other studies of PET change in global climate models (e.g., Rind et al. 1990) have directly used an internal land-model field that is also called “potential evapotranspiration.” However, this field is (quite confusingly) not the same concept as what we have been discussing: it is what would instantaneously start evaporating if the surface were to be suddenly wettened, without any chance to cool down the skin temperature Ts and establish energy balance (3) with Rn. In other words, this field is directly computed using the bulk LH Eq. (2), where rs is still the well-watered “open” stomatal resistance but Ts is now the actual skin-temperature output of the model instead of the well-watered skin temperature used above, which is often much cooler. Indeed, Rind et al. (1990) (and references therein) found that this model “PET” achieved summertime climatological values averaged over the United States of ~40 mm day−1 in the climate models of their day. [The observed summertime PET maxima in, e.g., Hartmann (1994) are almost an order of magnitude lower.] So this quantity, while interesting perhaps, is not the object of our study (and also it is not publicly archived by any of the GCMs in the current phase of the Coupled Model Intercomparison Project).
Therefore, in this study we use (5) to quantify and understand the PET response to future greenhouse warming.
The Penman–Monteith potential evapotranspiration (5) is usually many times larger in magnitude during the day than at night because of both Rn and the vapor pressure deficit. Thus, daytime climate changes may affect time-integrated PET much more than nighttime climate changes, so it is desirable to examine diurnally resolved climate and PET. In the recent fifth phase of the Coupled Model Intercomparison Project (CMIP5) (Taylor et al. 2012), subdaily surface output is conveniently accessible for the first time, at 3-hourly resolution. Sixteen of the CMIP5 global climate models archive all of the necessary information (surface energy budget terms and near-surface temperature, moisture, and wind) at this resolution for years 2081–99 in the business-as-usual representative concentration pathway 8.5 (RCP8.5) scenario and 1981–99 in the historical scenario. However, in three of these, the meteorological fields are given at, say, 10 m above the soil surface, instead of 10 m above the canopy top (M. Watanabe 2013, personal communication), making them inapplicable to (1) and (2) [and thus (5)] in forest areas. So, we use the remaining 13 models, which we list in Table 1 along with any model-specific exceptions to our procedures. We use output from run 1 (“r1i1p1” in CMIP5 filenames) only.
A prominent version of (5) is the recent American Society of Civil Engineers (ASCE) standardized reference evapotranspiration equation (Allen et al. 2005), which was explicitly developed for the purpose of standardizing the computation of reference or potential ET for all users. The development included the systematic intercomparison and testing of numerous operationally used Penman-type methods. Our full method closely based on Allen et al. is given in section a of the appendix. Briefly, we fix CH and rs as universal constants corresponding to “alfalfa” values as specified by Allen et al. (2005), with CH ≈ 5.7 × 10−3, and rs varying between 30 s m−1 (day) and 200 s m−1 (night). (We will see in section 5 that our conclusions are not very sensitive to these vegetation parameters.) We also compute (Rn − G) as LH + SH (3) because the models do not output G, and we let e*, Δ, and ρa depend on Ta as specified in Allen et al.
Using these procedures and values, for each of the 13 CMIP5 models in Table 1 we compute Penman–Monteith PET (5) for every model grid cell that is at least 80% ice-free land and for every 3-h interval in the 19-yr epochs 1981–99 and 2081–99 (except those that fall on 29 February in models that use the full Gregorian calendar). For each interval in the calendar year, we average over the 19 years to obtain a diurnally and annually varying PET climatology of each epoch. Averaging over the calendar then gives annual-mean climatologies of PET. These are shown for 1981–99 in Fig. 1 along with their multimodel means and appear quite reasonable with higher modeled PET in sunnier, lower-RH, and/or warmer locations. As an additional reality check, Fig. 2 plots these against the corresponding model climatologies of actual ET; each dot is one grid cell. In almost all of the models, our computed PET is a fairly clean, efficient upper bound on the model’s actual ET, as expected from the definition. That is, the most well-watered model grid cells are actually evapotranspiring at rates quite close to our independently computed PET. This success is a rather pleasant surprise considering the very different origins of the two quantities, the models’ use of full Monin–Obukhov surface layer dynamics for CH, and the potentially large contrast between ASCE-standard alfalfa and the vegetation specified in the model grid cells.
3. Model results
a. Full PET change
For each of the 13 models and for the multimodel mean, Fig. 3 maps the raw percentage change in climatological annual-mean PET (5) between the 1981–99 and 2081–99 epochs. At each location PET always or almost always increases; that is, ambient conditions become more evaporative with greenhouse warming. This more careful calculation confirms the similar results of Burke et al. (2006), Dai (2013), and Feng and Fu (2013), who quantified this future PET increase only for the mean (or for a single model), and did not resolve the diurnal cycle. In some models a few largely high-latitude regions do see PET decreases or little change in PET, but these are quite localized, and even in these places most models (and the mean) show increases in PET.
Furthermore, the magnitude of the projected PET increases is usually in the low double digits of percent, on the order of 10%–45%. In many models certain northern and/or mountainous locations see more than this, but over very broad swaths of land these sorts of values are typical. For the multimodel mean, the first row of Table 2 summarizes this by averaging the percentage change values over various latitude bands (and subsequent rows similarly average subsequent figures). The magnitudes in Fig. 3 agree well with those in Feng and Fu (2013) despite the differing methods. They are also comparable to change magnitudes for annual precipitation P (e.g., Meehl et al. 2007). This further confirms the importance of using P/PET or similar diagnostics when thinking about the land aridity response to global warming, instead of just P (and/or actual evapotranspiration E, which often contains the same information as P as discussed in section 1a).
In most models and in the mean, there is also a clear tendency toward greater percentage increases in PET at higher latitudes (as alluded to above and seen in Table 2), and again Feng and Fu (2013) obtain a similar structure. As far as we know, this basic property has not been explicitly noted in the literature before. (We will see in section 4 that the main reason for this is not Arctic amplification of warming.)
Yet despite all of these broad commonalities, the models also disagree a great deal, on both the detailed spatial patterns and the overall magnitude. We will see how these disagreements arise from differences in the climate changes projected by the models.
b. PET changes due to individual factors
Figures 4, 5, 6, and 7 show the percent changes in climatological annual-mean PET (5) that result from perturbing (Rn − G), Ta, RH, and |u| one at a time to 2081–99 levels while keeping the other variables at 1981–99 levels, as explained in detail in section b of the appendix. One can immediately see here and in Table 2 that the always-positive PET change owing to the Ta increase (Fig. 5) dominates the other factors in most locations, and explains most of the overall 10%–45% magnitude in Fig. 3. This is why PET increases are so much more common than decreases. Again, the physical mechanisms here are widening of the vapor pressure deficit by constant-RH increases in e*, and lowering of the saturated Bowen ratio by increases in Δ (plus isobaric lowering of ρa to a small extent). RH also changes, but the resulting PET changes (Fig. 6) are of both signs, inconsistent from model to model, very weakly positive in the multimodel mean, and only sporadically (nowhere, in the mean) negative enough to cancel the Ta-induced increases in Fig. 5. This validates the constant-RH baseline idea and justifies our decision to think of the vapor pressure deficit as e*(1 − RH) rather than the more customary (e* − ea). (An alternative null assumption of constant vapor pressure deficit would imply systematically increasing RH, which we do not see.) However, the RH-driven changes can still be very important locally in some models, explaining the East African PET decrease in BNU-ESM in Fig. 3, for example.
PET changes owing to the surface energy supply (Rn − G) (Fig. 4) are also usually positive, confirming the physical intuition from section 1b. However, with modal values of less than 10% (e.g., Table 2) they are generally of secondary importance to the Clausius–Clapeyron-driven changes (Fig. 5) just described. This was not clear a priori—in fact, some studies in the literature had used radiation changes alone to infer PET changes, as discussed above in section 1b. As with RH, though, some models have localized regions where radiation-induced change becomes dominant, such as the Amazon Basin in MRI-CGCM3 (and several other models) or the Tibetan Plateau in INM-CM4 (cf. Figs. 3–5).
In contrast, PET responses to |u| changes (Fig. 7) are only rarely important compared to the other changes. In the multimodel mean and in some individual models (the two BCC models, CNRM-CM5, and INM-CM4), they are hardly noticeable, usually no larger than ±5%. Like the RH responses (Fig. 6) they have no strongly preferred sign, although decreases are perhaps slightly more common than increases. This is all in stark contrast to the dominant “stilling” role posited for |u| in the putative recent PET declines, discussed in section 1b.
Finally, subtracting the sum of these attributed pieces (Figs. 4–7) from the full PET change (Fig. 3) gives the residual PET change due to nonlinearities, covariance changes, and changes in neglected inputs such as ps. This residual is shown in Fig. 8 and is quite weak (0%–10%) compared to the Ta-driven or even (Rn − G)-driven changes, though it is usually positive. (The GFDL-CM3 residual at high northern latitudes is a major exception to both of these statements, perhaps because the changes there in Figs. 3–6 are all so large.) In any case, we can clearly claim success in our attribution exercise since the residuals are much smaller than the full changes in Fig. 3 and are close to zero for the multimodel mean.
Having now examined all of the pieces, we can see that the constant-RH PET response to temperature change (Fig. 5) not only explains the general positivity and low-double-digit magnitude of the full PET change, but is also largely responsible for the high-latitude amplification noted in the previous subsection. The response to (Rn − G) (Fig. 4) is also polar amplified, but the temperature response still seems to contain most of the latitudinal contrast shown in Fig. 3, as can be clearly seen in Table 2. As for the intermodel disagreement in PET change, responsibility seems to lie with almost all of the terms, but disagreement in the Ta-driven term alone is still large, especially in the overall magnitude. [This makes sense given the well-known disagreement between global climate models on the magnitude of warming in response to an emissions scenario, that is, transient climate sensitivity (e.g., Meehl et al. 2007).]
Therefore, we now attempt a detailed quantitative understanding of the structure and magnitude of this model PET response to ambient temperature change as depicted in Fig. 5.
4. Analytic scaling for the PET response to temperature
a. Basic idea
How, exactly, is Penman–Monteith PET (5) sensitive to Ta with all else constant? First, one can note that in the numerator of (5) both the aerodynamic term and the radiative term increase roughly like Clausius–Clapeyron (C-C) with Ta at constant RH because e*(T) is a roughly exponential function, so Δ := de*/dT has roughly the same fractional rate of increase with T as e* does. More precisely [and using the empirical form from Allen et al. (2005) and the appendix for consistency]:
where T is in degrees Celsius and e* in pascals. At earthlike temperatures de*/(e*dT) is around 6%–7% °C−1 but 2/(T + 237.3) is only 0.7%–0.8% °C−1, so (8) means that dΔ/(ΔdT) is not far from de*/(e*dT) at all. [These values still hold using the physical C-C equation in place of the empirical (7).] So, all else constant, we can expect the entire numerator of PET (5) to increase at a C-C-like rate with warming of Ta, regardless of the relative importance of the radiative and aerodynamic terms. This is why we did not further split the response to Ta into responses to e* and Δ in section 3b above.
However, the denominator of (5) cannot necessarily increase so fast: although Δ increases at roughly C-C as demonstrated, γ(1 + rsCH|u|) does not depend on Ta at all. This term stops the denominator from fractionally increasing as fast as the numerator and apparently is the key reason why PET always increases with Ta (Fig. 5) despite the ambiguous sign of the Δ-driven response discussed in section 1b. If not for the presence of γ(1 + rsCH|u|), the denominator would increase about as fast as the numerator, and PET might not be very sensitive to Ta at all.
b. Derivation and exposition of the scaling
To quantify all of this, we now take the relative partial derivative of (5) with respect to Ta, repeatedly using the rules
where fa := a/(a + b) and fb := b/(a + b), and
plus the chain rule, to yield
Here frad is the fraction of the numerator of (5) made up by the radiative term, as in (9). Similarly, faero is the fraction of the numerator made up by the aerodynamic term, and fΔ is the fraction of the denominator of (5) made up by Δ.
We then use (8) to write de*/(e*dT) in terms of dΔ/(ΔdT):
Using frad + faero = 1 and the ideal-gas-law formula for ρa, this reduces to
the main equation that we will use to understand the constant-RH PET response to Ta as depicted in Fig. 5.
The first term within the brackets in (13) tells the story laid out in the previous subsection: the numerator of PET (5) scales like C-C [dΔ/(ΔdT) · 1] or about 5–6% °C−1, but the denominator Δ + γ(1 + rsCH|u|) scales closer and closer to C-C the more important Δ is in it [−dΔ/(ΔdT) · fΔ], weakening the net response. Since Δ is an increasing function of Ta, this cancellation should occur more (fΔ should be larger and the denominator should be more C-C-like) in warmer base climates, so the percentage sensitivity of PET to Ta should be less in warmer base climates. We will see in section 4d that this explains the polar-amplified response pattern in Fig. 5. (Similarly, the sensitivity should be greater in windier climates, in which fΔ is reduced.)
The second term within the brackets in (13) contains the small, miscellaneous departures from the above: the 0.7%–0.8% °C−1 discrepancy between the scalings of e* and Δ in the numerator, and the −0.3%–0.4% °C−1 isobaric dependence of air density on temperature. The partial cancellation between these two effects makes the net even smaller, ~0.4% °C−1 at the most since faero can only range between 0 and 1. Therefore, from here on we define ϵ(T) := 2/(T + 237.3) − 1/(T + 273.15) and write
c. From instantaneous to annual-mean scaling
Our Eq. (14) may be a theory for PET sensitivity at a particular instant. However, the results from section 3 and Fig. 5 that we wish to understand are about annually averaged PET. So, to test (14) it is not immediately clear what inputs should be used. For example, we might use the annual-mean warming . (From here on, an overbar will denote the annual mean.) However, if the warming in some place is, say, 6°C at night but 2°C during the day, then using the mean value of 4°C will overestimate the response because the vast majority of PET is concentrated during the day when the warming is only 2°C. So, we need to carefully consider the scaling of the annual mean, , in addition to the instantaneous PET considered earlier in this section.
The relative change in turns out to be the PET-weighted average of the relative change in instantaneous PET. This is because, again, the more PET is concentrated at a particular time, the more a percentage change in PET at that time matters to the percentage change in . Mathematically,
where from here on a double overbar denotes a PET-weighted annual average, for any variable a. So, using (14),
d. Testing the annual-mean scaling
To test this scaling theory (17), we compute , , and so forth for each model grid cell. For the base-state variables PET, fΔ, and Ta, we just use diurnally and annually varying 1981–99 climatologies computed as in section 2, not the full 19-yr time series. Similarly, for the change dTa, we use the same smoothed diurnally and annually varying climatological difference that we used to produce Fig. 5, as detailed in section b of the appendix. (Note that turns out to simply be the fraction of annual-total PET that comes from the aerodynamic term, so we compute aerodynamic and radiative separately, and directly use this fraction.)
Since dΔ/(ΔdT) is not that dependent on temperature and is small, the main sensitivity “wild card” in (17) should be , the fraction of the denominator of (5) made up by Δ at high-PET times of day and year. The determines whether the denominator will keep up with the numerator’s Clausius–Clapeyron pace and curtail the PET increase with warming, or lag behind it and allow a large PET increase.
So, in Fig. 9 we map for each model, as well as the multimodel mean of (summarized in Table 2, as above). One can see that it dramatically varies from as low as ~0.25 in the cool-summer climates of the coastal high latitudes to ~0.75 in the warm climates of the tropics. Apparently the strong dependence of Δ on temperature is in control of this fraction, even though it also depends on quantities in the denominator’s other term (|u| and our day length–dependent rs). Indeed, Fig. 10 shows that the PET-weighted (i.e., daytime, warm-season) basic-state temperature has a strikingly similar spatial pattern to this , often even at very fine spatial scales. Essentially, in cool, low-Δ climates the denominator of (5) is mainly made up of γ(1 + rsCH|u|), which stays fixed with Ta and lets (5) increase, while in warm climates it is dominated by Δ, which scales like C-C and cancels most of the numerator’s attempt to increase PET.
Figure 11 then maps the entire bracketed term from (17), that is, our scaling estimate of the percentage sensitivity of to PET-weighted warming. As guessed, its pattern is nearly the same as that of (Fig. 9) [and thus (Fig. 10)], varying from around 1.5% °C−1 over large areas of the planet’s warm, high- tropics to nearly C-C in the coolest-summer regions where is small and the numerator in (5) can increase nearly unopposed. The models agree on all of these fields much more than they agree on the gross response to Ta change depicted in Fig. 5. This is not surprising since these are only based on properties of the modeled 1981–99 base climates, which can be tuned to match observations.
On the other hand, the PET-weighted projected warming [the other factor in (17)] might vary considerably from model to model since the models do not agree on the warming response to a given greenhouse-gas forcing scenario (e.g., Meehl et al. 2007). Figure 12 maps for each model and for the mean, confirming that the spread in modeled warming is much larger than the spread in estimated sensitivity to that warming (Fig. 11). Taking the end members, over land seems to be almost three times stronger in GFDL-CM3 than in INM-CM4! Thus it appears that the main reason for the intermodel spread in the magnitude of the change due to warming (Fig. 5), noted in section 3b, is indeed the intermodel spread in the warming itself.
We are also now in a position to evaluate the source of the high-latitude amplification of the percentage change pattern in Fig. 5, and thus in Fig. 3. Figure 12 shows that the PET-weighted warming is, indeed, strongly Arctic amplified in some models (e.g., BNU-ESM and GFDL-CM3). However, in many other models this pattern is absent, even though it is well known that the Arctic amplification of the annual-mean warming is robust across climate models (e.g., Meehl et al. 2007). For example, in ACCESS1.0 maximizes in midlatitude North America and Europe and in the Amazon Basin, and in GFDL-ESM2G and GFDL-ESM2M maximizes in the subtropics. So does not consistently show high-latitude amplification, and in the multimodel mean any such amplification is quite weak (Fig. 12; Table 2). This is probably because high-latitude warming amplification is more of a cold-season than a warm-season phenomenon (Meehl et al. 2007), while a PET-weighted mean is largely over the warm season. In contrast, the sensitivity factor in (17), depicted in Fig. 11, shows strong and systematic high-latitude amplification because of the strong control of (Fig. 9) by the basic-state temperature (Fig. 10), as discussed above. Thus it appears that (Figs. 5 and 3) is polar amplified not because the warming is polar amplified, but rather largely because colder climates with Δ less important in the denominator of PET are inherently more sensitive (cf. the last two lines of Table 2).
Finally, we can confirm this picture by evaluating (17) and comparing to the model responses to Ta changes in Fig. 5. Before displaying the result, we need to note that, if the sensitivity factor in (17) is, e.g., 4% °C−1 and the projected warming is 9°C, the expected change should be noticeably larger than 36% because 1.049 ≈ exp(0.36) > 1.36. To account for this simple nonlinearity, we exponentiate (17) and subtract 1 to arrive at our final scaling guess for what Fig. 5 should look like.
This estimate is shown in Fig. 13 and is strikingly close to the model response in Fig. 5. In fact, the summary values in Table 2 differ from the actual values on the line above by only about +1% (of the basic state, about 10% of the changes). Thus, we can claim success in understanding the magnitude, structure, and intermodel spread in Fig. 5. The low double-digit percent magnitude of comes from the mid-single-digit °C greenhouse warming (Fig. 12) times the sub-Clausius–Clapeyron, 1%–4.5% °C−1 sensitivity of (5) at constant RH (Fig. 11 and sections 4a,b). The structure of comes mainly from the structure of the base-climate temperature (Fig. 10) via (Fig. 9) and the sensitivity (and also somewhat from the structure of the warming). The intermodel spread comes from the intermodel spread in the warming.
5. Sensitivity of results to imposed vegetation
One might wonder whether the above holds for parameter choices in (5) other than the ones presented in section 2 and the appendix. In particular, the transfer coefficient CH and bulk stomatal resistance rs could potentially modulate the Ta-independent term γ(1 + rsCH|u|) in the denominator of (5), and therefore alter and the bracketed sensitivity in (17). So, we also compute results using a few alternative choices for these two parameters.
We first examine the effect of setting rs ≡ 0, that is, neglecting the relatively small but appreciable stomatal resistance of well-watered transpiring leaves, as in many formulations of Penman–Monteith PET including those used by Burke et al. (2006) and Dai (2013), as well as in the case of pan evaporation. This gives an expression more in the spirit of Penman (1948) than Monteith (1981): the denominator of (5) simply becomes Δ + γ. This choice should systematically increase fΔ and thus reduce the percentage change in [by (17)], taking it even further from Clausius–Clapeyron. Indeed, the range of shifts upward, to roughly 0.4–0.85 (not shown). However, the original range in Fig. 9 was about 0.25–0.75, so this is a quantitative but not a qualitative increase. The spatial pattern of hardly changes, except for losing some finescale structure owing to the loss of |u| dependence.
Figure 14 shows the percentage changes in PET from changing Ta in this case. Comparison with the analogous Fig. 5 shows that setting rs ≡ 0 indeed weakens the response, making single-digit- percentage values somewhat more common and values > 30% less common, but the patterns are very close. The at-a-point differences between the two figures are much less than the spatial and model-to-model variations within each figure, and the summary statistics in Table 2 differ by only about 2%–3% (of the basic state).
We also examine a “smooth” version of (5), in which the 0.5-m vegetation height h and thus the roughness lengths zom and zoh in (A1) are reduced by a factor of 10, setting h to a grasslike 5 cm and halving CH from ≈5.7 × 10−3 to ≈2.8 × 10−3. [The Penman–Monteith formulations used in Burke et al. (2006), Dai (2013), and Feng and Fu (2013) also assume a smoother surface.] This, too, shifts the range of only slightly upward, to roughly 0.35–0.8, with a very similar spatial pattern to the original in Fig. 9. So, the percentage change in ends up looking almost identical to Fig. 5, but slightly (several percent) weaker (not shown). Again, the parameter-induced alterations in are much less than the spatial and model-to-model variation.
Also, in the no-resistance case, adding this smooth vegetation would not appreciably lower the results any further because, in that case, CH does not even appear in the denominator of (5) and thus can no longer affect . Thus, the effects are not additive—the no-resistance case gives a strict upper bound on and an effective lower bound on warming-induced (Fig. 14).
Finally, we examine a “rough”, forestlike PET in which h, zom, and zoh are increased by a factor of 10, setting h to 5 m and tripling CH to ≈1.7 × 10−2, a very large value. In this case, the range of falls to roughly 0.15–0.7, again with a very similar spatial pattern to Fig. 9. Then from Ta becomes somewhat larger than shown in Fig. 5, with values of 35%–40% or more becoming more widespread. Again, there is little qualitative or pattern change; the overall story is the same. In summary, widely different choices of vegetation parameters do not alter the big picture presented in sections 3 and 4 above.
There is also the question of whether rs, like (Rn − G), Ta, RH, and |u| (and ps), should have been treated as changing between the two epochs rather than staying fixed. After all, the carbon dioxide increase that causes greenhouse warming may also cause individual plant stomata to close (e.g., Sellers et al. 1996). However, there is still very large uncertainty about the bulk vegetation changes that will occur in concert with this, much larger than the uncertainty in the climate response to carbon dioxide (Huntingford et al. 2013). Almost nothing is known about this bulk response. Furthermore, the percentage sensitivity of Penman–Monteith PET (5) to a percentage change in rs turns out to depend very strongly on the vegetation parameters rs and CH, in contrast to the much weaker dependence just presented in the case of sensitivity to Ta. Therefore, in this study we decided to only scale the PET response to climate change, and not the response to carbon-dioxide-induced plant physiological change.
6. Summary and discussion
Potential evapotranspiration (PET), the rate at which surface water evaporates if available in a given climate, has been projected to increase with future greenhouse warming in most or all locations, driving strong global trends toward drought (e.g., Dai 2013; Burke et al. 2006) and/or aridity (Feng and Fu 2013). In this study, we systematically analyzed the projected response of the Penman–Monteith equation (5), the fundamental physical quantification of PET used by those studies. We found that, at least in the 13 modern global climate models listed in Table 1, the main reason for the projected PET increase is the warming itself (Fig. 5), not the greenhouse-driven increase in surface net radiation (Fig. 4). The warming causes the PET increase by widening the vapor pressure deficit e*(1 − RH) corresponding to a given relative humidity RH, and/or by increasing the local slope Δ := de*/dT of the Clausius–Clapeyron curve, which governs the partitioning between sensible and latent heat fluxes. Changes in RH are not of any strongly preferred sign and are not large enough to alter this.
The magnitude of the projected annual-mean PET increase between 1981–99 and a business-as-usual 2081–99 scenario is usually a low double-digit percentage (Figs. 5 and 3, Table 2), comparable to projections for local precipitation. This is because the numerator of Eq. (5) increases like Clausius–Clapeyron [5–6% °C−1] with constant-RH warming, but in the denominator only the first term Δ increases similarly, while the second term stays fixed. Thus, the net response of (5) to warming is sub-Clausius–Clapeyron, usually about 1.5–4% °C−1 (Fig. 11). The higher values are found in cooler climates where Δ is smaller and thus less important in the denominator of (5) [i.e., in (17) is smaller], and the lower values are found in warmer climates, explaining the strongly polar-amplified change pattern. Since the projected PET-weighted-mean warming for this scenario tends to be in the single digits of °C in most places (Fig. 12), the gross percentage response of PET to warming ends up in the lower double digits (Figs. 5 and 13). Large disagreement between models on the exact amount of warming produces similar disagreement on the total PET response. (The smaller but appreciable radiation- and RH-driven PET change components shown in Figs. 4 and 6 also vary widely between models, adding to the disagreement.)
A key, further advantage of our scaling approach (17) is that a climate model is not even needed for a user to locally compute the sensitivity of PET to future warming. All variables within the square brackets in (17) can be computed during routine calculations of observed present-day Penman–Monteith PET. For example, the values of fΔ = Δ/[Δ + γ(1 + rsCH|u|)] and PET can be noted at each calculation time step and averaged over several years of data collection to obtain seasonally and/or diurnally resolved climatologies, which can then be used to find . If it turns out that can be accurately estimated straight from and , then the computation will be even simpler, as there will be no need to archive short-term values of fΔ. So, whether the sensitivities plotted in Fig. 11 contain model biases is not actually that important for the practical use of (17).
We would also like to briefly give a more qualitative, physical explanation for why PET is less sensitive to Ta in warmer base climates. First, consider a climate cold enough that LH is unimportant in (3), even under well-watered conditions, and the dominant balance is between SH and (Rn − G). In this climate, fixing (Rn − G) effectively fixes SH, which fixes (Ts − Ta) by (1). Now, if we rewrite (2) with the substitutions introduced later in section 1b,
we can see that LH will be able to increase with temperature at a Clausius–Clapeyron rate, driven by Δ and e*. Everything else in (18) is fixed by assumption. However, as the climate warms and well-watered LH becomes appreciable, the evaporation will start to cool Ts relative to Ta and limit the fractional increase of (18). [Eventually, energy conservation (3) will start to severely limit the increase in well-watered LH, since (Rn − G) is fixed here, and (Ts − Ta) and SH can only go so negative owing to constraints involving the wet-bulb depression associated with our fixed RH.]
We also note that the PET percentage responses to changes in (Rn − G), RH, and |u|, depicted in Figs. 4, 6, and 7, can also be analytically scaled in the manner demonstrated for Ta in section 4, with similar levels of success. However, the modeled changes in these variables (for input to these scalings) are not as well understood as the modeled warming dTa, so these scalings do not provide as much understanding.
Finally, we are still interested in under what conditions or assumptions this large systematic PET increase with climate warming actually implies a systematic drying out of the land, as suggested by much of the work cited in section 1. To this end, we also have work in progress testing the sensitivity of modeled soil moisture to large changes in global temperature across a very wide range of continental geographies, forcing mechanisms, and land and atmospheric modeling choices.
The lead author would like to thank A. Dai for providing results on drought projections using different PET methods, and A. Swann for a key suggestion on section 4b. The authors also acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table 1) for producing and making available their model output. For CMIP the U.S. 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 Portals. This work was supported by NSF Grants AGS-0846641 and AGS-0936059.
a. Parameter and procedural choices for the Penman–Monteith equation
Allen et al. (2005) provide parameters for two different reference vegetation types: short clipped grass, and alfalfa (with the expectation that “crop coefficients” will be determined for conversion of the resulting potential evapotranspiration output to values suitable for other vegetation). We use the alfalfa values, reasoning that natural vegetation is closer to alfalfa in roughness and leafiness than it is to short clipped grass. Similarly, procedures are standardized separately for hourly and for daily calculation time steps; we use the hourly procedures on the 3-hourly model intervals. For meteorological variables, the model output is given as synoptic “snapshots” every 3 h, so for each interval we average the initial and final values of Ta, specific humidity qa, and |u| to estimate 3-h means, analogous to the hour means used by Allen et al. [Note that the raw output includes the wind components u and υ but not the speed |u|, so |u| has to be computed as at each snapshot before this averaging step. Also, two of the models (see Table 1) give u and υ on a grid staggered by one-half the spacing in latitude and longitude from the main grid used for all the other variables, so we compute u at each main grid point as the mean of u at the four surrounding wind grid points, and similar for υ, before this computation of |u|.]
With these choices of time step and vegetation type, the ASCE standardized procedures for variables in (5), and our few departures from them, are given as follows. A constant ps is hydrostatically estimated from the surface elevation, but for simplicity we directly use the 3-hourly ps output from the model, averaged like Ta and qa above. The e* is computed from the (3)-h-mean Ta using the empirical form e*(T) = 610.8 exp[17.27T/(T + 237.3)], where e* is in pascals and T is in degrees Celsius, and Δ using its derivative. Also, ρa is computed from the dry-air ideal gas law, using the (3)-h-mean Ta multiplied by 1.01 to account for virtual effects. A number of standardized methods are given to compute ea from measurements; we directly use the 3-h-mean model qa above, multiplying by ps/ε to convert the units (nearly identical to their method 1). RH can then be computed as ea/e*. (In a few models this RH can occasionally slightly exceed 1, presumably due to interpolation; in these cases we set RH = 1 to avoid unphysical negative values of the aerodynamic term.) The Lυ is idealized as a constant 2.45 × 106 J kg−1. A field estimation method for Rn and a simple parameterization of G are given, but we simply compute (Rn − G) from the model-output actual turbulent heat fluxes SH and LH using (3), which is still valid. These fluxes are already provided as 3-h means over our intervals, so there is no need for averaging. Then rs is set at 30 s m−1 (“open”) during the day and 200 s m−1 (“closed”) at night, where “day” and “night” are defined as Rn > 0 and Rn < 0. We use (Rn − G) > 0 and (Rn − G) < 0 instead; this is justified since Allen et al. (2005) parameterize G as a small positive fraction of Rn.
For the transfer coefficient CH, the standardized choice is the neutral, log-layer form,
where k is the von Kármán constant, zw is the height of the wind speed measurements, zom is the momentum roughness length, zh is the height of the temperature and humidity measurements, zoh is the scalar roughness length, and d is the zero-plane displacement. Allen et al. do not attempt to justify this choice, but one could argue that the great majority of PET is in warmer seasons or climates during the daytime, when the surface layer is either neutral or convective. For most wind speeds the Monin–Obukhov correction to CH for convective conditions is much smaller than for stable conditions, so the worst of the potential problems are avoided. In any case, the standardized values for (A1) are as follows: k is set to 0.41, and zw and zh are each set to 2 m, although we use 10 m for zw to match the height of the model wind output. If h is the assumed vegetation height (0.5 m for our standard alfalfa choice), zom is set to 0.123h, and zoh to 0.0123h. Finally, d is set to 0.08 m on the assumption that the weather measurements are taken over clipped grass, but we conservatively set d = 0 as it is not clear what exactly the model output heights are measured relative to. With these choices, CH works out to ≈5.7 × 10−3.
b. Determining the PET responses to individual variables
We would like to isolate the PET changes owing to changes in the individual inputs (Rn − G), Ta, RH, and |u|. However, we cannot simply give (5) the 2081–99 time series for one of these and the 1981–99 time series for all other variables because the differing synoptic histories of the two epochs would destroy any interinput correlations other than the diurnal and annual cycles, adding an artificial change to the result. So, for each of these four inputs, we compute diurnally and annually varying climatologies for each model (as for PET), further smooth them with a 7-day running mean that respects the diurnal cycle, difference the two epochs (divide them, in the case of |u|), and perturb each year of the 1981–99 input time series by this diurnally and annually varying difference (factor), creating an input time series with the climatological properties of 2081–99 but the synoptic history of 1981–99. These can then be used one at a time in (5) to isolate the responses to (Rn − G), Ta, RH, and |u|. [When we perturb (Rn − G), we still use the original 1981–99 (Rn − G) series to define day and night for setting rs. Global warming may accomplish many feats, but it certainly will not transmute night into day! Consistent with this, when computing the 2081–99 PET in section 2, we subtract our diurnally and annually varying climatological difference from each year of the 2081–99 (Rn − G) series before it is used to define night and day.]