Abstract

Archived global measurements of water loss from evaporation pans constitute an important indirect measure of evaporative flux. Historical data from evaporation pans shows a decreasing trend over the last half century, but the relationship between pan evaporation and moisture-limited terrestrial evaporation is complex, leading to ambiguities in the interpretation of these data. Under energy-limited conditions, pan evaporation Epan and moisture-limited terrestrial evaporation E increase or decrease together, whereas in moisture-limited conditions these fluxes form a complementary relation (CR) in which increases in one rate accompany decreases in the other rate. This has led to debate about the meaning of the observed trends in the context of changing climate. Here a two-dimensional numerical model of a wet pan in a drying landscape is used to demonstrate that, over a wide range of realistic atmospheric and surface conditions, the influence that changes in E have on Epan 1) are complementary and linear, 2) do not depend upon surface wind speed, and 3) are strikingly asymmetrical, in that a unit decrease in E causes approximately a fivefold increase in Epan, as found in a recent analysis by Kahler and Brutsaert of daily evaporation from U.S. grasslands. Previous attempts to explain the CR have been based on one-dimensional diffusion and energy balance arguments, leading to analytic solutions based on Penman-type bulk difference equations. However, without acknowledging the spatially complex humidity and temperature fields around the pan and, specifically, how these fields change as the contrast between the wet pan and the drying land surface increases, such integrated bulk difference equations are a priori incomplete (they ignore important divergence terms), and thus these explanations must be considered physically incomplete. Results of this study improve the theoretical foundation of the CR, thus increasing the reliability with which it can be applied to estimate water balance and to understand the pan evaporation record of climate change.

1. Introduction

The global trend of decreasing pan evaporation Epan over the last half century has been interpreted as evidence for both a slowing (less terrestrial evapotranspiration E; Peterson et al. 1995) and accelerating (more E; Brutsaert 2006; Brutsaert and Parlange 1998) hydrologic cycle. Although decreases in Epan in energy-limited environments with constant rainfall may correspond to reductions in E (Roderick and Farquhar 2004), the influence of E on Epan in moisture-limited scenarios with variable precipitation (i.e., moisture availability for E) results in an inverse relationship between the two in which global dimming and E may increase together (Brutsaert 2006). Specifically, Brutsaert and Parlange’s (1998) interpretation of the decreasing Epan record in moisture-limited regions relies on the complementary relationship (CR; Bouchet 1963; Brutsaert 2006; Brutsaert and Parlange 1998; Morton 1965) hypothesis, that is, increases in E from progressively wetter land surfaces are accompanied by decreases in Epan until the two rates converge at a moisture-unlimited rate (Fig. 1; Bouchet 1963; Morton 1965). Despite its extensive validation over hundreds of watersheds worldwide (Bouchet 1963; Brutsaert and Stricker 1979; Crago and Crowley 2005; Diodato and Bellocchi 2007; Hobbins et al. 2001; Kahler and Brutsaert 2006; Liu et al. 2006; Morton 1983; Ozdogan and Salvucci 2004; Parlange and Katul 1992; Ramírez et al. 2005; Xu et al. 2006; Yang et al. 2007), a satisfactory explanation for the CR has remained elusive and controversial (Brutsaert 2005; Lhomme and Guilioni 2006; Pettijohn and Salvucci 2006). Here a new model is presented that illustrates how changes in regional evaporation influence pan evaporation through alteration of vapor and temperature fields above the pan. Results of the present study are in agreement with recent field experiments and indicate that the complementary relation between terrestrial evaporation and pan evaporation is robust and that the proportionality constant between positive changes in terrestrial evaporation and negative changes in pan evaporation is approximately 5, indicating that pans are sensitive instruments for measuring changes in evaporative flux.

Fig. 1.

The CR hypothesis between potential (pan) evaporation Epan and moisture-stressed terrestrial evaporation E.

Fig. 1.

The CR hypothesis between potential (pan) evaporation Epan and moisture-stressed terrestrial evaporation E.

Mathematically, Bouchet’s (1963) CR hypothesis can be expressed as (Brutsaert and Parlange 1998)

 
formula

In Eq. (1), E0 is the equilibrium evaporation for a wet surface (the convergence point on Fig. 1), and b is a constant of proportionality related to heat transfer efficiency between the pan and the environment (Bouchet 1963; Kahler and Brutsaert 2006; Morton 1965). To arrive at Eq. (1), Bouchet (1963) and Morton (1965), and others (Granger 1989; Szilagyi 2001) who have attempted to derive the CR from first principles, have made a series of strong assumptions regarding the role of energy balance and diffusion processes on evaporation from land and from evaporation pans. Various explanations have included assumptions that are mutually inconsistent, for example, according to Szilagyi (2001), the surface temperature of the evaporation pan remains constant as E decreases but the humidity of the air decreases, increasing the gradient over the pan and thus increasing Epan. Alternately, in Granger’s (1989) derivation of CR, complementarity results from an increase in surface temperature while screen height (i.e., 2 m) humidity remains constant, again increasing the humidity gradient and Epan, but for a mechanistically different reason. Lhomme and Guilioni (2006) suggest that both Granger’s (1989) and Szilagyi’s (2001) proofs of CR rely on physically questionable and unsupported assumptions.

