A calibration-free, robust estimation of monthly land surface evapotranspiration rates for continental-scale hydrology

In the past several decades there have been continuous scientific and societal discussion about the consequences of climate change on the flora and fauna of the Earth and on the human society as a whole. There is a wide-based consensus that the global hydrologic cycle is one of the most important bio-physical processes shaping the future of the globe through its precipitation, evaporation and/or transpiration (called evapotranspiration (ET) when combined) components leading to changes in frequency, magnitude and/or duration of flooding, droughts, and severe weather. The physical process of evaporation (whether the source is the free water surface, bare soil or stomata of the plants) is especially important since without it there is no precipitation on Earth, and also, the unusually high latent heat of vaporization plays an important role in very effectively cooling the Earth's surface. The role of evaporation/ET is even more fundamental since without the well-known water-vapor amplification, CO₂ loading of the atmosphere would have much less serious potential consequences because water vapor is a significantly more potent greenhouse gas than CO₂ itself, and to make it worse, it cannot be regulated, having about two-thirds of the Earth's surface constantly covered with oceans. Considering the vital role of ET in regulating local and global climate, it is somewhat surprising that meteorologists, climatologists, even hydrologists, only estimate it at a regional/continental/global scale, and depending on the approach, sometimes in a surprisingly crude way (i.e. when estimated mean annual land ET rates exceed precipitation rates at a regional scale without any obvious physical reason). This is so because large-scale direct measurement of ET is still not available, it can only be derived from properties of the evaporating surface (including its energy balance), its vegetation (for land surfaces) and the overpassing air.

One of the most important applications of estimated ET rates can be found in weather and climate modeling that require as input, the sensible and latent heat fluxes, calculated by so-called Land Surface Models (LSMs), such as employed by, for example, the European Center for Medium-Range Weather Forecast (ECMWF 2007) or by the National Centers for Environmental Prediction (NCEP), the latter producing the North American Regional Reanalysis (NARR) data set (Mesinger et al. 2006) that contains the so-derived ET rates.

McMahon et al. (2013), in an authoritative study, evaluated existing ET estimation methods that rely solely on widely available standard meteorological data and thus having the distinct advantage of going back in time no other methods, for example satellite-measurement based, can match, thus making it possible to detect any long-term changes in the hydrologic cycle possibly connected to, for example climate, and/or land use/cover change. Their conclusion was that the so-called complementary relationship (CR) based ET estimation methods were the most reliable ones for obtaining actual ET rates at a monthly time-step.

The CR of evaporation (Bouchet 1963) is based on the idea that at a certain time-scale, typically several days or longer, the lower atmosphere will adjust to the moisture status of the land surface (Crago et al. 2017) through the exchange of latent heat (i.e. ET) at the land-atmosphere interface, linking the two together. Note that at a daily time step any weather front passing the study area can destroy this dynamic balance between the air and the underlying surface by bringing air masses to the area with a moisture content unrelated to that of the land, therefore daily application of the CR is generally not recommended (Morton et al. 1985). Here the land surface is assumed to be homogeneous at a regional scale, i.e. at a few km and onward, in horizontal distance. The CR has the unique ability of providing regional ET rates with an unsurpassed accuracy (McMahon et al. 2013) among methods that require only basic atmospheric (air temperature, humidity, wind) and/or radiation data without knowledge of the land-use/land-cover type, surface temperature or soil-moisture status.

Different versions of the CR have been employed in a wide range of hydrological studies (see e.g. Morton 1975; Byrne et al. 1988; Granger & Gray 1990; Hobbins et al. 2001; Ozdogan & Salvucci 2004; Xu & Chen 2005; Xu & Singh 2005; Brutsaert 2006; Yang et al. 2011; Gao et al. 2012; Haj el Tahir et al. 2012; Ma et al. 2015) since its inception by Bouchet (1963). The recent generalization of the CR by Brutsaert (2015) has sparked a rush of renewed interest in the method, as demonstrated lately by Crago et al. (2016, 2017), Ma & Zhang (2017), and Szilagyi et al. (2016, 2017). For a background on the CR and snapshots of its important evolutionary steps through the past decades see Brutsaert & Stricker (1979), Morton (1983), Ramirez et al. (2005), Kahler & Brutsaert (2006), Huntington et al. (2011), Han et al. (2012), Szilagyi (2014, 2015), Brutsaert (2015), Crago et al. (2016, 2017), and Szilagyi et al. (2016, 2017).