Below it is found that by explicitly acknowledging—and modeling—the process as multidimensional, which has not to our knowledge been attempted, a priori assumptions about how the pan influences, and is influenced by, the air around it are avoided. It is then demonstrated that Eq. (1) follows directly from well-known and agreed upon physical principles of boundary layer meteorology.

2. Model physics

To solve for the coupled 2D steady-state atmospheric humidity and thermal field in the atmospheric dynamic sublayer (ASL; Fig. 2), the COMSOL Multiphysics 3.3 finite-element software package is used. The governing equation for water vapor being transported by advection and diffusion at steady state is

 
formula

In Eq. (2), ρ is the vapor concentration (kg m−3), D is the eddy diffusivity (m2 s−1) (described below), and u (m s−1) is the logarithmic mean wind profile characteristic of a neutral atmospheric sublayer (Brutsaert 2005):

 
formula

In Eq. (3), u* is the friction velocity (m s−1); κ is the von Kármán constant (=0.41); z is the height above the ground surface (m); zd is the zero-plane displacement height (m), which can be estimated (Dingman 2002) as 0.7zveg, where zveg is the surface vegetation height; and z0 is the surface roughness height, which can be estimated (Dingman 2002) as 0.1zveg. The corresponding diffusivity profile in the ASL is the following linear function of height (Dingman 2002):

 
formula

The above advection–diffusion equation for water vapor density ρ is then coupled with the following governing equation for steady-state heat transport:

 
formula

In Eq. (5), T is the temperature (K), ρa is the air density (kg m−3), Cp is the heat capacity (J kg−1 K−1), and k is the turbulent diffusivity for turbulent sensible heat transport (W m−1 K−1):

 
formula

The representation of diffusion in Eqs. (4) and (6) assumes that the eddy diffusivity for momentum transfer is the same for water vapor and sensible heat transport, a common assumption in evaporation studies for mean, neutral atmospheric stability conditions (Dingman 2002).

Fig. 2.

Conceptual diagram of the imposed 2D geometry and boundary conditions used in the COMSOL finite-element software for the steady-state vapor density ρ and temperature T solution within the ASL. The horizontal length scale Lx is taken as 50 m, and the depth of the ASL Lz is varied between 50 and 200 m. The potential (pan) evaporation Epan is situated along 24.5 m ≤ x ≤ 25.7 m and is represented by the ρ and T solution found iteratively in the initial equilibrium evaporation model in which Epan = E = E0 and SHpan = SH = SH0, where surface energy balance is held SH = RN − E.

Fig. 2.

Conceptual diagram of the imposed 2D geometry and boundary conditions used in the COMSOL finite-element software for the steady-state vapor density ρ and temperature T solution within the ASL. The horizontal length scale Lx is taken as 50 m, and the depth of the ASL Lz is varied between 50 and 200 m. The potential (pan) evaporation Epan is situated along 24.5 m ≤ x ≤ 25.7 m and is represented by the ρ and T solution found iteratively in the initial equilibrium evaporation model in which Epan = E = E0 and SHpan = SH = SH0, where surface energy balance is held SH = RN − E.

3. Numerical experiments

To investigate the validity of the CR [Eq. (1)] for moisture-limited environments and to explore the relation between the parameter b and characteristics of the ASL, a set of numerical experiments designed to model an evaporation pan in an increasingly dry landscape is carried out. The geometry of the model space is illustrated in Fig. 2. For all models, Lx = 50 m is assumed to provide adequate horizontal fetch; Lz is varied between 50 and 200 m to study the sensitivity of the CR to typical range in ASL depth. The evaporation pan boundary condition is defined along 24.5 m ≤ x ≤ 25.7 m, representative of a class-A evaporation pan diameter of ∼1.2 m.

First, the solution for the equilibrium evaporation condition E0 is obtained, that is, the convergence point of the CR (shown in Fig. 1) at which E = CpEpan = E0, where Cp is an empirical pan coefficient equal to the ratio of evaporation from a large body of water (as estimated by E0) to Epan. Whereas Cp can be calibrated to each site and typically has a value of 0.7 for elevated class-A evaporation pans, we assume a Cp equal to unity following Kahler and Brutsaert (2006). To solve for the equilibrium condition, an upper humidity boundary condition (Fig. 2) is imposed and the upper temperature boundary condition is adjusted iteratively (using COMSOL with MATLAB) until the lowest-level humidity and temperature, solved for under the imposed flux boundary conditions on evaporation E and sensible heat flux SH, yield a relative humidity of 1 [as occurs in the laminar sublayer of air (order of 1 mm) over a saturated surface (Penman 1948)]. Next, the effect of a drying landscape is simulated by incrementally decreasing the bottom evaporation flux boundary conditions while holding the humidity and temperature of the model pan (i.e., z = 0 m and 24.5 m ≤ x ≤ 25.7 m, representative of a class-A evaporation pan diameter of ∼1.2 m) and top boundary (i.e., z = Lz and 0 ≤ xLx) at the values found for the E = E0 case (see Fig. 2). For each incrementally lowered E, surface energy balance (RN − G = E + SH, where RN is the net available radiation at the land surface, and G is ground or pan heat flux, assumed to be negligible) was maintained by imposing equal incremental increases in the SH (see Fig. 2). As E was lowered (and SH increased), the modeled humidity decreased and the temperature increased, causing the system to dry out and the relative humidity near the surface to drop below 1. Pan evaporation Epan is then estimated for each moisture-case solution by integrating the flux of water vapor from the model pan (24.5 m ≤ x ≤ 25.7 m); SH from the pan was similarly found. This process was repeated for a range of realistic atmospheric conditions (e.g., 50 m ≤ Lz ≤ 200 m, characteristic of the typical range of ASL height) and surface roughness conditions (e.g., 0.2 m ≤ zveg ≤ 2 m).