Recently, Szilagyi et al. (2017), building on the latest developments in CR research (Brutsaert 2015; Crago et al. 2016), published a calibration-free version of the generalized CR and tested it with 30-year normals (1981–2010) of the input variables at a monthly basis. They claimed that their formulation would yield ET estimates on a par with current LSM outputs, such as employed in NARR, but this claim was never substantiated by the authors. Also, it is not clear if this calibration-free version of the CR indeed yields superior estimates to an already existing calibration-free version of the CR by Morton et al. (1985), called the WREVAP model. Finally, it is not known either if the model-independent calculation of the Priestley–Taylor (PT) parameter (Priestley & Taylor 1972) value proposed by Szilagyi et al. (2017) would similarly result in optimal ET estimates during a continuous simulation of the monthly ET rates, such as presented here. All these unanswered questions motivated the present work, but perhaps the strongest justification for it is to demonstrate that the Szilagyi et al. (2017) calibration-free formulation of the CR yields ET estimates superior to those produced by the LSM of NCEP and published in NARR, thus in the future it could help the calibration and verification of such LSMs employed by, for example, NCEP and ECMWF, and in doing so making them more accurate via an improved parameterization which would result in better weather predictions and ultimately in more realistic climate change scenarios.

The WREVAP ET rates can be obtained with the help of the original FORTRAN code of Morton et al. (1985), downloadable (together with help files and examples) from this author's personal website (http://snr.unl.edu/szilagyi/szilagyi.htm) or via a MATLAB script version of the FORTRAN code from the supplementary material of McMahon et al. (2013). An attractive feature of WREVAP, besides being calibration-free, is that it does not require wind data, however, it does require the mean annual precipitation rate, unlike other formulations of the CR approach.

The 32-km resolution monthly surface net radiation and wind speed values of NARR as well as its LSM ET values for 1979–2015 can be downloaded from the NOAA website (www.esrl.noaa.gov/psd/data/gridded/data.narr.html). See Mesinger et al. (2006) for more information of their version of a LSM employed by NARR.

Equation (1) was applied in a continuous monthly time step for the 37-year period of 1979–2015 over the conterminous United States employing the 4-km spatial resolution Parameter-Elevation Regressions on Independent Slopes Model (PRISM, Daly et al. 1994) air and dew-point temperature data. The 32-km NARR surface net radiation and 10-m wind data were linearly interpolated onto the PRISM grid employing a power transformation (Brutsaert 1982) of the 10-m wind (u₁₀) values into u₂ = u₁₀(2/10)^1/7, required by (2) and (5).

Note that the value of the α parameter should depend on the actual data employed (i.e. one cannot automatically expect to have the same α value between cases derived from variables differing in spatial resolution and/or method of measurements, for example as in the case for PRISM and NARR data and any arbitrary combinations of them) but not on the CR model chosen as the above procedure is completely model independent.

The ET estimates were validated by PRISM precipitation (P) and United States Geological Survey (USGS) Hydrologic Unit Code level-6 (HUC6) runoff (Q) data over the conterminous United States. The conterminous USA has a diverse climate ranging from desert/semi-desert (in the south-west) to Mediterranean (California) to oceanic (north-west coast) to humid continental (the eastern half) to humid subtropical (Gulf region) to tropical (southern Florida). The mean annual precipitation varies accordingly from 100 to 300 mm (desert-semi-desert) to 1,000 mm (eastern half) to over 2,000 mm (north-west coast), with a spatial average of 790 mm (Figure 1). The 4-km spatial resolution PRISM precipitation data is considered as the best available gridded precipitation product over the contiguous United States. See Daly et al. (2008) for details on data quality.

Validation of the ET estimates involved a simplified water-balance (Morton 1983; Ramirez et al. 2005; Szilagyi 2015; Szilagyi et al. 2016, 2017) for the long-term mean annual value of the watershed-averaged ET rate, ET_b, as P – Q for each HUC6 watershed (Figure 2(a) and 2(b)), employing 37-year averages of the monthly P (spatially averaged over the catchment) and annual Q values, against which the estimates of (1) could be compared. There were three catchments (identified in red in Figure 2(a)), out of the altogether 334 HUC6 watersheds where the USGS runoff values were unreliable and presented themselves as clear outliers not only in the regression plot of the ET vs ET_b values, but among the Q values of the neighboring watersheds, thus the ET_b values of these catchments were overwritten by the estimates of (1).

The simplified water balance assumes that the mean annual value of the change in stored water volume, ΔS, over the watershed (in the absence of trends in S) is negligible when the averaging period is long. This is certainly true for a stationary process of zero mean since the sample mean approximates the theoretical mean value with increasing length of the sampling period, as the basic rule of statistics (e.g. Brockwell & Davis 1987). To illustrate this point, mean annual values of ΔS were calculated for the available 15 years (2002–2016) of water storage (S) data provided by the Gravity Recovery and Climate Experiment (GRACE) at a 0.5 degree spatial resolution (Tapley et al. 2004) over the conterminous USA, accessible from ftp://podaac.jpl.nasa.gov/allData/tellus/L3/mascon/RL05/JPL/CRI/netcdf/.

For obtaining annual ΔS values, the original S values derived for certain days of the year were interpolated by a spline method for each day of the 2002–2016 period. The annual changes, ΔS, were calculated by taking the difference in the daily S values for the first day (i.e. January 1) of each consecutive year. A mean annual ΔS value was obtained for each GRACE cell by averaging the resulting 15 ΔS values. Figure 3 displays a histogram of the mean annual GRACE-derived ΔS values. The spatial average of the mean annual ΔS values is indeed very close to zero, −0.66 mm, with a standard deviation of a mere 8.8 mm, which means (by assuming a normal distribution) that about 95% of the mean annual ΔS values are within −18 and 18 mm, which is indeed negligible (less than 4%) in comparison with the mean annual ET value of 540 mm (Figure 4). Naturally, for the 37-year period of this study, instead of the 15-year period of GRACE data, this interval of ΔS values would be even smaller.

The long-term mean annual ET estimates of (1) for the 1979–2015 period are displayed in Figure 4 at the 4-km resolution of the PRISM data. ET rates are smallest along the California-Southern Nevada and Western Arizona border, respectively, dropping to the lowest value of about 50 mm yr⁻¹ near Death Valley, California. The highest values, in excess of 1,300 mm yr⁻¹, are found along the coast of the Gulf of Mexico, due to the presence of lagoons, marshy areas, as well as the available high net surface radiation fluxes. When compared with the distribution of mean annual PRISM precipitation rates (Figure 1), the effect of large-scale irrigation (Szilagyi et al. 2011) in the most heavily irrigated state, i.e. Nebraska (third state north of Texas), is clearly discernible along the Platte Valley, running west-to-east in the lower third of the state, as a westward protruding bulge of elevated ET rates, denoted by greenish to light blue colors in Figure 4, in opposition of the eastward shift of the region of decreased precipitation (marked by blue colors in Figure 1) as one moves from south to north.

Figure 5 depicts the long-term mean ET to precipitation ratios. A value larger than unity can only be expected over water-rich regions of wetlands and lakes (see for example the Great Salt Lake in Utah and Lake Okeechobee in Florida), or over heavily irrigated extensive areas in an otherwise dry climate where this ratio is already close to unity, before the effects of irrigation (Szilagyi et al. 2017) is accounted for. Such areas are found in almost every state of the West. This does not mean that all regions with a larger than unity value in Figure 5 are such water-rich areas. The model tends to overestimate ET rates in the basins of mountainous regions with varied topography. The exact reason for this overestimation is not known. Possible explanations may include the following. (a) In these basins air humidity is enhanced by evaporation from the surrounding high-elevation areas and mountain sides where precipitation rates are more abundant thus creating a partial disconnect between air humidity and the underlying land in the basin, somewhat similar to coastal areas of sea-land air circulation. (b) In most of these basins irrigation agriculture is significant (Brown & Pervez 2014) leading to elevated ET rates in excess of precipitation. (c) In these typically high-elevation mountain basins a significant portion of precipitation derives from snow. Unshielded rain-gauges are known to have catch-deficiencies as high as 90% (Dingman 2015), meaning that only about 10% of the snow-rate is measured, especially in windy conditions. (d) A combination of all these factors. Overall, ET rates in excess of P are only found in less than 11% of the PRISM cells.