Last, the 1D Penman (1948)Brutsaert (1982) apparent potential evaporation EPen-Brut was calculated at incremental reference heights (i.e., z = 0.5: 0.5: 5.0 m) upwind (x = 0 m) of the model Epan as

 
formula

In Eq. (7), Δ is the slope of the saturation vapor pressure curve (Pa K−1), γ is the psychrometric constant (Pa K−1), and EA,DP is the following “drying power of the air” term (W m−2; Brutsaert 1982):

 
formula

In Eq. (8), λυ is the latent heat of vaporization (J kg−1), qa* is the saturated specific humidity (kg kg−1) at reference height z, qa is the specific humidity (kg kg−1), and ga is the aerodynamic conductance for water vapor transport in a neutral atmosphere (m s−1):

 
formula

4. Discussion of results

Our results from the 2D model indicate that evaporation from a pan Epan increased proportionally to decreased regional evaporation E. These modeling results agree with Bouchet’s (1963) heuristic hypothesis and subsequent validation studies over hundreds of watersheds and climates worldwide on multiple time and length scales (Bouchet 1963; Brutsaert and Stricker 1979; Crago and Crowley 2005; Diodato and Bellocchi 2007; Hobbins et al. 2001; Kahler and Brutsaert 2006; Liu et al. 2006; Morton 1983; Ozdogan and Salvucci 2004; Parlange and Katul 1992; Ramírez et al. 2005; Xu et al. 2006; Yang et al. 2007).

The relation between Epan and E was found to be linear (as shown in Fig. 3) in which four suites of model results for different Lz and zveg combinations are presented. Each circle shown in Fig. 3 depicts the E/E0 and SH/SH0 boundary conditions and Epan/E0 and SHpan/SH0 model solutions, found by integrating the flux solution over the model evaporation pan shown in Fig. 2. This linearity between Epan and E (and SHpan and SH) implies that b in Eq. (1), which Brutsaert (2005) has interpreted as the fraction of liberated sensible heat that becomes available to increased pan evaporation, is constant—that is, the value of b is the same at the wet end (E close to E0) and the dry end (E close to 0) of the CR. Our model results agree with Bouchet’s (1963) hypothesized linearity, that is, ∂Epan/∂E = −b. The reason the relationship is linear is that the governing Eqs. (2) and (5) are linear, and thus the solutions ρ(x, z) and T(x, z) can be expressed as a superposition of the E = 0 case and the E = E0 case. Further, because the pan flux depends upon the gradient of ρ between the pan surface and the air above it, Epan can also be expressed as a superposition (weighted sum) of the two end-member solutions.

Fig. 3.

The inherent linearity and asymmetry of the CR over a range of realistic atmospheric turbulent surface sublayer heights and surface roughness geometries: (a) linear complementarity between Epan and E, where b ranges from 3.3 to 6.2 and (b) “mirrored” complementarity in sensible heat flux between pan (SHpan) and environmental SH.

Fig. 3.

The inherent linearity and asymmetry of the CR over a range of realistic atmospheric turbulent surface sublayer heights and surface roughness geometries: (a) linear complementarity between Epan and E, where b ranges from 3.3 to 6.2 and (b) “mirrored” complementarity in sensible heat flux between pan (SHpan) and environmental SH.

The sensitivity of Epan to E was found to be significantly greater than 1, varying between 3.3 and 6.2 and yielding a striking asymmetry in the complementary relationship. Moreover, the model results indicate that the asymmetry of CR (via b) is only sensitive to the simple geometric factors, zveg and Lz. The inherent linearity of the CR has perhaps been obscured in practice by use of various moisture indices (precipitation, point measurements in soil moisture at a single depth, etc.) for the x axis, which typically leads to a concave, nonlinear relation, as plotted in Fig. 1.

Because the dependence of Epan on E and also the dependence of SHpan on SH are linear, the energy balance of the pan is achieved with a constant pan temperature and humidity. Specifically, the increase in Epan caused by the drier air that results from the lowered regional evaporation is exactly equal to the decrease in SHpan caused by the hotter air that results from the increased regional sensible heat flux. This is not an assumption of the model but rather is a result found in the numerical solutions: the energy balance of the pan is only achieved if the pan temperature is held constant. This is in direct contradiction to the assumption made by Granger (1989) in his explanation of the CR but supports the deduction of Morton (1983) and Szilagyi (2001) that the temperature of the evaporimeter will not change in response to changes in E.

The CR was found to further simplify in that the sensitivity of Epan to E, that is, ∂Epan/∂E = −b, was independent of wind speed. Changes in wind speed increase advection and turbulent diffusion proportionally and thus their effects cancel. As the advection of dry air over a pan has been a primary focus in theoretical investigations of CR using one-dimensional equations (Brutsaert and Stricker 1979; Parlange and Katul 1992), the findings of our study add useful insight into the pan evaporation signature of advection of mass and energy from an increasingly arid environment. Furthermore, as some studies focusing on decreasing pan evaporation trends have found changes in wind speed to correspond with changes in Epan (Ozdogan and Salvucci 2004; Ozdogan et al. 2006; Rayner 2007; Roderick et al. 2007; Zhang et al. 2007), this preservation of CR with changing wind is a significant finding. In summary, for any set of model conditions, it was found that, while changes in wind speed can be a direct driver of E and Epan per Eq. (2), the overall wind speed does not affect the b parameter in Eq. (1); that is, the sensitivity of Epan to E is not altered by wind speed.