The spatial distribution of the long-term mean annual (1)-derived ET to ET_b ratios is displayed in Figure 6(a). In about 80% of the catchments (Figure 6(b)) the (1)-derived ET estimates are within 20% (turquoise or green colors) of the water-balance derived ET_b values. The largest underestimation takes place in Arizona (white and pink colors), while the largest overestimation (brown color) is found in Washington State, the desert region of California, and over the Great Lakes. In the latter case the ET_b values represent ET rates found on the islands of the lakes, where air humidity and the resulting ET estimate are greatly influenced by the surrounding open water surface, resulting in significant overestimation of the actual island ET rates. In desert climates, on the other hand, where precipitation rates are low (i.e. around 100–150 mm yr⁻¹ in Southern California, see Figure 1), even a 50–60 mm yr⁻¹ overestimation of the actual ET rates may yield an ET to ET_b ratio in excess of 1.6, as seen in Figure 6(a).

The (1)-derived HUC6-averaged mean annual ET estimates are plotted against the simplified water balance derived values, ET_b, in Figure 7. The ET estimates of (1) scatter nearly perfectly around the 1:1 line, with a best-fit line of ET = 0.97ET_b + 22, quite remarkable from a calibration-free method, having an R² value of 0.87, RMSE of 89 mm yr⁻¹ and a mean absolute error of only about 5 mm yr⁻¹ (Table 1). As a result, these estimates outperform earlier calibrated versions of the CR by Brutsaert (2015), Szilagyi (2015), Szilagyi et al. (2016), and Crago et al. (2016) using 30-year monthly normals (1981–2010) of the input variables across the contiguous USA.

In comparison, the WREVAP ET rates significantly overestimate the water-balance derived long-term mean annual values by 124 mm yr⁻¹ on average and yield a best-fit line slope of only 0.88 (Figure 8(a) and 8(b)). The overestimation becomes even more pronounced (i.e. 146 mm yr⁻¹) for the NARR LSM values (Figure 9(a) and 9(b)). See Table 1 for the performance measures of the three ET estimates.

The CR approach of regional evaporation works poorly near sudden discontinuities of the moisture status of the environment (e.g. near shorelines or large-scale irrigation in an otherwise dry environment) where the moisture content of the air may be significantly decoupled from that of the underlying surface (Morton 1983), as seen in the largest ET to P ratios of Figure 5 along the Pacific Coast of southern California. Where such discontinuities are absent, however, the current calibration-free formulation of the CR works remarkably well in comparison with earlier versions of the CR approach and a widely available LSM ET product.

The CR of evaporation is a useful practical tool in land surface ET mapping over extensive areas as it requires only basic, typically widely available meteorological/radiation data often with a long temporal coverage, such as the here employed PRISM and NARR databases. The current calibration-free version is best suited for a continental scale as its model-independent estimation of the PT parameter value requires periodically or permanently wet areas. The present method can play a crucial role in water-balance investigations where large-scale modeling of latent heat fluxes is of prime importance and where a water-balance based calibration/validation is not possible due to missing or unreliable streamflow and/or precipitation data. It could also become a potentially useful tool for independent validation and/or calibration of LSM-predicted latent heat fluxes (over land areas) that form the backbone of climate modelling, see for example the NARR reanalysis product. To the best knowledge of the author, no such attempt has yet been made. The expected improvement in LSM-derived ET fluxes could eventually lead to better global circulation and/or climate models that would lead us to a deeper understanding of the earth-atmosphere-climate system and eventually to improved climate-change scenarios.

The author is grateful to the Editor and two reviewers for their valuable comments that greatly improved the original manuscript.