Simulations with other boundary conditions (not shown) indicate that the dependence of Epan on E does not depend on the upper boundary conditions (Ttop and ρtop), the overall RN, or the partitioning of RN into SH and E. Furthermore, the dependence of the CR on roughness height zveg and depth of the ASL Lz was found to be minimal, as illustrated in Fig. 3.

The results of this study illustrate that the dependence of Epan on E arises from multi-dimensional diffusion. The multidimensional nature of interaction between an evaporation pan and the environment is illustrated in the series of model solutions in Fig. 4. The four sets of panels highlight the ρ and T fields in a small area around the pan as the system progressively dries out, given Lz = 100 m and zveg = 0.2 m. Note that the gradients of temperature over the pan reverse under moderate drying conditions, leading to negative Bowen ratios and locally advected sensible heat (Table 1). In other words, the source of energy for pan evaporation is both net radiation and sensible heat, whereas for homogeneous landscapes the more typical relation is that both the sensible heat and latent heat fluxes cool the surface as it is heated by net radiation. Such two-dimensional dependence cannot be derived with one-dimensional diffusion equations, as has been attempted in the past (Granger 1989; Morton 1965; Szilagyi 2001) with Penman (1948)-based equations for open water evaporation.

Fig. 4.

The ρ and T solution surfaces around a pan, illustrating the importance of accounting for 2D mass and energy transfer involved in the CR, in contrast to previous 1D explanations. Regional moisture availability decreases from (a) to (d), as reflected in the decreasing imposed bottom boundary condition E from (a) E = 200 W m−2 to (d) 50 W m−2, as labeled.

Fig. 4.

The ρ and T solution surfaces around a pan, illustrating the importance of accounting for 2D mass and energy transfer involved in the CR, in contrast to previous 1D explanations. Regional moisture availability decreases from (a) to (d), as reflected in the decreasing imposed bottom boundary condition E from (a) E = 200 W m−2 to (d) 50 W m−2, as labeled.

Table 1.

Regional and pan fluxes during imposed drydown used in Fig. 4 (Lz = 100 m; zveg = 0.2 m). Here, RNpan = CpanEpan + SHpan.

Regional and pan fluxes during imposed drydown used in Fig. 4 (Lz = 100 m; zveg = 0.2 m). Here, RNpan = CpanEpan + SHpan.
Regional and pan fluxes during imposed drydown used in Fig. 4 (Lz = 100 m; zveg = 0.2 m). Here, RNpan = CpanEpan + SHpan.

Under 1D diffusion, evaporation rates can be expressed in terms of vertically integrated bulk difference equations, that is, E = −D(z)∂ρ/∂z is integrated to E = D+(ρ2ρ1)/(z2z1), where D+ is

 
formula

This step is performed in all common evaporation equations (Dingman 2002), including Penman, Dalton, Bowen, and Penman–Monteith. In a horizontally inhomogeneous vapor field such as occurs directly above a pan in a dry field, however, vertical integration of the diffusion equation must include the horizontal divergence of the vapor flux, which cannot be known a priori (i.e., without solving a model like that in Fig. 2). It is suggested that this is the fundamental reason for failure of these equations to represent pan evaporation and for previous problems with CR derivations. As a direct example, Fig. 5 demonstrates how the 1D EPen-Brut flux, calculated at various reference heights above the surface, upwind of Epan, significantly underestimates the Epan flux itself. This underestimation increases as z increases because of the logarithmic nature of the driving temperature and humidity gradients. Therefore, utilization of the standard z = 2 m in potential evaporation estimates for use in CR studies would drastically underestimate b, leading to poor indirect estimates of regional E using observed Epan per Eq. (1).

Fig. 5.

Underestimation of Epan using 1D EPen-Brut flux calculated at various reference heights zm above the land surface (Lz = 100 m; zveg = 0.2 m) caused by the failure of such 1D equations to represent the multidimensional advection-diffusion environment. As such, previous CR analyses using such 1D Penman-based equations for potential evaporation Epan have underestimated the true sensitivity of Epan to E.

Fig. 5.

Underestimation of Epan using 1D EPen-Brut flux calculated at various reference heights zm above the land surface (Lz = 100 m; zveg = 0.2 m) caused by the failure of such 1D equations to represent the multidimensional advection-diffusion environment. As such, previous CR analyses using such 1D Penman-based equations for potential evaporation Epan have underestimated the true sensitivity of Epan to E.

As an example of the difficulties associated with intrepreting observed pan evaporation trends, Rayner (2007) optimizes the Thom et al. (1981) class-A pan evaporation model for pan evaporation in Australia from 1975 to 2004 and finds only moderate correlations (r2 ≈ 0.45–0.65) between daily observed and measured pan evaporation anomalies (implying even lower r2 in the correlation between observed and modeled Epan fluxes with systematic bias). Part of the difficulty in such pan evaporation correlations may be attributed to the fact that, while Thom et al. (1981) intend driving parameters such as humidity to be taken “at the pan rim height,” such measurements are typically only available at the standard 2-m reference height. As seen in Fig. 4, the conditions at a height of even 1 m would differ significantly from those at the immediately pan edge and would depend on the prevailing wind direction. Ultimately, these mixing processes determine the gradients of humidity and temperature directly over the pan and thus the pan evaporation rate. These processes are fundamentally—at a minimum—two (if not three) dimensional. Therefore, even though 1D equations can be calibrated to explain a significant amount of pan evaporation variance, as conducted by Roderick et al. (2007) for monthly pan evaporation rates across a wide variety of climates using a suitably modified combination equation, such empirical equations cannot be used to successfully understand the multidimensional physics of an evaporation pan in a drying landscape. The purpose of this study is not to dismiss such valuable studies but rather to resolve the poorly understood CR hypothesis and ultimately understand the implications that twentieth-century pan evaporation trends have on actual E.

As an example of the problems that neglecting horizontal divergence of water vapor causes in explanations of CR, Szilagyi’s (2001) derivation concludes that CR is perfectly symmetrical (b = 1), whereas our full model (e.g., see Fig. 3) demonstrates it to be in the range of 3–6. Because he bases his derivation on the one-dimensional diffusion equation, he is forced to make an assumption regarding blending geometry (i.e., the vertical, 1D vapor pressure profile gradients for the region and the juxtaposed 1D, vertical vapor pressure profile between the hypothetical evaporimeter and an assumed blending height), which he expresses through his choice of a “lower reference level at a height where originally the vapor pressure is half of that at the evaporimeter level.” This choice leads him to find b = 1. Our results indicate that this choice is not arbitrary but rather is a product of the solutions of Eqs. (2) and (5). That this reference level should not be arbitrary is also pointed out by Lhomme and Guilioni (2006).

5. Comparison with field data

As demonstrated in Figs. 6a and 6b, the asymmetry of CR seen in our model results are very comparable to the Kahler and Brutsaert (2006) presentation of observed E and Epan relationships measured at the Konza Prairie Research Natural Area and the Little Washita River basin. The E data used by Kahler and Brutsaert (2006) for the Konza Prairie were measured at Station 40 (1987) and Station 944 (1989) near Manhattan, Kansas, during the First International Satellite Land Surface Climatology Project (ISLSCP) Field Experiment (FIFE) during the summers of 1987–89 (Sellers et al. 1992), available from Oak Ridge Nation Laboratory (http://daac.ornl.gov/FIFE/FIFE_Home.html). The Epan data were measured by the U.S. Corps of Engineers at a mowed, nonirrigated field at the Turkey Creek Reservoir, some 16.7 km away at 337° bearing from Station 40 (Kahler and Brutsaert 2006), available from the National Climatic Data Center (http://lwf.ncdc.noaa.gov/oa/ncdc.html). Next, Kahler and Brutsaert (2006) used the 1993–2003 (May–September) E data for the Little Washita basin, measured at the Southern Great Plains (SGP) Station 4 near Cement, Oklahoma, operated by the Atmospheric Radiation Measurement (ARM) Program of the U.S. Department of Energy Office of Biological and Environmental Research Environmental Sciences Division (Kustas et al. 1996; Stannard et al. 1993). The Epan was measured in a nonirrigated, mowed, Bermuda grass lawn at the Oklahoma State University Research Station, located approximately 18.3 km from the ARM Station 4 at a 52° bearing (Kahler and Brutsaert 2006).

Fig. 6.

Comparison of CR results found using the 2D model in this study with the Kahler and Brutsaert (2006) analysis of normalized data from the (a) Konza Prairie and (b) Little Washita basin. Fluxes are normalized as Epan+ = (CpEpan)/E0 (circles), E+ = E/E0 (triangles), and EMI = E/(CpEpan). This figure is reproduced/modified from Figs. 4 and 5 of Kahler and Brutsaert (2006).

Fig. 6.

Comparison of CR results found using the 2D model in this study with the Kahler and Brutsaert (2006) analysis of normalized data from the (a) Konza Prairie and (b) Little Washita basin. Fluxes are normalized as Epan+ = (CpEpan)/E0 (circles), E+ = E/E0 (triangles), and EMI = E/(CpEpan). This figure is reproduced/modified from Figs. 4 and 5 of Kahler and Brutsaert (2006).

In Fig. 6, the Kahler and Brutsaert (2006) dimensionless x axis moisture-index variable is used, that is, EMI = E/(CpEpan), where Cp is the proportionality constant and is assumed to equal 1; similarly, the fluxes are normalized on their y axis, where Epan becomes Epan+ = (CpEpan)/E0 and E is normalized as E+ = E/E0. As shown, our model results over a range of physical conditions agree well with their analysis of measured pan data and measured land surface evaporation. Kahler and Brutsaert (2006) found a b = 4.33 to best fit the Konza Prairie data and a b = 6.88 to best fit the Little Washita basin. A similar range of b values, between 3.3 and 6.2, was found to result from a realistic range of model conditions.

In Fig. 6, the model with Lz = 200 m and zveg = 0.2 m best fits (out of the four models presented) the Konza Prairie data while the model with Lz = 200 m and zveg = 2.0 m best fits the Little Washita data. The model zveg = 0.2 m for the Konza Prairie corresponds well to the range of reported canopy heights (e.g., average zveg = 0.28 and 0.42 m for sites 26 and 16, respectively) for the grassland sites (Stewart and Verma 1992) while exceeding the local mowed-grass canopy height surrounding the evaporation pan itself; the model Lz = 200 m is approximately 100 m higher than the log mean ASL height reported by Sugita and Brutsaert (1990) from radiosonde profiles. Although the model zveg = 2.0 m for the Little Washita site exceeds the range of surface roughness heights for bare soil, corn, wheat, alfalfa, and rangeland (Shiebe et al. 1993), Jacobs and Brutsaert (1998) find a regional momentum roughness height (z0m) equal to 0.2 m from the intensive radiosounding program described in Schiebe et al. (1993). Given that z0m is typically 0.1zveg, the best-fit model zveg = 2.0 m for this site is not unreasonable when considering topographic effects of regional surface roughness. Furthermore, atmospheric profiles of potential temperature and specific humidity taken on 10 and 12 June 1992 during an intensive radiosounding program by W. Brutsaert and N. Dias, as described in Schiebe et al. (1993), indicate that Lz ranges diurnally from approximately 50 to 175 m. While it is encouraging that our simplified 2D model of surface and atmospheric conditions at the two sites yields comparable pan evaporation trends, care should be taken because of the model assumptions made. Specifically, considering that the model upper temperature and humidity boundary conditions are set to force a driving gradient (and are not calibrated to actual temperature and humidity profiles at the sites), it is not logical to infer too much from direct zveg and Lz comparisons, apart from first-order sensitivity analyses. Furthermore, the 2D model used in this study assumes that horizontal heterogeneities (in roughness or other surface characteristics) that are likely to occur over scales of kilometers, as would be found between E and Epan measurements used by Kahler and Brutsaert (2006), are insignificant. The good fit between modeled and measured Epan shown in Fig. 6 suggests that the results of this study may be applicable over longer horizontal spatial scales, despite the Lx = 50 m assumption. In other words, whereas b may change slightly because of horizontal differences in Lz and zveg, the general complementarity between Epan and E holds regardless of boundary conditions (and horizontal variations thereof), given variable moisture availability.

Note that, while Kahler and Brutsaert (2006) found b ∼ 5 at the daily scale, others (Brutsaert and Stricker 1979; Morton 1983; Parlange and Katul 1992; Yang et al. 2006) have found a b closer to 1, as originally suggested by Bouchet (1963) and Morton (1965). In addressing this issue, Kahler and Brutsaert (2006) note that many previous studies average data over long periods (from weeks to months) and that using such data in turbulence calculations can lead to large biases. Furthermore, a review of these works indicates that in most cases Epan was not measured but rather was modeled using Penman (1948)-based equations, which, as demonstrated above, do not represent the true multidimensional, coupled advection–diffusion physics of an evaporation pan in a drying landscape and, therefore, may not provide the best theoretical cornerstone for inferences into E trends.

While the relationship between terrestrial and pan fluxes in the energy-limited scenario is not explored in this study, the model representation thereof would be similar to the equilibrium solution described previously, but for variable RN conditions. As such, the equal regional and pan fluxes would simply shift vertically on Fig. 1 at the moisture-unlimited convergence point.

6. Summary and conclusions

In summary, the following conclusions have been reached, under the assumptions made herein to model an evaporation pan in a drying environment:

  1. The relation between E and Epan is linear. This implies that b in Eq. (1), which Bouchet (1963) has interpreted as the fraction of liberated sensible heat that becomes available to increased pan evaporation, is constant (i.e., does not change with E).

  2. The sensitivity of Epan to E, that is, ∂Epan/∂E = −b, does not depend upon wind speed. Changes in wind speed increase advection and turbulent diffusion [through Eqs. (3) and (4)] proportionally and influence E and Epan identically, and thus the wind effect cancels out.

  3. The dependence of Epan on E does not depend on the upper boundary conditions (Ttop and ρtop) or the overall available RN. Dependence on a realistic range of roughness height zveg and depth of the ASL Lz is minimal, as seen in Fig. 3.

  4. The sensitivity of Epan to E is significantly greater than 1. It varies in a moderately small range around 5, yielding a striking asymmetry in the complementary relationship, as observed by Kahler and Brutsaert (2006).

  5. The dependence of Epan on E arises from two-dimensional diffusion. The dependence cannot be derived with one-dimensional diffusion equations, as has been attempted in the past (Granger 1989; Morton 1965; Szilagyi 2001).

  6. Because the dependence of Epan on E and also the dependence of SHpan on SH are linear, the energy balance of the pan is achieved with a constant pan temperature and humidity. This is in direct contradiction to the assumption made by Granger (1989) in his explanation of CR.

Our physical model of the CR in a simple 2D advective–diffusive environment demonstrates how complementarity between E and Epan in a moisture-limited environment results from the fact that the water vapor and heat leaving (or entering) the pan diffuse both laterally and vertically into the surrounding environment. The physical representation of evaporation from a pan illustrated in Fig. 2 contains the minimal set of assumptions required to explain the widely observed CR phenomenon. Improvements to the model for future studies might include an investigation of how complementarity is influenced by 1) the addition of a third dimension, 2) non-steady-state dynamics, 3) stable and unstable atmospheric conditions, and 4) different pan evaporation designs and thus different Cp magnitudes.

Consistent with recent observations (Kahler and Brutsaert 2006), our results imply that, given a moisture-limited environment with variable rainfall (and thus moisture availability), a unit increase in land surface evaporation could explain an approximately fivefold decrease in Epan. Our analysis improves on the theoretical explanation of the CR, thus lending support to the interpretation (Brutsaert 2006; Brutsaert and Parlange 1998) that measured decreases in pan evaporation over the twentieth century in moisture-limited environments could be the result of increases in terrestrial evaporation corresponding to well-documented increases in terrestrial rainfall in many regions of the world for the last 50–100 years (Folland et al. 2001).

Acknowledgments

This work was supported by NASA Grants NAG5-11338 and 11695 and NSF Grant EAR0233643.

REFERENCES

REFERENCES
Bouchet
,
R. J.
,
1963
:
Évapotranspiration réelle et potentielle signification climatique.
Evaporation, IAHS Publication 62, IAHS, 134–142
.
Brutsaert
,
W.
,
1982
:
Evaporation into the Atmosphere: Theory, History, and Applications.
Kluwer Academic, 299 pp
.
Brutsaert
,
W.
,
2005
:
Hydrology: An Introduction.
1st ed. Cambridge University Press, 605 pp
.
Brutsaert
,
W.
,
2006
:
Indications of increasing land surface evaporation during the second half of the 20th century.
Geophys. Res. Lett.
,
33
,
L20403
.
doi:10.1029/2006GL027532
.
Brutsaert
,
W.
, and
H.
Stricker
,
1979
:
Advection-aridity approach to estimate actual regional evapotranspiration.
Water Resour. Res.
,
15
,
443
450
.
Brutsaert
,
W.
, and
M. B.
Parlange
,
1998
:
Hydrologic cycle explains the evaporation paradox.
Nature
,
396
,
30
.
Crago
,
R.
, and
R.
Crowley
,
2005
:
Complementary relationships for near-instantaneous evaporation.
J. Hydrol.
,
300
,
199
211
.
Dingman
,
S. L.
,
2002
:
Physical Hydrology.
2nd ed. Prentice Hall, 646 pp
.
Diodato
,
N.
, and
G.
Bellocchi
,
2007
:
Modeling reference evapotranspiration over complex terrains from minimum climatological data.
Water Resour. Res.
,
43
,
W05444
.
doi:10.1029/2006WR005405
.
Folland
,
C. K.
, and
Coauthors
,
2001
:
Observed climate variability and change.
Climate Change 2001: The Scientific Basis, J. T. Houghton et al., Eds., Cambridge University Press, 99–181
.
Granger
,
R. J.
,
1989
:
A complementary relationship approach for evaporation from nonsaturated surfaces.
J. Hydrol.
,
111
,
31
38
.
Hobbins
,
M. T.
,
J. A.
Ramirez
, and
T. C.
Brown
,
2001
:
The complementary relationship in estimation of regional evapotranspiration: An enhanced advection-aridity model.
Water Resour. Res.
,
37
,
1389
1403
.
Jacobs
,
J. M.
, and
W.
Brutsaert
,
1998
:
Momentum roughness and view-angle dependent heat roughness at the Southern Great Plains Atmospheric Research Measurement test-site.
J. Hydrol.
,
211
,
61
68
.
Kahler
,
D. M.
, and
W.
Brutsaert
,
2006
:
Complementary relationship between daily evaporation in the environment and pan evaporation.
Water Resour. Res.
,
42
,
W05413
.
doi:10.1029/2005WR004541
.
Kustas
,
W. P.
,
D. I.
Stannard
, and
K. J.
Allwine
,
1996
:
Variability in surface energy flux partitioning during Washita ’92: Resulting effects on Penman–Monteith and Priestley–Taylor parameters.
Agric. For. Meteor.
,
82
,
171
193
.
Lhomme
,
J. P.
, and
L.
Guilioni
,
2006
:
Comments on some articles about the complementary relationship.
J. Hydrol.
,
323
,
1
3
.
Liu
,
S. M.
,
R.
Sun
,
Z. P.
Sun
,
X. O.
Li
, and
C. M.
Liu
,
2006
:
Evaluation of three complementary relationship approaches for evapotranspiration over the Yellow River basin.
Hydrol. Processes
,
20
,
2347
2361
.
Morton
,
F. I.
,
1965
:
Potential evaporation and river basin evaporation.
J. Hydraul. Div. Amer. Soc. Civ. Eng.
,
91
,
67
97
.
Morton
,
F. I.
,
1983
:
Operational estimates of areal evapotranspiration and their significance to the science and practice of hydrology.
J. Hydrol.
,
66
,
1
76
.
Ozdogan
,
M.
, and
G. D.
Salvucci
,
2004
:
Irrigation-induced changes in potential evapotranspiration in southeastern Turkey: Test and application of Bouchet’s complementary hypothesis.
Water Resour. Res.
,
40
,
W04301
.
doi:10.1029/2003WR002822
.
Ozdogan
,
M.
,
G. D.
Salvucci
, and
B. T.
Anderson
,
2006
:
Examination of the Bouchet–Morton complementary relationship using a mesoscale climate model and observations under a progressive irrigation scenario.
J. Hydrometeor.
,
7
,
235
251
.
Parlange
,
M. B.
, and
G. G.
Katul
,
1992
:
An advection–aridity evaporation model.
Water Resour. Res.
,
28
,
127
132
.
Penman
,
H. L.
,
1948
:
Natural evaporation from open water, bare soil and grass.
Philos. Trans. Roy. Soc. London
,
A193
,
120
145
.
Peterson
,
T. C.
,
V. S.
Golubev
, and
P. Ya
Groisman
,
1995
:
Evaporation losing its strength.
Nature
,
377
,
687
688
.
Pettijohn
,
J. C.
, and
G. D.
Salvucci
,
2006
:
Impact of an unstressed canopy conductance on the Bouchet–Morton complementary relationship.
Water Resour. Res.
,
42
,
W09418
.
doi:10.1029/2005WR004385
.
Ramírez
,
J. A.
,
M. T.
Hobbins
, and
T. C.
Brown
,
2005
:
Observational evidence of the complementary relationship in regional evaporation lends strong support for Bouchet’s hypothesis.
Geophys. Res. Lett.
,
32
,
L15401
.
doi:10.1029/2005GL023549
.
Rayner
,
D. P.
,
2007
:
Wind run changes: The dominant factor affecting pan evaporation trends in Australia.
J. Climate
,
20
,
3379
3394
.
Roderick
,
M. L.
, and
G. D.
Farquhar
,
2004
:
Changes in Australian pan evaporation from 1970 to 2002.
Int. J. Climatol.
,
24
,
1077
1090
.
Roderick
,
M. L.
,
L. D.
Rotstayn
,
G. D.
Farquhar
, and
M. T.
Hobbins
,
2007
:
On the attribution of changing pan evaporation.
Geophys. Res. Lett.
,
34
,
L17403
.
doi:10.1029/2007GL031166
.
Schiebe
,
F. R.
,
J. W.
Naney
,
P. B.
Allen
, and
T. J.
Jackson
,
1993
:
Little Washita watershed.
Hydrology Data Report: Washita ’92, U.S. Department of Agriculture Agriculture Research Services Rep. NAWQL 93-1, II-1–II-10. [Available online at http://hydrolab.arsusda.gov/washita92/datarpt.htm.]
.
Sellers
,
P. J.
,
F. G.
Hall
,
G.
Asrar
,
D. E.
Strebel
, and
R. E.
Murphy
,
1992
:
An overview of the First International Satellite Land Surface Climatology Project (ISLSCP) Field Experiment (FIFE).
J. Geophys. Res.
,
97
,
18345
18371
.
Stannard
,
D. I.
,
W. P.
Kustas
,
K. J.
Allwine
, and
D. E.
Anderson
,
1993
:
Micrometeorological data collection.
Hydrology Data Report: Washita ’92, U.S. Department of Agriculture Agriculture Research Service Rep. NAWQL 93-1, V-1–V-19. [Available online at http://hydrolab.arsusda.gov/washita92/datarpt.htm.]
.
Stewart
,
J. B.
, and
S. B.
Verma
,
1992
:
Comparison of surface fluxes and conductances at two contrasting sites within the FIFE area.
J. Geophys. Res.
,
97
,
18623
18628
.
Sugita
,
M.
, and
W.
Brutsaert
,
1990
:
Wind velocity measurements in the neutral boundary layer above hilly prairie.
J. Geophys. Res.
,
95
,
7617
7624
.
Szilagyi
,
J.
,
2001
:
On Bouchet’s complementary hypothesis.
J. Hydrol.
,
246
,
155
158
.
Thom
,
A. S.
,
J. L.
Thony
, and
M.
Vauclin
,
1981
:
On the proper employment of evaporation pans and atmometers in estimating potential transpiration.
Quart. J. Roy. Meteor. Soc.
,
107
,
711
736
.
Xu
,
C. Y.
,
L. B.
Gong
,
T.
Jiang
,
D. L.
Chen
, and
V. P.
Singh
,
2006
:
Analysis of spatial distribution and temporal trend of reference evapotranspiration and pan evaporation in Changjiang (Yangtze River) catchment.
J. Hydrol.
,
327
,
81
93
.
Yang
,
D. W.
,
F. B.
Sun
,
Z. T.
Liu
,
Z. T.
Cong
, and
Z. D.
Lei
,
2006
:
Interpreting the complementary relationship in non-humid environments based on the Budyko and Penman hypotheses.
Geophys. Res. Lett.
,
33
,
L18402
.
doi:10.1029/2006GL027657
.
Yang
,
D. W.
,
F. B.
Sun
,
Z. Y.
Liu
,
Z. T.
Cong
,
G. H.
Ni
, and
Z. D.
Lei
,
2007
:
Analyzing spatial and temporal variability of annual water-energy balance in nonhumid regions of China using the Budyko hypothesis.
Water Resour. Res.
,
43
,
W04426
.
doi:10.1029/2006WR005224
.
Zhang
,
Y.
,
C.
Liu
,
Y.
Tang
, and
Y.
Yang
,
2007
:
Trends in pan evaporation and reference and actual evapotranspiration across the Tibetan Plateau.
J. Geophys. Res.
,
112
,
D12110
.
doi:10.1029/2006JD008161
.

Footnotes

Corresponding author address: J. Cory Pettijohn, Department of Forest Ecosystems and Society, Oregon State University, Corvallis, OR 97331. Email: cory.pettijohn@oregonstate.edu