Forestry, Faculty of
Non UBC
DSpace
Hydrology 4 (1): 5 (2017)
Hember, Robbie A.
Coops, Nicholas C.
Spittlehouse, David L.
Given the widespread ecological implications that would accompany any significant change in evaporative demand of the atmosphere, this study investigated spatial and temporal variation in several accepted expressions of potential evaporation (PE). The study focussed on forest regions of North America, with 1 km-resolution spatial coverage and a monthly time step, from 1951–2014. We considered Penman’s model (EPen), the Priestley–Taylor model (EPT), ‘reference’ rates based on the Penman–Monteith model for grasslands (ERG), and reference rates for forests that are moderately coupled (ERFu) and well coupled (ERFc) to the atmosphere. To give context to the models, we also considered a statistical fit (EPanFit) to measurements of pan evaporation (EPan). We documented how each model compared with EPan, differences in attribution of variance in PE to specific driving factors, mean spatial patterns, and time trends from 1951–2014. The models did not agree strongly on the sensitivity to underlying drivers, zonal variation of PE, or on the magnitude of trends from 1951–2014. Sensitivity to vapour pressure deficit (Da) differed among models, being absent from EPT and strongest in ERFc. Time trends in reference rates derived from the Penman–Monteith equation were highly sensitive to how aerodynamic conductance was set. To the extent that EPanFit accurately reflects the sensitivity of PE to Da over land surfaces, future trends in PE based on the Priestley–Taylor model may underestimate increasing evaporative demand, while reference rates for forests, that assume strong canopy-atmosphere coupling in the Penman–Monteith model, may overestimate increasing evaporative demand. The resulting historical database, covering the spectrum of different models of PE applied in modern studies, can serve to further investigate biosphere-hydroclimate relationships across North America.
hydrologyArticleSpatial and Temporal Variability of PotentialEvaporation across North American ForestsRobbie A. Hember 1,2,*, Nicholas C. Coops 1 and David L. Spittlehouse 31 Faculty of Forestry, University of British Columbia, 2424 Main Mall, Vancouver, BC V6T 1Z4, Canada;nicholas.coops@ubc.ca2 Pacific Forestry Centre, Canadian Forest Service, Natural Resources Canada, 506 West Burnside Road,Victoria, BC V8Z 1M5, Canada3 Competitiveness and Innovation Branch, Ministry of Forests, Lands and Natural Resources Operations,Victoria, BC V8W 9C2, Canada; dave.spittlehouse@gov.bc.ca* Correspondence: robert.hember@canada.ca; Tel.: +1-250-298-2339Academic Editor: Matthew McCabeReceived: 18 October 2016; Accepted: 4 January 2017; Published: 11 January 2017Abstract: Given the widespread ecological implications that would accompany any significant changein evaporative demand of the atmosphere, this study investigated spatial and temporal variationin several accepted expressions of potential evaporation (PE). The study focussed on forest regionsof North America, with 1 km-resolution spatial coverage and a monthly time step, from 1951–2014.We considered Penman’s model (EPen), the Priestley–Taylor model (EPT), ‘reference’ rates based onthe Penman–Monteith model for grasslands (ERG), and reference rates for forests that are moderatelycoupled (ERFu) and well coupled (ERFc) to the atmosphere. To give context to the models, we alsoconsidered a statistical fit (EPanFit) to measurements of pan evaporation (EPan). We documentedhow each model compared with EPan, differences in attribution of variance in PE to specific drivingfactors, mean spatial patterns, and time trends from 1951–2014. The models did not agree stronglyon the sensitivity to underlying drivers, zonal variation of PE, or on the magnitude of trends from1951–2014. Sensitivity to vapour pressure deficit (Da) differed among models, being absent from EPTand strongest in ERFc. Time trends in reference rates derived from the Penman–Monteith equationwere highly sensitive to how aerodynamic conductance was set. To the extent that EPanFit accuratelyreflects the sensitivity of PE to Da over land surfaces, future trends in PE based on the Priestley–Taylormodel may underestimate increasing evaporative demand, while reference rates for forests, thatassume strong canopy-atmosphere coupling in the Penman–Monteith model, may overestimateincreasing evaporative demand. The resulting historical database, covering the spectrum of differentmodels of PE applied in modern studies, can serve to further investigate biosphere-hydroclimaterelationships across North America.Keywords: potential evaporation; evaporative demand; pan evaporation; climate trends1. IntroductionThe evaporating power of the atmosphere, or simply evaporative demand, exerts a primarycontrol on the distribution, abundance, and diversity of life on land [1–3]. On one hand,the atmospheric demand for water allows plants to draw water above the surface without exertingvast amounts of energy. On the other, fluctuations in evaporative demand simultaneously pose a riskto plants if soil water supply and transport efficiency cannot meet peak evaporative demand [4–6].These compensating factors appear to put forest biomes at risk under anticipated climate change [7–11].While declining precipitation must impose major constraints on productivity of trees across theforest-grassland ecotone, there is a growing appreciation that evaporative demand also stronglyHydrology 2017, 4, 5; doi:10.3390/hydrology4010005 www.mdpi.com/journal/hydrologyHydrology 2017, 4, 5 2 of 24constrains productivity of humid forest biomes [12–14]. For these reasons, evaporative demand is animportant flux in the surface water balance and a primary determinant of land climates [1,3,15,16].Potential evaporation (PE) is a common measure of evaporative demand and is defined as therate of evaporation that would occur from a surface with unlimited water supply and no resistance totransfer of water. It is, therefore, determined by the energy available to change the phase of water fromliquid to gas [17–20]. Although PE has a conceptual definition, it is not easy to measure directly aboveland surfaces [21]. Inferences can be made from measurements of evaporation from lakes, lysimeters,evaporation pans, or models. For example, Class A pan evaporation was measured continuously atmany meteorological stations across Canada and the U.S. and has been used to assess climatologicalvariability of PE [22,23].Despite extensive investigation of evaporation from pans, lakes, and lysimeters, empirical studieshave not resulted in a universal approach to predict PE. Some of the earliest models were basedon temperature and were developed with the limited availability of radiation, humidity, and windmeasurements in mind [16,24,25]. Other models were simultaneously developed to predict PE basedon a broader number of driving variables, including net radiation (Rn), temperature (Ta), vapourpressure deficit (Da), and wind speed (U) [26,27]. Later, Penman’s model was modified to produce thePenman–Monteith equation, which accounts for varying surface conductance to water transfer [28].While the Penman–Monteith model incorporates an aerodynamic component, other models weredeveloped that omit an aerodynamic component. For example, it was later argued that at regionalscales PE could be approximated by multiplying the equilibrium rate of evaporation by a constantvalue, α = 1.26 [29]. Although empirical studies suggest that the overall magnitude of PE using thisPriestley–Taylor coefficient of α = 1.26 will bring estimates within closer proximity to models thatinclude an aerodynamic component, spatial and temporal variability in estimates will only reflectvariation in Ta and Rn. In forests, α must be modified [30,31], yet it is still a useful concept if it isrecognized that there is no universal value of α. Lastly, other studies express PE based solely on afunction of Rn [32], or solely based on Da [33].The above modelling approaches all continue to be used broadly in hydrological and ecologicalresearch. For example, temperature-based models are applied in relatively basic surface waterbalance models [34]. Use of Penman’s model has also been applied in ecological studies [2].The Penman–Monteith model is widely applied [35–38]. Other studies represent potentialevapotranspiraiton based on the Priestley–Taylor model [30,39], or close equivalents based on Rn andTa [25,40], or Rn alone [32], or Da alone [41,42]. Despite its importance in modelling the surface waterbalance and subsequent application in ecological research, differences in the variation of evaporativedemand resulting from these approaches, including long-term trends at the regional scale, are notwell documented.Historical and projected changes in PE may arise from natural climate variability, and changingatmospheric greenhouse gas and aerosol concentrations. However, the response of PE to regionalclimate change involves scale- and surface-dependent relationships with multiple meteorologicaldrivers that are not always directly measured. As such, continuous high temporal and spatialresolution maps of PE and projections of future change required to study many ecological processesremain relatively scarce. For regions with long-term continuous measurements of pan evaporation,some studies report a decrease in pan evaporation [43–45]. Other studies resort to use of models tounderstand trends. For example, Penman’s model suggested a decrease in historical PE across someregions of India [43], but also projected long-term future increases of PE. Across Australia, trends in PEinferred from a temperature-based model and Class-A pan evaporation did not agree in direction [46].Similarly, future projections of PE based on several models, including the Priestley–Taylor model andthe Penman–Monteith model, concluded that differences among models partially led to differences inthe direction of future trends in aridity index [47]. Trends in PE based on the Penman–Monteithequation generally indicate increasing trends in evaporative demand [7,48–50]. Positive futuretrends in PE are largely attributed to increasing saturation vapour pressure that accompanies globalHydrology 2017, 4, 5 3 of 24warming [48]. Comparison of historical trends in drought using Penman’s model with previousstudies that applied Thornthwaite’s temperature-based model of PE emphasized that the latter mightexaggerate increasing trends in aridity [51]. Conversely, time trends in global land aridity of CoupledModel Intercomparison Project Phase 5 (CMIP5) model simulations based on estimates of PE inferredfrom just net radiation were commonly not significant [32].There is considerable demand for a spatial and temporal index of change in evaporative demandto understand how forests may have responded over time to changing atmospheric conditions [52,53].Consequently, we investigated the spatial and temporal variability of PE across North Americausing six different model-based estimates of PE. To put the magnitude and underlying mechanics ofmodel-based estimates of PE into context, we compared predictions against measurements of panevaporation (EPan) for 28 sites across Canada and the conterminous U.S.. One of the six estimatesof PE was defined by a statistical model fit to EPan measurements, while the other five estimatesreflected approaches that are widely used in applied research, but that require no calibration. First,we evaluated the agreement between models and EPan. Second, we evaluated the variance in PE thatis attributed to each meteorological driver in each of the models. Lastly, to understand spatial andtemporal variation in the model estimates, we developed spatially-continuous monthly mean values ofthe driving variables covering North America from 1951–2014. Maps were then generated to documentlong-term averages (‘climatologies’) of PE, and historical time trends in PE from 1951–2014.2. Data and Methods2.1. Pan Evaporation MeasurementsDaily total pan evaporation (EPan) measurements from Class A pans were collected across NorthAmerica. Data collected at stations in Canada were requested and received from Environment Canadaduring March 2014. Data collected at stations in the U.S. were downloaded from the National Oceanicand Atmospheric Administration during March 2016. Any stations with available daily measurementsbetween 1950 and 2014 were considered.In Canada, large periods of the pan evaporation measurements were accompanied bysimultaneous measurements of pan water temperature, air temperature (Ta), and wind speed (U).Separate databases of daily mean downward global shortwave radiation (Rsgd) and dew-pointtemperature were also available. Daily mean vapour pressure deficit of the air (Da) was calculated fromair temperature and dew-point temperature. From the initial databases, ten stations in Canada hadmultiple years with overlapping measurements of EPan and meteorological variables (all of Rsgd, Ta, Da,and U). For all U.S. stations, daily mean values of Rsgd were downloaded from two National RenewableEnergy Laboratory (NREL) databases [54]. All other meteorological variables were downloaded fromthe Global Summary of the Day (GSOD) database [55]. Relatively few U.S. stations simultaneouslymeasured pan evaporation and meteorological drivers. By manually cross-referencing betweendatabases, we identified eighteen stations that collected pan evaporation measurements and were alsowithin close proximity to stations with the necessary meteorological measurements (a distance of lessthan 50 km in most cases). Daily values of Rsgd, Ta, Da, and U for each station with pan evaporationmeasurements were then estimated from inverse-distance weighting interpolation of measurementsfrom the ten nearest stations. As a quality control measure, multivariate regression was used to relatepan evaporation against meteorological variables at each candidate station. Despite reliance on off-sitemeteorological variables in the U.S., we did not observe a significant difference in the skill or behaviourof the models between Canada and the U.S. (not shown). Combining the stations from Canada andU.S. produced a total sample of 28 stations, each consisting of multiple years of measurements foranalysis (Table 1). Stations were relatively evenly distributed across North America, with the exceptionof Mexico (Figure 1). Although the emphasis was on documenting variability of PE for forest biomesof North America, we chose to analyze all 28 available stations.Hydrology 2017, 4, 5 4 of 24Table 1. Station information: Identifiers (ID) from Environment Canada or NOAA, geographic coordinates, mean annual temperature (MAT) and mean annualprecipitation (MAP) for the 1971–2000 base period, number of days with measurements of pan evaporation (EPan) used in the analysis (n), and time period ofcontaining EPan measurements.Station ID Station Name Lat. (◦N) Long. (◦W) Elev. (m a.s.l.) MAT (◦C) MAP (mm·year−1) N (d) Period2101300 Whitehorse A 60.71 135.07 706 −0.6 274 2787 1974–19962202800 Norman Wells 65.28 126.80 73 −4.6 294 2573 1964–19842403500 Resolute Cars 74.72 94.99 66 −15.2 127 2231 1962–2004301222F Edmonton Stony Plain 53.55 114.11 766 2.7 530 2654 1988–20045023222 Winnipeg Richarson Int’l A 49.92 97.23 239 3.2 545 5422 1962–19965060600 Churchill A 58.74 94.06 29 −6.9 385 3672 1963–20007025250 Pierre Elliott Trudau Intl A 45.47 73.75 36 6.3 1059 1173 1987–19947093GJ3 La Grande IV A 53.76 73.68 306 −5 732 955 1987–19957095480 Nitchequon 53.20 70.90 536 −7.6 845 2261 1966–19857117825 Schefferville A 54.80 66.82 522 −7.7 772 818 1962–1969040232 Antioch Pump Plant 3 37.98 121.73 18 15.7 346 7966 1955–1978050834 Bonny Dam 2 NE 39.63 102.18 1113 10.8 466 7601 1951–2006090181 Allatoona Dam 2 34.17 84.73 299 16 1347 10,578 1952–2006100448 Arrowrock Dam 43.62 115.92 988 8.7 569 6630 1951–1998151080 Buckhorn Lake 37.35 83.38 287 13 1199 6049 1961–1994200089 Alberta Ford For Cen 46.65 88.48 400 3.2 867 246 1969–1971241044 Bozeman Montana SU 45.67 111.05 1482 5.9 487 3401 1951–1969260800 Beowawe 49 S U of N Ranch 39.90 116.58 1751 8.4 263 2208 1981–1998290417 Animas 3 ESE 31.95 108.82 1348 15.3 305 3427 1967–1979310375 Aurora 6 N 35.38 76.78 6 16.3 1273 10,845 1973–2014340184 Altus Dam 34.88 99.30 467 14.8 735 7927 1979–2014391076 Brookings 2 NE 44.32 96.77 503 6 634 9470 1953–2014410174 Alpine 30.35 103.67 1360 16.5 438 3550 1971–1982410613 Beaumont Resrch Ctr 30.07 94.28 9 19.3 1471 12,181 1979–2014410639 Beeville 5 NE 28.45 97.70 76 20.6 832 9254 1979–2014420336 Arches National Park HQS 38.62 109.62 1259 14 250 3666 1980–2004450587 Bellingham 3 SSW 48.72 122.52 9 10.3 964 4947 1985–2014481000 Boysen Dam 43.42 108.20 1467 7.6 245 3772 1951–1976Hydrology 2017, 4, 5 5 of 24l , , Figure 1. Location of 28 stations with pan evaporation measurements where it was possible to reconstruct meteorological variables (either from measurements at the local climate station, or from climate stations within 50 km of the pan evaporation measurements). The map also identifies the six forest ecozones for which mean trend estimates of potential evaporation from model predictions were reported. A pan coefficient, kpan, was applied to measurements of evaporation from Class A pans to account for heat transfer through the sides and bottom of the pans. Early studies in the U.S. suggested a coefficient of kpan = 0.70 [23] and this was also adopted in estimates of open-water evaporation [22]. Another study for pan measurements in the conterminous U.S. reports a value of kpan = 0.77 [45], while irrigation manuals report a value of kpan = 0.80 for moderate wind and moderate humidity [56]. In this study, we applied the mid-range value, kpan = 0.77. 2.2. Models of Potential Evaporation Although temperature-based models may have practical value, they have limitations that have been addressed in detail elsewhere [21,25,57]; thus we restrict focus on methods of estimation that explicitly consider a radiative component. Penman’s model [27] combines radiative and aerodynamic drivers of evaporation: ܧ୔ୣ୬ =ݏ(ܴ୬ − ܩ) + ߛ ݂(ܷ)(݁ୱ − ݁ୟ)ݏ + ߛ1ߣ୐ (1) where s is the slope of the saturation vapour pressure-over-temperature relationship for air (kPa·K−1), Rn is net radiation (W·m−2), G is the ground heat flux (W·m−2), γ is the psychrometric constant (kPa·K−1), es is saturation vapour pressure (kPa), ea is actual vapour pressure (kPa), f(U) is the wind function, here given by f(U) = 1.38 U + 1.31 [17], and ߣL is the latent heat of vapourization (J·kg−1). The equilibrium rate of evaporation was calculated according to [1,17,19,20]: ܧ୉୕ =ݏ(ܴ୬ − ܩ)ݏ + ߛ1ߣ୐ (2) Figure 1. Location of 28 stations with pan evaporation measurements where it was possible toreconstruct meteorological variables (either from measurements at the local climate station, or fromclimate stations within 50 km of the pan evaporation measurements). The map also identifies thesix forest ecozones for which mean trend estimates of potential evaporation from model predictionswere reported.pan coefficient, kpan, was applied to measurements of evaporation from Class A pans to accountfor heat transfer through the sides an bottom of the pans. Early studies in the U.S. suggested acoefficient of kpan = 0.70 [23] and this as als te i esti tes f e - ter e r ti [22].t r st for pan measurements in the conterminous U.S. reports a value of kpan = 0.77 [45],while irrigation manuals report a value of kpan = 0.80 for moderate wind and moderate humidity [56].In this study, w applied the mid-r ng value, kpan = 0.77.2.2. odels of Potential EvaporationAlthough te perature-based odels ay have practical value, they have li itations that havebeen addressed in detail else here [21,25,57]; thus e restrict focus on ethods of esti ation thatexplicitly consider a radiative co ponent. Pen an’s odel [27] co bines radiative and aerodyna icdrivers of evaporation:EPen =s(Rn − G) + γ f (U)(es − ea)s + γ1λL(1)where s is the slope of the saturation vapour pressure-over-temperature relationship for air (kPa·K−1),Rn is net radiation ( ·m−2), G is the ground heat flux ( ·m−2), γ is the psychrometric constant(kPa·K−1), es is saturation vapour pressure (kPa), ea is actual vapour pressure (kPa), f (U) is the windfunction, here given by f (U) = 1.38 U + 1.31 [17], and λL is the latent heat of vapourization (J·kg−1).The equilibrium rate of evaporation was calculated according to [1,17,19,20]:EEQ =s(Rn G)s + γ1λL(2)Hydrology 2017, 4, 5 6 of 24where G, s, γ, and λL are as defined above. The Priestley–Taylor model of potential evaporation wasthen calculated from [29]:EPT = αEEQ. (3)Above wet surfaces, they argued that α = 1.26 yields a strong predictor of potential evaporationat the regional scale.A ‘reference’ rate of evaporation (ERef) is commonly defined for a specific type of surfacevegetation. It is generally defined as an intermediate variable between PE and actualevapotranspiration. In essence, it indicates the maximum rate of actual evapotranspiration [19],or the expected rate of actual evapotranspiration at maximum hydraulic conductivity of plants andsoil. For example, the reference evapotranspiration for a hypothetical surface of vegetation with acharacteristic height and surface conductance can be calculated from the Penman–Monteith equationas [35]:ERef =s(Rn − G) + ρacpga(es − ea)s + γ(1 + gags) 1λL(4)where ρa is the density of air (kg·m−3), and cp is the specific heat of the air (J·kg−1·K−1), gais aerodynamic conductance (m·s−1), and gs is surface conductance (m·s−1). In reality, surfaceconductance is highly variable across forest types and environments, yet here we demonstrate thebehaviour of predictions from the Penman–Monteith equation when ga and gs are set as constants.For example, for a reference grassland, conductance values are commonly set to ga = 0.010 m·s−1 andgs = 0. 014 m·s−1. In so doing, the reference evapotranspiration rate specifically for grasslands (ERG)represents the expected rate of actual evaporation in the absence of water stress.Values of aerodynamic and surface conductance for the reference grassland are not representativeof forest canopies [58,59]. We, therefore, considered application of the Penman–Monteith model tosimulate reference values for forest canopies. Variation in forest canopy structure affects the conversionof sensible heat from the surrounding air into latent heat [60–63]. As surface roughness increases, sodoes the potential extraction of sensible heat and the surface is said to be increasingly ‘coupled’ tothe atmosphere. Conversely, as surface roughness decreases, less sensible heat can be extracted fromthe atmosphere and the surface is said to be increasingly ‘uncoupled’ from the atmosphere [28,60,64].We considered two parameterizations, including forest canopies that are moderately coupled tothe atmosphere (ERFu) and those that are well coupled to the atmosphere (ERFc). For canopies thatare moderately coupled to the atmosphere, we set aerodynamic conductance as ga = 0.058 m·s−1,consistent with average conditions recently reported for coniferous, deciduous, and mixed-woodforests using eddy-covariance systems [31]. For canopies that are well coupled to the atmosphere, weset aerodynamic conductance as ga = 0.150 m·s−1, more consistent with the value commonly assumedin models [59,65]. In choosing to focus on these two levels of aerodynamic conductance for forestcanopies, it was not intended to advocate either one, but rather to understand the implications of eachlevel. In both estimates for reference forests, we applied an equivalent value of surface conductance,gs = 0. 010 m·s−1, reflecting moderately high use of water under optimal conditions [58].As a compliment to the PE models, we additionally developed a statistical model that was fitagainst daily measurements of EPan. This estimate was developed by applying a linear mixed-effects(LME) model with a combination fixed and random effects at the level of each station. We assumedthat daily values of EPan would be related to four variables, including downward global shortwaveradiation, air temperature, vapour pressure deficit, and wind speed, yielding:EPanFit = (β0 + b0, i) + (β1 + b1,i)Rsgd + (β2 + b2,i)Ta + (β3 + b3,i)Da + (β4 + b4, i)U (5)where β0...β4 are fixed effects and b0,i...b4,i, are random effects applied to i = 1...n stations. To the extentthat the assumed drivers and the assumption of additive linear effects can explain PE, this estimatealso reflects the upper limit of how well models should perform. In all cases, the response variable,Hydrology 2017, 4, 5 7 of 24EPan, and the predictor variables, Rsgd, Ta, Da, and U, were standardized (i.e., ‘z-scored’) by subtractingthe sample mean and dividing by the sample standard deviation. The resulting coefficients expressthe magnitude of change in the response variable (in terms of number of standard deviations) that isimposed by an increase of one standard deviation in each predictor variable. All models were fittedin R (R Development Core Team, 2015) using the lmer function from lme4 package [66]. As interestin EPanFit was to describe the behavior of the calibrated model rather than as a predictive tool, weevaluated the performance of Equation (5) based simply on the marginal coefficient of determination(i.e., the variance explained by the fixed effects alone) and the standard error of the coefficients.For all models, values of PE were converted to report flux densities in units of mm·day−1. Unlessotherwise specified, we focused on reporting mean daily PE over a consistent period during the warmseason when most forest regions of North America are above freezing, defined as the months fromMay to September.2.3. Comparison between Models and Pan EvaporationComparison between models and pan evaporation was first performed on a daily time scale byregressing model predictions of daily total PE with daily total EPan. For each comparison, we producedscatterplots and reported the coefficient of determination (R2), the model efficiency coefficient (MEC),the root mean squared error (RMSE), and the mean error (ME) to describe overall model performance.Model comparisons were then performed on a monthly time scale, as the ability to predictmonthly mean PE has broader practical value in applied research. Monthly mean values of daily totalpan evaporation were calculated from the daily values, only considering months with no missingdaily measurements. Monthly mean values of the predictor variables were calculated at each station.To understand the general integrity of the measurements, monthly measurements of pan evaporationwere plotted against individual monthly mean values for each of the driving variables. Scatterplotsand performance statistics were reported for predictions of monthly mean potential evaporation basedon EPen, EPT, ERG, ERFu, ERFc, and EPanFit.2.4. Analysis of Variance AttributionThe relative importance of the drivers was inferred from estimates of sensitivity assessed firston a daily time scale. For the fitted model, the sensitivity of drivers was simply inferred from thefixed effects coefficients of the LME model. For example, the sensitivity of daily EPan to daily Rsgdwas inferred from β1 in Equation (5). To produce comparable sensitivities from the models of PE,we fitted LME models (identical to Equation (5)) to daily values of EPen and EPT, ERG, ERFu, and ERFc.By comparing coefficients, we were able to infer whether attribution in the models was consistent withthat of statistical fits to pan evaporation measurements. We additionally performed the same analysis,only for monthly mean values. That is, we re-fit LME models to monthly mean EPan to produce asecond (monthly) version of EPanFit, as well as sets of coefficients for monthly mean EPen, EPT, ERG,ERFu, and ERFc. The coefficients were compared to understand any differences in attribution betweenmethods. As the coefficients were developed at daily and monthly time scales, the analysis was able toassess whether attribution was scale-dependent.2.5. Model ApplicationTo apply the models across North America, monthly climate data were assembled for a 1 kmgrid covering North America between 1951 and 2014. Where possible, the database was developedby compiling existing datasets. Other variables, developed from raw data sources, are describedin more detail below. For each variable, processing and storage of grids was performed separatelyfor long-term monthly normals and time series of monthly anomalies for 1951–2014. The formerreflect twelve grids indicating the mean monthly values for January to December for the base period1971–2000 for each variable. The latter (anomalies) indicate the monthly deviations from the 1971–2000normal for all months during 1951–2014.Hydrology 2017, 4, 5 8 of 24Transient monthly variation in Rsgd was taken from the National Centers for EnvironmentalPredition/National Center for Atmosphere Research (NCEP/NCAR) global reanalysis [67]. Similarto the methodology we applied previously [68], a first approximation of monthly normal Rsgd wasdeveloped by combining the synoptic and zonal variation captured by values from the North AmericanRegional Reanalysis [69] with topographically-driven variation captured by values from the ArcGISSolar Radiation Tool [70]. In the first step, estimates of Rsgd from North American Regional Reanalysis(NARR) were downloaded and resampled from the original resolution of 32 km down to a resolutionof 1 km using cubic interpolation. In the second step, potential Rsgd was calculated using the ArcGISSolar Radiation Tool with default settings (eight zenith divisions, eight azimuth divisions, uniformsky diffuse model type, a diffuse proportion of 0.3 and a transmittivity of 0.5) for each month at 1 kmresolution. The tool was run with a digital elevation model developed by the U.S. Geological Surveythe HYDRO1k project (https://lta.cr.usgs.gov/HYDRO1K). The ArcGIS Solar Radiation Tool doesnot take latitude as an input argument so the function was executed for 0.5◦-wide zonal bands with0.15◦ overlapping buffers on the north and south sides so as to avoid edge effects. In the third step,the two datasets were combined to produce a first approximation of Rsgd by correcting the ArcGISSolar Radiation Tool estimates so that they were equivalent to the NARR estimates at the original scaleof the NARR dataset (i.e., 32 km). This was achieved by applying a two-dimensional 32-km averagefilter (“fspecial” and “filter2” functions, The Mathworks Matlab 7.9.0) to the difference between NARRand ArcGIS estimates and adding it to the estimates of potential radiation from ArcGIS. The intendedoutcome was an unbiased estimate of the initial NARR dataset with added subgrid-scale variationcaused by topographic effects. That is, the ArcGIS Solar Radiation Tool accounts for expected variationin irradiance caused by earth–sun geometry and shading by neighbourhood terrain, while NARRestimates express the synoptic-scale variation.To expand on previous accuracy assessment [68], we compared the first approximation withmeasurements of monthly Rsgd across North America provided by Environment Canada and from theNational Solar Radiation Data Base for the United States produced by the National Renewable EnergyLaboratory [54]. The comparison indicated that initial NARR estimates systematically overestimatedRsgd by 15%–30% depending on month of the year. We made the assumption that the estimatesfrom station networks were closer to the true values and developed a second approximation of Rsgdthat corrected the first approximation to better match the station values. This was achieved in twosteps: First, the initial monthly biases between NARR and the station observations (an average offsetof 1.73 MJ·day−1) was subtracted from the first approximation. Second, the remaining residual errors ateach station were interpolated to the grid using a polynomial surface fit with spatial variation definedby the errorgram and solved using singular value decomposition [71,72]. The resulting continuous gridof differences suggested that NARR dataset overestimated Rsgd in the southwest U.S. and large regionsof central Canada and underestimated it in the northeast conterminous U.S. The second (and final)approximation of Rsgd was then developed by adding the residual error field to the first approximation.To validate the second approximation of Rsgd, nearest grid cell estimates were comparedwith 25 independent Fluxnet tower observations. Tower observations spanned different years between1998 and 2010 so a proper validation of the base period was not possible. Deviations from thelong-term mean from the NARR dataset, corresponding with the period of comparison with eachtower were, therefore, added to the second approximation grid to ensure the comparison was forestimates representative of the same period. Removing regional biases between the first approximationand Environment Canada/NREL station networks increased the level of explained variance from75% to 85% and reduced the root-mean squared error (RMSE) from 3.96 to 0.97 MJ·m−2·day−1.Mean relative error between the second approximation and the tower observations was consistentlybelow 10% for each month of year.A first approximation of monthly normal Ta was compiled from the ClimateNA dataset [73,74].Transient variability in monthly temperature was derived from three sources, including EnvironmentCanada [75], the U.S. Historical Climate Station Network (USHCN) [76] for the conterminous U.S.Hydrology 2017, 4, 5 9 of 24and from the Global Summary of the Day (GSOD) product [55] for Alaska. For Canada, daily climaterecords were extracted for 697 stations that included all stations listed in Environment Canada'sreference climate station network, all stations in the Adjusted and Homogenized Canadian ClimateData archive [77], and some additional remote stations with long records. Canadian station timeseries were scanned for inhomogeneities by eye. Non-climatic steps (that were obviously not natural)were identified and the shorter of the resulting time series segments were adjusted to eliminate thedifference. Seven stations were removed in the process due to poor record quality. Gaps in the originalrecords were filled using regression analysis to relate each station record with that of four nearestneighbours. Monthly mean air temperature was then calculated from the gap-filled records of dailyminimum and maximum air temperatures.Transient monthly variation of Da was then derived from resampling of the Climate ResearchUnit estimates of vapour pressure [78] and saturation vapour pressure (ea) derived from Ta. A firstapproximation of long-term mean monthly Da was developed based on a function of minimumand maximum temperature [79,80], with the temperature data coming from ClimateNA [74].As with radiation, the first approximation was compared with measurements compiled from theGlobal Summary of the Day dataset [55]. Environment Canada measurements were also addedafter noticing that not all Environment Canada data appears to have been entered into the GSODdatabase. The difference between nearest grid cells of the first approximation and stations was thencalculated and several outliers were censored from the station record. The differences were theninterpolated to produce continuous maps of the error, showing considerable differences over thesouthern conterminous U.S. Adding the difference field to the first approximation to produce thesecond approximation led to improved agreement with independent observations of vapour pressuredeficit at Fluxnet towers. The level of explained variance in warm-season mean vapour pressure deficitincreased from 54% to 76% and the RMSE decreased from 1.79 to 1.40 hPa.Values of s, γ, and λL were calculated as functions of air temperature. The above preparationof climate input grids included Rsgd, yet models of PE required input grids of Rn. Net radiationwas derived from incident shortwave radiation in two steps. First, we converted incident shortwaveradiation to net shortwave radiation to account for differences in mean albedo of various forest surfacetypes. This was achieved by producing a lookup table linking surface albedo values for land coverclasses. Land cover classes were derived from the 2002 version of the North American Land CoverCharacteristics map produced by the National Center for Earth Resources Observation and Science,United States Geological Survey [81]. Average surface albedo values for each land cover class weredrawn from a previous synthesis of values [82]. Net shortwave radiation (Rsgn) was then calculatedby multiplying Rsgd by (1-albedo). In a second step, we developed linear mixed effects models thatrelated monthly mean observations of net shortwave radiation (W·m−2) to observations of monthlymean net radiation (W·m−2). In the model, we included a random effect for the intercept and slope.Data were compiled from 14 forested Fluxnet station sites across North America from the Oak RidgeNational Laboratory Fluxnet Database [83] (Sites codes: CA-Ca1, CA-Ca2, CA-Ca3, US-Ha1, US-MRf,US-Me2, US-NR1, CA-Qfo, CA-Oas, CA-Obs, CA-Ojp, US-MMS, US-Wrc). This gave a sample size of1054 monthly observations that were representative of deciduous and evergreen forest stands fromboreal and temperate forest biomes. The global model yielded a relationship, Rn = 0.830Rsgn − 31.364.The slope and intercept had standard errors of 0.006 and 2.551, respectively, and were both significant(p < 0.001). The coefficient of determination, excluding random effects, was R2 = 0.95, and the rootmean squared error was 13.5 W·m−2. We additionally tested whether differences existed betweendeciduous and evergreen stands, or between boreal and temperate forest biomes. We found thatparameters were not statistically different between boreal and temperate stands (p = 0.23), but wedid find that parameters differed between deciduous and evergreen stands (p = 0.009). The modelfor deciduous stands yielded a relationship, Rn = 0.882Rsgn − 39.500. The slope and intercept hadstandard errors of 0.011 and 2.300, respectively, and were both significant (p < 0.001). The model forevergreen stands yielded a relationship, Rn = 0.817Rsgn − 29.512. The slope and intercept had standardHydrology 2017, 4, 5 10 of 24errors of 0.031 and 5.894, respectively, and were both significant (p < 0.001). Parameters for deciduousand evergreen forest were linked to those corresponding land cover classes in the USGS land coverclassification map, while parameters from the global model were applied in all other land cover classes,including mixed forest.Although Class A evaporation pans and natural land surfaces may exhibit diurnal and seasonalcycles in heat capacity, for simplicity, we assumed that the ground heat flux was zero at a decadaltime scale. Evaluation of the long-term mean ground heat flux from 12 of the aforementioned Fluxnetsites located in forest ecosystems suggested that this assumption may lead to an overestimate in netradiation of 0.11 percent.We summarized long-term warm-season mean values of PE from each method by forestregion according to ecozones identified in the unified Level 1 ecozones of Canada and the U.S.produced by the Commission for Environmental Cooperation [84]. We additionally summarizedthe zonal variation in warm-season mean values over the range of latitude in North America.We summarize time trends by forest region. There was one station in the western temperate forestecozone, two stations in the Montane Cordillera ecozone, one station in the western boreal ecozone,two stations in the Boreal Plain ecozone, four stations in the eastern boreal ecozone, and six stationsin the eastern temperate forest ecozone. There is a high degree of uncertainty in values of ga andgs at this scale. We therefore investigated the sensitivity of trends to the parameterization of thePenman–Monteith model by calculating time trends in the resulting values of PE for a broad range invalues of ga and gs and then applied multivariate regression analysis to derive a general relationshipbetween time trends in PE and values of ga and gs used in the Penman–Monteith model.3. Results3.1. Daily Model ComparisonAll candidate predictor variables in the LME model fit to daily EPan were significant (p < 0.001) andfixed effects coefficients ranged from 0.402 for Rsgd to 0.126 for U (Table 2). Variance inflation factorsfor the predictor variables were all well below 10.0 so we assumed multicollinearity did not confoundinterpretation of the coefficients. Large systematic differences, perhaps reflecting measurement errorand instrumentation, was implied by substantial variation in the random effect on the intercept,ranging from −0.632 to 0.429 (Table 2). Relationships between models of PE and EPan exhibitedlow precision (Figure 2). Coefficients of determination were similar across models. EPen exhibiteda model efficiency coefficient of 0.27, EPT exhibited a model efficiency coefficient of 0.33, ERFu andERFc both exhibited model efficiency coefficients close to zero (i.e., no better than applying the meanas a predictor), while EPanFit exhibited a model efficiency coefficient of 0.40. In terms of overall bias,EPen exhibited a mean error of 0.5 mm·day−1 (10%) (i.e., EPen overestimated EPan by 10 percent), ERG,ERFu, and ERFc exhibited larger biases (i.e., mean errors of 0.3, −0.1, and 0.1 mm·day−1 (7%, −3% and2%), respectively, while EPT and EPanFit both exhibited a mean error of 0.0 mm·day−1 (−1%).Table 2. Fixed and random effects for parameters in the linear mixed effects model of daily panevaporation measurements. The marginal coefficient of determination (i.e., without applicationof the random effects) was R2 = 0.41). The mean (observation minus prediction) error was−0.04 ± 0.01 mm·day−1 (−0.96 ± 0.31%). The root mean squared error was 1.88 mm·day−1.Intercept Rsgd Ta Da UFixed effectsEstimate 0.033 0.402 0.209 0.156 0.126Standard error 0.060 0.032 0.023 0.026 0.031Hydrology 2017, 4, 5 11 of 24Table 2. Cont.Intercept Rsgd Ta Da URandom effectsWhitehorse A −0.030 −0.008 −0.179 0.149 −0.020Norman Wells 0.295 0.062 −0.046 0.038 0.263Resolute Cars 0.230 −0.064 0.013 0.208 0.029Edmonton Stony Plain 0.406 0.091 −0.059 0.001 0.434Winnipeg Richarson Int’l A 0.299 0.079 0.039 −0.016 0.211Churchill A 0.101 −0.093 0.027 0.156 −0.020Pierre Elliott Trudau Intl A −0.020 0.052 0.003 −0.057 −0.075La Grande IV A 0.379 0.057 0.072 0.039 0.297Nitchequon −0.187 −0.036 0.012 −0.032 −0.075Schefferville A 0.136 0.034 0.048 0.055 0.237Antioch Pump Plant 3 0.082 0.275 −0.040 −0.164 −0.144Bonny Dam 2 NE 0.394 0.021 0.096 −0.102 −0.006Allatoona Dam 2 −0.632 −0.262 −0.034 −0.033 −0.031Arrowrock Dam −0.052 0.105 0.028 −0.208 −0.181Buckhorn Lake −0.565 −0.211 −0.069 −0.086 −0.035Alberta Ford For Cen −0.611 −0.166 −0.028 −0.007 −0.111Bozeman Montana SU −0.208 0.043 −0.099 −0.029 −0.228Beowawe 49 S U of N Ranch 0.049 0.010 0.161 −0.096 −0.011Animas 3 ESE 0.456 0.376 −0.058 −0.124 −0.073Aurora 6 N −0.077 −0.162 0.051 0.148 0.003Altus Dam −0.264 −0.268 0.083 0.219 0.070Brookings 2 NE 0.136 −0.120 0.106 −0.104 −0.110Alpine −0.137 −0.029 −0.307 0.177 −0.009Beaumont Resrch Ctr 0.037 0.121 −0.087 −0.126 −0.121Beeville 5 NE −0.299 −0.136 0.062 0.073 0.042Arches National Park HQS 0.429 0.378 −0.078 −0.066 −0.189Bellingham 3 SSW −0.352 −0.155 0.056 0.181 0.038Boysen Dam 0.006 0.008 0.225 −0.195 −0.1833.2. Monthly Model ComparisonMonthly variability of EPan ranged from 1.0 to 11.1 mm·day−1, with no major outliers (Figure 3).Of the individual meteorological drivers, the relationship was strongest between EPan and Rsgd andthis relationship showed some indication of curvature (Figure 3a). Relationships with Ta and Da wereweaker (Figure 3b,c), while there was no relationship between monthly EPan and U (Figure 3d).Consistent with fits previously reported at the daily scale, all candidate predictor variables inthe LME model fit to monthly EPan were significant (p < 0.001) (monthly LME model coefficientsnot shown). The marginal coefficient of determination for the monthly fit to pan evaporation wasR2 = 0.68, with a mean (observation minus prediction) error of −0.16 ± 0.06 mm·day−1, and a rootmean squared error of 0.93 mm·day−1.Models listed in decreasing order of performance (i.e., agreement with monthly EPanmeasurements based on model efficiency coefficients), were EPanFit (MEC = 0.67), EPT (MEC = 0.64),EPen (MEC = 0.53), ERG (MEC = 0.31), ERFu (MEC = 0.10), and ERFc (MEC = −0.20) (Figure 4). The rateof evaporation for a reference forests with canopies that are moderately coupled to the atmospheresubstantially underestimated EPan (Figure 4d). By design, these values were not expected to beunbiased expressions of PE, but we plotted them here to compare precision and explained variancefor context. EPanFit slightly outperformed EPT, explaining 68% of the variance, with a mean errorof −3% (Figure 4f). It is interesting to note that, despite slightly greater overall performance of EPanFitand EPT, these models still had limited capacity to reproduce the highest monthly values of EPan above7.0 mm·day−1 (Figure 4b,f).Hydrology 2017, 4, 5 12 of 24Hydrology 2017, 4, 5 12 of 24 Figure 2. Comparison between observations of daily total pan evaporation (EPan) (after applying the pan coefficient of 0.77) and models of potential evaporation based on: (a) Penman’s model (EPen); (b) the Priestley–Taylor model (EPT); (c) reference rate for grassland based on the Penman–Monteith model (d) reference rate for forests canopies that are moderately coupled to the atmosphere (ERFu); (e) reference rates for forest canopies that are well coupled to the atmosphere (ERFc); (f) linear mixed-effects model fit to daily pan evaporation measurements (EPanFit). Linear curves indicate the ordinary least squares fits. Statistics include: coefficient of determination (R2), model efficiency coefficient (MEC, values of 0.0 or <0.0 indicate equivalent, or less predictive power than applying the mean of observations, respectively, while a value of 1.0 indicates a perfect predictor), root mean squared error (RMSE, mm·day−1), and mean error (ME, mm·day−1 or %). Figure 2. Comparison between observations of daily total pan evaporation (EPan) (after applying thepan coefficient of 0.77) and models of potential evaporation based on: (a) Penman’s model (EPen);(b) the Priestley–Taylor model (EPT); (c) reference rate for grassland based on the Penman–Monteithmodel (d) reference rate for forests canopies that are moderately coupled to the atmosphere (ERFu);(e) referenc rates for forest canopies that are w ll coupled to the atmosphere (ERFc); (f) linearmix d-effects mod l fit to daily pan evaporation measurements (EPanFit). Linear curves indicatethe ordinary least squares fits. Statistics include: coefficient of determination (R2), model efficiencycoefficient (MEC, values of 0.0 or <0.0 indicate equivalent, or less predictive power than applyingthe mean of observations, respectively, while a value of 1.0 indicates a perfect predictor), root meansquared error (RMSE, mm·day−1), and mean error (ME, mm·day−1 or %).Hydrology 2017, 4, 5 13 of 24Hydrology 2017, 4, 5 13 of 24 Figure 3. Relationships between monthly mean pan evaporation (EPan) and climate variables, including incident global shortwave radiation (Rsgd), air temperature (Ta), vapour pressure deficit (Da), and wind speed (U). Figure 4. Comparison between observations of monthly mean pan evaporation (EPan) and model predictions of potential evaporation (PE) based on: (a) Penman’s model; (b) the Priestley–Taylor model; (c) reference rates for grasslands based on the Penman–Monteith model; (d) reference rates for forests based on the Penman–Monteith model assuming the canopy is moderately coupled to the atmosphere; (e) reference rates for forests based on the Penman–Monteith model assuming the canopy is well coupled to the atmosphere; (f) linear mixed-effects model fit to monthly pan evaporation (EPanFit). Figure 3. Relationships betw en monthly mean pan evaporation (EPan) and climate variabl s, includingincident global shortwave radiation (Rsgd), air temperature (Ta), vapour ressu deficit (Da), andwi speed (U).Hydrology 2017, 4, 5 13 of 24 Figure 3. Relationships betw en monthly mean pan evaporation (EPan) nd cl mate variables, including incident global shortwave radiation (Rsgd), air temperature (Ta), vapour pre sure deficit (Da), and win sp ed (U). Figure 4. Comparison betw en observations of monthly mean pan evaporation (EPan) and model predictions of potential evaporation (PE) based on: (a) Penman’s model; (b) the Priestley–Taylor model; (c) reference rates for gra slands based on the Penman–Monteith model; (d) reference rates for forests based on the Penman–Monteith model a suming the canopy is moderately coupled to the atmosphere; (e) reference rates for forests based on the Penman–Monteith model a suming the canopy is well coupled to the atmosphere; (f) linear mixed-e fects model fit to monthly pan evaporation (EPanFit). Figure 4. Comparison between observati f o t l e Panpredict ons of poten ial evaporation (PE) based on: (a) Penma ’s model; (b) the Priestley–Taylor model;(c) ref ren e rates for gra lands based on the Pe man–Monteith model; (d) r ference rat s fo forestsbased on the Penman–Mo teith m del assuming the canopy is moderately c upled to the atmosphere;(e) referenc rat s for forests ba ed on the Penman–Monteith odel assuming th canopy is wellcoupled to the atm sph re; (f) lin r ixed- ffects model fit to monthly pan evaporation (EPanFit).Hydrology 2017, 4, 5 14 of 243.3. Attribution of Variance in PE to Meteorological DriversIn decreasing order of sensitivity, daily variability of EPan was controlled by Rsgd, Ta, Da, and U(Figure 5). On a monthly time scale, the controls of Ta and Da on EPan dissipated and were compensatedby increases in sensitivity to Rsgd, while control by U was no longer significant.Hydrology 2017, 4, 5 14 of 24 3.3. At ribution of Variance in PE to M teorological ri In decreasing order of sensitivity, daily variability of EPan was controlled by Rsgd, Ta, Da, and U (Figure 5). On a monthly time scale, the controls of Ta and Da on EPan dissipated and were compensated by increases in sensitivity to Rsgd, while control by U was no longer significant. Figure 5. The sensitivity of potential evaporation (PE) to driving variables on a daily and monthly time scale. Left-to-right (red bars): Fixed effects from linear mixed-effects model fits to daily pan evaporation (EPanFit); fixed effects from model fits to daily predictions of PE from Penman’s model; fixed effects from model fits to daily predictions of PE from the Priestley–Taylor model; fits to daily reference evaporation from grasslands based on the Penman–Monteith model; fits to daily reference evaporation from forest canopies that are moderately coupled to the atmosphere based on the Penman–Monteith model; fits to daily reference evaporation from forest canopies that are well coupled to the atmosphere based on the Penman–Monteith model. Blue bars then repeat through the same estimates, only for fits to monthly instead of daily time scale. On a daily time scale, EPen was most strongly driven by Rsgd and secondarily by Da and U, while Ta exerted a low positive effect. There were no significant differences in sensitivity between daily and monthly time scale for EPen. EPT exhibited the strongest sensitivity to Rsgd, with secondary dependence on Ta. On a monthly time scale, control on EPT shifted slightly from Rsgd to Ta. Despite differences in aerodyne mic, the behaviour of ERFu and ERFc only showed modest differences in sensitivity to Da. Both exhibited strong dependence on Da, followed by Rsgd and weak positive dependence on Ta. Hence, none of the models perfectly matched the variance attribution of monthly EPanFit; sensitivity of monthly EPT to Ta was too high, while sensitivities of monthly EPen, ERFu, and ERFc to Da were too high. 3.4. Spatial and Temporal Variability Differences in the spatial variation of PE were evident from visual inspection of the mean warm-season values over the 1971–2000 base period (Figure 6). Patterns for Penman’s model and the reference grassland PE were similar, with the exception that EPen showed slightly greater values of PE across regions of southern Canada (Figure 6a,c). The Priestley–Taylor model indicated relatively low PE in arid regions of the southwest U.S. compared to other models and also showed slightly greater levels of PE across southern Canada relative to ERG (Figure 6b). Reference forest values showed more prominent zonal gradients and the differences between spatial patterns of PE for moderately- and well-coupled canopies were much less than that between reference forest rates and other models (Figure 6d,e). All models indicated decreases in warm-season mean daily PE as latitude increased from 25 to 80 degrees north (Figure 7). Zonal gradients were greatest for estimates of reference evaporation using the Penman–Monteith equation, and similar among other models. At the scale of ecozones, historical values of warm-season mean daily PE ranged from 2.4 to 4.5 mm·day−1 (Table 3). Mean historical values of warm-season PE ranged from 2.9 mm·day−1 (ERG) to 3.8 mm·day−1 (EPanFit) (Table 3). Figure 5. The sensitivity of potential evaporation (PE) to driving variables on a daily and monthly timescale. Left-to-right (red bars): Fixed effects from linear ixed-effects model fits to daily pan eva oration(EPanFit); fixed effects from model fits to daily predictions of PE from Penman’s model; fixe effectsfrom mo l fits to daily predictions of PE fr m the Priestley–Taylor model; fits to daily referenceevaporation from grasslands based on the Penman–Monteith model; fits to daily reference evaporationfrom forest canopies that are moderately coupled to the atmosphere based on the Penman–Monteithmodel; fits to daily reference evaporation from forest canopies that are well coupled to the atmospherebased on the Penman–Monteith model. Blue bars then repeat through the same estimates, only for fitsto monthly instead of daily time scale.On a daily time cale, EPen was most strongly driven by Rsgd and secondarily by Da and U,while Ta exerted a low positive effect. There were no significant differences i sensitivity betweendaily and monthly time scale for EPen. EPT exhibited the strongest sensitivity to Rsgd, with secondarydependence on Ta. On a monthly time scale, control on EPT shifted slightly from Rsgd to Ta. Despitedifferences in aerodyne mic, the behaviour of ERFu and ERFc only showed modest differences insensitivity to Da. Both exhibited strong dependence on Da, followed by Rsgd and weak positivedependence on Ta. Hence, none of the models perfectly matched the variance attribution of monthlyEPanFit; sensitivity of monthly EPT to Ta was too high, while sensitivities of monthly EPen, ERFu,and ERFc to Da were too high.3.4. Spatial and Temporal VariabilityDiffe ences in the patial vari tion of PE were vid nt from visual inspe tion of h meanwarm-season values over the 1971–2000 base period (Figure 6). Patterns for Penman’s model and thereference grassland PE were similar, with the exception that EPen showed slightly greater values ofPE across regions of southern Canada (Figure 6a,c). The Priestley–Taylor model indicated relativelylow PE in arid regions of the southwest U.S. compared to other models and also showed slightlygreater levels of PE across southern Canada relative to ERG (Figure 6b). Reference forest values showedmore prominent zonal gradients and the differences between spatial patterns of PE for moderately-and well-coupled canopies were much less than that between r ference forest rates and other models(Figure 6d,e).All models indicated decreases in warm-season ean daily PE as latitude increased from 25 to80 degrees north (Figure 7). Zonal gradients were greatest for estimates of reference evaporation usingthe Penman–Monteith equation, and similar among other models. At the scale of ecozones, historicalvalues of warm-season mean daily PE ranged from 2.4 to 4.5 mm·day−1 (Table 3). Mean historicalvalues of warm-season PE ranged from 2.9 mm·day−1 (ERG) to 3.8 mm·day−1 (EPanFit) (Table 3).Hydrology 2017, 4, 5 15 of 24 Figure 6. Comparison between long-term (1971–2000) mean daily potential evaporation (PE) overthe warm season (May–September). Units of PE are expressed as mm·day−1: (a) Penman’s model(EPen); (b) Priestley–Taylor model (EPT); (c) reference grassland (ERG) based on the Penman–Monteithmodel; (d) reference forest with moderately coupled canopy (ERFu) based on the Penman–Monteithmodel; (e) reference forest with coupled canopy (ERFc) based on the Penman–Monteith model; (f) linearmixed-effects model fit to monthly mean pan evaporation (EPanFit).Hydrology 2017, 4, 5 16 of 24Hydrology 2017, 4, 5 16 of 24 Figure 7. Zonal mean warm-season (May–September) potential evaporation from various models over the land surface of North America. Table 3. Ecozone-average mean warm-season (May–September) variables over the period 1951–2014: Incident global solar radiation, Rsgd (MJ·m−2·day−1), air temperature, Ta (°C), vapour pressure, ea (hPa), vapour pressure deficit, Da (hPa), Penman’s model, EPen (mm·day−1), the Priestley–Taylor model, EPT (mm·day−1), the Penman–Monteith model for a reference grassland, ERG (mm·day−1), the Penman–Monteith model for a reference forest with moderately coupled canopy, ERFu (mm·day−1), and reference forest with well coupled canopy, ERFc (mm·day−1), and statistical fits to pan evaporation measurements EPanFit (mm·day−1). Rsgd Ta ea Da EPen EPT ERG ERFu ERFc EPanFit Boreal, west 16.0 9.6 8.0 4.3 2.7 2.5 2.4 2.6 2.7 3.2 Boreal Plains 17.7 11.7 9.4 4.7 3.1 3.0 2.8 2.8 2.9 3.6 Boreal, east 17.2 11.6 10.4 3.8 2.9 2.9 2.5 2.5 2.4 3.4 Temperate, west 18.1 12.2 10.8 4.0 3.1 3.1 2.7 2.5 2.5 3.6 Montane Cordillera 21.4 11.8 8.0 6.3 4.0 3.8 3.6 3.7 3.9 4.2 Temperate, east 19.5 20.7 19.1 6.2 3.9 4.0 3.6 4.1 3.8 4.5 Mean 18.3 12.9 11.0 4.9 3.3 3.2 2.9 3.0 3.1 3.8 On average over forest ecozones, there was little change in Rsgd from 1951–2014 (Table 4). According to NCEP/NCAR global reanalysis, Rsgd increased considerably over the boreal forest ecozones, and decreased considerably over eastern temperate forest ecozone. The station-based reconstruction of warm-season air temperature suggested an increase of 1.21 °C from 1951–2014 (Table 4). A large positive effect of warming on saturation vapour pressure led to a net increase in vapour pressure deficit of 0.81 hPa despite concomitant increase in vapour pressure of 0.38 hPa. No attempt was made to analyze trends in U. Continuous high-resolution reconstructions of warm-season (May–September) mean PE based on values of EPen, EPT, ERG, ERFu, ERFc, and EPanFit suggested a range between negative and positive time trends from 1951–2014 (Figure 8). The magnitude of time trends were irregularly distributed, but similarities in the pattern were apparent among models. Negative trends were consistently concentrated over the east-central U.S. and pockets of industrialized areas across Canada and Alaska. On average across forest ecozones, the long-term change in estimates of potential evaporation (inferred from least-squares fits of potential evaporation against time for each model) ranged from just 0.08 mm·day−1 for the Priestley–Taylor model, to 0.45 mm·day−1 for the reference forest that is well coupled to the atmosphere (Table 4). Average long-term changes were intermediate for Figure 7. Zonal mean warm-season (May–September) potential evaporation from various models overthe land surface of North America.Table 3. Ecozone-average mean warm-season (May–September) variables over the period 1951–2014:Incident global solar radiation, Rsgd (MJ·m−2·day−1), air temperature, Ta (◦C), vapour pressure, ea(hPa), vapour pressure deficit, Da (hPa), Penman’s model, EPen (mm·day−1), the Priestley–Taylormodel, EPT (mm·day−1), the Penman–Monteith model for a reference grassland, ERG (mm·day−1), thePenman–Monteith model for a reference forest with moderately coupled canopy, ERFu (mm·day−1),and reference forest with well coupled canopy, ERFc (mm·day−1), and statistical fits to pan evaporationmeasurements EPanFit (mm·day−1).Rsgd Ta ea Da EPen EPT ERG ERFu ERFc EPanFitBoreal, west 16.0 9.6 8.0 4.3 2.7 2.5 2.4 2.6 2.7 3.2Boreal Plains 17.7 11.7 9.4 4.7 3.1 3.0 2.8 2.8 2.9 3.6Boreal, east 17.2 11.6 10.4 3.8 2.9 2.9 2.5 2.5 2.4 3.4Temperate, west 18.1 12.2 10.8 4.0 3.1 3.1 2.7 2.5 2.5 3.6Montane Cordillera 21.4 11.8 8.0 6.3 4.0 3.8 3.6 3.7 3.9 4.2Temperate, east 19.5 20.7 19.1 6.2 3.9 4.0 3.6 4.1 3.8 4.5Mean 18.3 12.9 11.0 4.9 3.3 3.2 2.9 3.0 3.1 3.8On average over forest ecozones, there was little change in Rsgd from 1951–2014 (Table 4).According to NCEP/NCAR global reanalysis, Rsgd increased considerably over the boreal forestecozones, and decreased considerably over eastern temperate forest ecozone. The station-basedreconstruction of warm-season air temperature suggested an increase of 1.21 ◦C from 1951–2014(Table 4). A large positive effect of warming on saturation vapour pressure led to a net increasein vapour pressure deficit of 0.81 hPa despite concomitant increase in vapour pressure of 0.38 hPa.No attempt was made to analyze trends in U.Continuous high-resolution reconstructions of warm-season (May–September) mean PE basedon values of EPen, EPT, ERG, ERFu, ERFc, and EPanFit suggested a range between negative and positivetime trends from 1951–2014 (Figure 8). The magnitude of time trends were irregularly distributed,but similarities in the pattern were apparent among models. Negative trends were consistentlyconcentrated over the east-central U.S. and pockets of industrialized areas across Canada and Alaska.On average across forest ecozones, the long-term change in estimates of potential evaporation(inferred from least-squares fits of potential evaporation against time for each model) ranged fromjust 0.08 mm·day−1 for the Priestley–Taylor model, to 0.45 mm·day−1 for the reference forest thatHydrology 2017, 4, 5 17 of 24is well coupled to the atmosphere (Table 4). Average long-term changes were intermediate forPenman’s model and the LME fit to pan evaporation measurements at 0.16 and 0.14 mm·day−1respectively, slightly higher for the reference grassland (0.19 mm·day−1), and high for reference forests(0.37 and 0.45 mm·day−1 for moderately- and well-coupled canopies, respectively).Based on a compilation of the time trends in warm-season PE for reference forests, using a range ofaerodynamic and surface conductance, a general relationship was produced from regression analysis(Figure 9). The analysis confirmed that time trends increased strongly with increasing aerodynamicconductance and weakly with increasing surface conductance.Table 4. Ecozone-average long-term changes in warm-season (May–September) mean daily variablesover the period 1951–2014: Incident global solar radiation, Rsgd (MJ·m−2·day−1), air temperature,Ta (◦C), vapour pressure deficit, Da (hPa), Penman’s model, EPen (mm·day−1), the Priestley–Taylormodel, EPT (mm·day−1), the Penman–Monteith model for a reference grassland, ERG (mm·day−1), thePenman–Monteith model for a reference forest with moderately-coupled canopy, ERFu (mm·day−1),and reference forest with well-coupled canopy, ERFc (mm·day−1), and statistical fits to pan evaporationmeasurements EPanFit (mm·day−1).Rsgd Ta ea Da EPen EPT ERG ERFu ERFc EPanFitBoreal, west 0.09 1.30 0.43 0.63 0.17 0.10 0.19 0.31 0.36 0.15Boreal Plains 0.08 1.33 0.20 1.01 0.22 0.10 0.26 0.48 0.57 0.18Boreal, east 0.40 1.49 0.72 0.65 0.23 0.18 0.24 0.34 0.38 0.21Temperate, west 0.09 1.11 0.37 0.74 0.17 0.10 0.20 0.34 0.41 0.14Montane Cordillera −0.04 1.14 0.19 0.95 0.18 0.09 0.22 0.43 0.52 0.14Temperate, east −0.58 0.86 0.35 0.89 0.00 −0.10 0.06 0.31 0.44 0.04Mean 0.01 1.21 0.38 0.81 0.16 0.08 0.19 0.37 0.45 0.144. DiscussionEcological studies commonly make use of basic expressions of evaporative demand. Here, weexplored the behavior of several common approaches to calculating PE. While such basic expressionsof PE may have strong practical value in studies that attempt to predict climate change impacts onterrestrial ecosystems across large spatial and temporal domains, further work is needed to improvethe representation of spatial variation in meteorological conditions above forests, canopy structureand physiology, and ultimately the prediction of actual evapotranspiration rather than PE or referencerates [52,53].We considered measurements of pan evaporation as one independent indicator of PE. While thereare likely key differences between PE above forests and evaporation pans, the comparison helpedplace the attribution of variance and spatial and temporal variability of PE predictions into context.Upon compiling the pan evaporation data from Canada and the U.S., we discovered that records wereonly rarely accompanied by on-site meteorological variables, so considerable effort was placed onderiving best available values from nearby climate stations, which constitutes an important sourceof error. Although varying conditions at the evaporation pans, and errors associated with site, wind,heating, splashing, likely led to different pan coefficients, we chose to restrict the scope of the study byapplying a constant pan coefficient and including a random effect on the intercept coefficient at thescale of individual pan evaporation stations in the LME models. Applying a mid-range value of thepan coefficient, k = 0.77, led to a mean EPanFit of 3.8 mm·day−1 across North American forests. As thisvalue was 0.5 mm·day−1 greater than the next highest value (EPen), any future consideration of thedataset should likely consider reducing k.Hydrology 2017, 4, 5 18 of 24 Figure 8. Long-term change in warm-season (May–September) mean potential evaporation (PE) from1951–2014 (mm·day−1): (a) Penman’s model (EPen); (b) Priestley–Taylor model (EPT); (c) referencegrassland (ERG); (d) reference forest with moderately coupled canopy (ERFu); (e) reference forestwith well coupled canopy (ERFc); (f) linear mixed-effects model fit to monthly mean pan evaporation(EPanFit). Values were calculated by multiplying the slope coefficient from regression relationshipbetween PE and time by the length of the study period (64 years).Hydrology 2017, 4, 5 19 of 24Hydrology 2017, 4, 5 19 of 24 Figure 9. Dependence of long-term changes in reference rates of evaporation on levels of aerodynamic conductance (ga) and surface conductance (gs) set in the Penman–Monteith equation. Trends are representative of mean warm-season daily reference evaporation over forest regions of North America (z-axis indicates change in mm·day−1 from 1951–2014). Varying sensitivity to Da was a defining difference among the most popular models of PE. The Priestley–Taylor model and the Penman–Monteith model parameterized for conditions above reference forest canopies reflect estimates with absent Da-sensitivity and strong Da-sensitivity, respectively. Any inferences about change in surface hydrology and consequent socio-economic implications over time are, therefore, going to strongly depend on which model is applied. The differentiation in long-term change among models was exacerbated by strong apparent positive trend in historical Da across North America [85]. Interestingly, the different dependence on Da among models also shaped the zonal variation in PE (Figure 7), where estimates of PE derived from reference forest rates actually exceeded estimates of PE with infinite surface conductance in Penman’s model. Specifically, estimates of ERFu and ERFc start to exceed EPen when Da > 18.5 hPa. On average across forests of North America, results indicate that PE increased from 1951–2014. If this is correct, increasing PE must increase the tendency for drought stress. Confidence in the claim of a positive trend is gained from the fact that all models consistently indicated positive trends. However, the claim of positive trends remains uncertain due to reliance on irradiance from NCEP/NCAR reanalysis, which has poorly constrained uncertainty, and by the inability to account for trends in wind speed. These sources of uncertainty will affect all models similarly. If we ignore uncertainty in irradiance and wind trends, it must still be recognized that the range of trend estimates among models of PE includes low estimates, which may have been exceeded by increasing trends in precipitation [86]. If we exclude consideration of EPT on the basis that it does not include dependence on Da, we estimate average long-term changes in warm season PE ranging from 0.14 to 0.45 mm·day−1. Hence, there is still large uncertainty in how much evaporative demand increased. Greater increases in ERFu and ERFc, owing to strong dependence on increasing Da (and the general pattern of increasing trends with increasing ga and gs) must impose substantial impacts on the level of water stress imposed upon plants, either through a preemptive regulation of gas exchange to avoid desiccation at the expense of carbon supply, or lax regulation of gas exchange at the expense of water transport efficiency and risk of hydraulic failure. Yet either way, the magnitude of those increases in evaporative demand (i.e., the forcing imposed on guard cells or hydraulic conductivity) is not supported by other models of PE, including a statistical fit to EPan. That is not to say, however, that the more modest changes in evaporative demand implied by EPT cannot impose severe water stress, as the sensitivity of plant health and productivity to changes in EPT, specifically, remains largely unexplored. e f l -t i r f l f aera) s)of mean warm-season daily refe nce evapor ti n over f rest regions of North America(z-axis indic tes change in mm·day−1 from 1951–20 4).While plenty of surface water balance research is performed at a sub-monthly time interval,it remains far more practical to simulate the water balance at a monthly time step in large-scale forestecosystem modelling. Analysis of variance attribution suggested that controls on EPan were scaledependent. On a daily time scale, EPan was significantly related to Rsgd, Ta, Da, and U, whereas ona monthly time scale EPan was only significantly related to Rsgd, Ta, and Da. We interpreted thisdifference to mean that systematic effects of advection (as an energy source over pans) and groundheat flux (as an energy sink) diminish as time scale increases. That is, on any given day, advectionand ground heat flux can enhance or suppress PE, while beyond this temporal scale, it cancels out.The lack of effect of U at the monthly time scale may suggest that accuracy in large-scale analysis ofspatial and temporal variation (including trend analysis) may not suffer from failure to represent wind.Varying sensitivity to Da was a defining difference among the most popular models of PE.The Priestley–Taylor model and the Penman–Monteith model parameterized for conditions abovereference forest canopies reflect estimates with absent Da-sensitivity and strong Da-sensitivity,respectively. Any inferences about change in surface hydrology and consequent socio-economicimplications over time are, therefore, going to strongly depend on which model is applied.The differentiation in long-term change among models was exacerbated by strong apparent positivetrend in historical Da across North America [85]. Interestingly, the different dependence on Da amongmodels also shaped the zonal variation in PE (Figure 7), where estimates of PE derived from referenceforest rates actually exceeded estimates of PE with infinite surface conductance in Penman’s model.Specifically, estimates of ERFu and ERFc start to exceed EPen when Da > 18.5 hPa.On average across forests of North America, results indicate that PE increased from 1951–2014.If this is correct, increasing PE must increase the tendency for drought stress. Confidence in the claim ofa positive trend is gained from the fact that all models consistently indicated positive trends. However,the claim of positive trends remains uncertain due to reliance on irradiance from NCEP/NCARreanalysis, which has poorly constrained uncertainty, and by the inability to account for trends inwind speed. These sources of uncertainty will affect all models similarly. If we ignore uncertainty inirradiance and wind trends, it must still be recognized that the range of trend estimates among modelsof PE includes low estimates, which may have been exceeded by increasing trends in precipitation [86].If we exclude consideration of EPT on the basis that it does not include dependence on Da, we estimateHydrology 2017, 4, 5 20 of 24average long-term changes in warm season PE ranging from 0.14 to 0.45 mm·day−1. Hence, there isstill large uncertainty in how much evaporative demand increased.Greater increases in ERFu and ERFc, owing to strong dependence on increasing Da (and the generalpattern of increasing trends with increasing ga and gs) must impose substantial impacts on the level ofwater stress imposed upon plants, either through a preemptive regulation of gas exchange to avoiddesiccation at the expense of carbon supply, or lax regulation of gas exchange at the expense of watertransport efficiency and risk of hydraulic failure. Yet either way, the magnitude of those increasesin evaporative demand (i.e., the forcing imposed on guard cells or hydraulic conductivity) is notsupported by other models of PE, including a statistical fit to EPan. That is not to say, however, that themore modest changes in evaporative demand implied by EPT cannot impose severe water stress, as thesensitivity of plant health and productivity to changes in EPT, specifically, remains largely unexplored.At the regional scale, all trends in PE were positive with the exception of trends in the easterntemperate ecozone, where EPT was negative and EPen was zero. Interestingly, this could be thesignature of high tropospheric aerosol forcing on surface climates over the east-central conterminousU.S. (Figure 8). Similar patterns exist for pockets of industrial activity across the continent. Furtherresearch is needed to understand the quantitative relationship between PE and tropospheric aerosols.If there is a link, it stands to reason that continued trends in PE may be strongly affected by the onsetof aerosol emission regulations that began late in our study period (i.e., 1990’s), and conversely by anynew or growing emission sources.5. ConclusionsWe investigated spatial and temporal variation of PE across forest regions of North America.While records of EPan were not long enough to perform direct trend analysis, we instead choseto produce a statistical model, EPanFit, that helped to place the magnitude of model predictions incontext and evaluate differences in the relative effect of the driving meteorological factors. Consistentwith previous model inter-comparison studies [1,87,88], we found substantial differences in PEamong models. Large-scale (zonal) differences were clearly demonstrated across North America.Long-term change in PE from 1951–2014 ranged from 0.08 to 0.45 mm·day−1 depending on whetherit was calculated with the Priestley–Taylor model, or with the Penman–Monteith model set for areference forest by making the common assumption that the canopy is well coupled to the atmosphere.A general analysis of variance attribution among the approaches helped to put the underlying reasonsfor divergent behavior in quantitative terms. We conclude that until a consensus emerges on whichexpression is more representative of reality in forests, studies that adopt such basic expressions of PEshould consider the range of different model predictions to represent uncertainty in trend magnitude.The need to forecast impacts of a changing climate forests, including wildfire risk, insect outbreakrisk, drought-related mortality, and health decline syndromes, indicates a need to improve monitoringof the key drivers of PE, including Rsgd and Da. This source of high resolution, continuous monthlyvariation of PE from various models can facilitate research on a wide range of applications that requireknowledge of the surface water balance.Acknowledgments: This research was supported by the Pacific Institute for Climate Solutions.Author Contributions: R.H., N.C., and D.S. designed the study. R.H. produced and analyzed results and wrotethe first draft. R.H., N.C. and D.S. provided editorial comments and suggestions.Conflicts of Interest: The authors declare no conflict of interest.References1. Fisher, J.B.; Whittaker, R.J.; Malhi, Y. ET come home: Potential evapotranspiration in geographical ecology.Glob. Ecol. Biogeogr. 2011, 20, 1–18. [CrossRef]2. Webb, W.L.; Lauenroth, W.K.; Szarek, S.R.; Kinerson, R.S. Primary production and abiotic controls in forests,grasslands, and desert ecosystems in the united states. Ecology 1983, 64, 134–151. [CrossRef]Hydrology 2017, 4, 5 21 of 243. Stephenson, N.L. Actual evapotranspiration and deficit: biologically meaningful correlates of vegetationdistribution across spatial scales. J. Biogeogr. 1998, 25, 855–870. [CrossRef]4. Sperry, J.S.; Tyree, M.T. Water-stress-induced xylem embolism in three species of conifers. Plant Cell Environ.1990, 13, 427–436. [CrossRef]5. Tyree, M.T.; Sperry, J.S. Do woody plants operate near the point of catastrophic xylem dysfunction caused bydynamic water stress? Answers from a model. Plant Physiol. 1988, 88, 574–580. [CrossRef] [PubMed]6. Choat, B.; Jansen, S.; Brodribb, T.J.; Cochard, H.; Delzon, S.; Bhaskar, R.; Bucci, S.J.; Feild, T.S.; Gleason, S.M.;Hacke, U.G.; et al. Global convergence in the vulnerability of forests to drought. Nature 2012, 491, 752–755.[CrossRef] [PubMed]7. Hember, R.A.; Kurz, W.A.; Coops, N.C. Relationships between individual-tree mortality and water-balancevariables indicate positive trends in water stress-induced tree mortality across North America. Glob. ChangeBiol. 2016. [CrossRef] [PubMed]8. Millar, C.I.; Stephenson, N.L. Temperate forest health in an era of emerging megadisturbance. Science 2015,349, 823–826. [CrossRef] [PubMed]9. McDowell, N.G.; Williams, A.P.; Xu, C.; Pockman, W.T.; Dickman, L.T.; Sevanto, S.; Pangle, R.; Limousin, J.;Plaut, J.; Mackay, D.S.; et al. Multi-scale predictions of massive conifer mortality due to chronic temperaturerise. Nat. Clim. Change 2016, 6, 295–300. [CrossRef]10. Martínez-Vilalta, J.; Lloret, F. Drought-induced vegetation shifts in terrestrial ecosystems: The key role ofregeneration dynamics. Glob. Planet. Change 2016, 144, 94–108. [CrossRef]11. Clark, J.S.; Iverson, L.; Woodall, C.W.; Allen, C.D.; Bell, D.M.; Bragg, D.C.; D’Amato, A.W.; Davis, F.W.;Hersh, M.H.; Ibanez, I.; et al. The impacts of increasing drought on forest dynamics, structure, andbiodiversity in the United States. Glob. Change Biol. 2016. [CrossRef] [PubMed]12. Sucoff, E. Water potential in red pine: Soil moisture, evapotranspiration, crown position. Ecology 1972, 53,681–686. [CrossRef]13. Sharkey, T.D. Transpiration-induced changes in the photosynthetic capacity of leaves. Planta 1984, 160,143–150. [CrossRef] [PubMed]14. Novick, K.A.; Ficklin, D.L.; Stoy, P.C.; Williams, C.A.; Bohrer, G.; Oishi, A.C.; Papuga, S.A.; Blanken, P.D.;Noormets, A.; Sulman, B.N.; et al. The increasing importance of atmospheric demand for ecosystem waterand carbon fluxes. Nat. Clim. Change 2016. [CrossRef]15. Sanderson, M. Drought in the Canadian Northwest. Geogr. Rev. 1948, 38, 289. [CrossRef]16. Thornthwaite, C.W. An approach toward a rational classification of climate. Geogr. Rev. 1948, 38, 55–94.[CrossRef]17. McMahon, T.A.; Peel, M.C.; Lowe, L.; Srikanthan, R.; McVicar, T.R. Estimating actual, potential, referencecrop and pan evaporation using standard meteorological data: A pragmatic synthesis. Hydrol. EarthSyst. Sci. 2013, 17, 1331–1363. [CrossRef]18. Oke, T.R. Boundary Layer Climates, 2nd ed.; Routledge: Cambridge, UK, 1987.19. Rosenberg, N.; Blad, B.; Verma, S. Microclimate, the Biological Environment, 2nd ed.; Jon Wiley and Sons Inc.:Toronto, ON, Canada, 1983.20. Slatyer, R.; McIlroy, I. Practical Microclimatology. Commonw. Sci. Ind. Res. Organ. 1961. [CrossRef]21. Mather, J.R.; Ambroziak, R.A. A Search for Understanding Potential Evapotranspiration. Geogr. Rev. 1986,76, 355–370. [CrossRef]22. Ferguson, H.L.; O’Neill, A.D.J.; Cork, H.F. Mean Evaporation over Canada. Water Resour. Res. 1970, 6,1618–1633. [CrossRef]23. Kohler, M.A.; Nordenson, T.J.; Fox, W.E. Evaporation from Pans and Lakes; US Department of Commerce:Washington, DC, USA, 1955.24. Thornthwaite, C.W.; Mather, J.R. The Water Balance; Laboratory of Climatology: Centerton, NJ, USA, 1955.25. Hargreaves, G.H.; Allen, R.G. History and Evaluation of Hargreaves Evapotranspiration Equation. J. Irrig.Drain. Eng. 2003, 129, 53–63. [CrossRef]26. Christiansen, J.E. Estimating Evaporation and Evpoatranspiration from Climate Data; Utah Water ResearchLaboratory, Utah State University: Logan, UT, USA, 1966.27. Penman, H.L. Natural evaporation from open water, bare soil and grass. Proc. R. Soc. Lond. 1948, 193,120–145. [CrossRef]28. Monteith, J.L. Evaporation and environment. Symp. Soc. Exp. Biol. 1965, 19, 205–234. [PubMed]Hydrology 2017, 4, 5 22 of 2429. Priestley, C.H.B.; Taylor, R.J. On the assessment of surface heat flux and evaporation using large-scaleparameters. Mon. Weather Rev. 1972, 100, 81–92. [CrossRef]30. Black, T.A. Evapotranspiration from douglas-fir stands exposed to soil-water deficits. Water Resour. Res.1979, 15, 164–170. [CrossRef]31. Brümmer, C.; Black, T.A.; Jassal, R.S.; Grant, N.J.; Spittlehouse, D.L.; Chen, B.; Nesic, Z.; Amiro, B.D.;Arain, M.A.; Barr, A.G.; et al. How climate and vegetation type influence evapotranspiration and wateruse efficiency in Canadian forest, peatland and grassland ecosystems. Agric. For. Meteorol. 2012, 153, 14–30.[CrossRef]32. Greve, P.; Seneviratne, S.I. Assessment of future changes in water availability and aridity. Geophys. Res. Lett.2015, 42. [CrossRef] [PubMed]33. Hogg, E.H. Temporal scaling of moisture and the forest-grassland boundary in western Canada. Agric. For.Meteorol. 1997, 84, 115–122. [CrossRef]34. Willmott, C.J.; Rowe, C.M.; Mintz, Y. Climatology of the terrestrial seasonal water cycle. J. Climatol. 1985, 5,589–606. [CrossRef]35. Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop Evapotranspiration: Guidelines for Computing Crop WaterRequirements, FAO Irrigation and Drainage Paper 56; Food and Agriculture Organization of the United Nations:Rome, Italy, 1998.36. Cleugh, H.A.; Leuning, R.; Mu, Q.Z.; Running, S.W. Regional evaporation estimates from flux tower andMODIS satellite data. Remote Sens. Environ. 2007, 106, 285–304. [CrossRef]37. Spittlehouse, D.L. Water availability, climate change and the growth of Douglas-fir in the Georgia Basin.Can. Water Resour. J. 2003, 28, 673–687. [CrossRef]38. Coops, N.; Coggins, S.; Kurz, W. Mapping the environmental limitations to growth of coastal Douglas-firstands on Vancouver Island, British Columbia. Tree Physiol. 2007, 27, 805–815. [CrossRef] [PubMed]39. Seely, B.; Welham, C.; Scoullar, K. Application of a hybrid forest growth model to evaluate climate changeimpacts on productivity, nutrient cycling and mortality in a montane forest ecosystem. PLoS ONE 2015, 10,e0135034. [CrossRef] [PubMed]40. King, D.A.; Bachelet, D.M.; Symstad, A.J.; Ferschweiler, K.; Hobbins, M. Estimation of potentialevapotranspiration from extraterrestrial radiation, air temperature and humidity to assess future climatechange effects on the vegetation of the Northern Great Plains, USA. Ecol. Model. 2015, 297, 86–97. [CrossRef]41. Hogg, E.H.; Barr, A.G.; Black, T.A. A simple soil moisture index for representing multi-year drought impactson aspen productivity in the western Canadian interior. Agric. For. Meteorol. 2013, 178–179, 173–182.[CrossRef]42. Wang, Y.; Hogg, E.H.; Price, D.T.; Edwards, J.; Williamson, T. Past and projected future changes in moistureconditions in the Canadian boreal forest. For. Chron. 2014, 90, 678–691. [CrossRef]43. Chattopadhyay, N.; Hulme, M. Evaporation and potential evapotranspiration in India under conditions ofrecent and future climate change. Agric. For. Meteorol. 1997, 87, 55–73. [CrossRef]44. Hobbins, M.; Ramirez, J.; Brown, T. Trends in pan evaporation and actual evapotranspiration across theconterminous US: Paradoxical or complementary? Geophys. Res. Lett. 2004, 31. [CrossRef]45. Linacre, E.T. Estimating U.S. Class A pan evaporation from few climate data. Water Int. 1994, 19, 5–14.[CrossRef]46. Hobbins, M.T.; Dai, A.; Roderick, M.L.; Farquhar, G.D. Revisiting the parameterization of potentialevaporation as a driver of long-term water balance trends. Geophys. Res. Lett. 2008, 35, L12403. [CrossRef]47. Kingston, D.G.; Todd, M.C.; Taylor, R.G.; Thompson, J.R.; Arnell, N.W. Uncertainty in the estimation ofpotential evapotranspiration under climate change. Geophys. Res. Lett. 2009, 36, L20403. [CrossRef]48. Scheff, J.; Frierson, D.M.W. Scaling potential evapotranspiration with greenhouse warming. J. Clim. 2013, 27,1539–1558. [CrossRef]49. Dai, A. Increasing drought under global warming in observations and models. Nat. Clim. Change 2013, 3,52–58. [CrossRef]50. Dai, A. Characteristics and trends in various forms of the Palmer Drought Severity Index during 1900–2008.J. Geophys. Res. Atmospheres 2011, 116. [CrossRef]51. Sheffield, J.; Wood, E.F.; Roderick, M.L. Little change in global drought over the past 60 years. Nature 2012,491, 435–438. [CrossRef] [PubMed]Hydrology 2017, 4, 5 23 of 2452. Amatya, D.; Tian, S.; Dai, Z.; Sun, G. Long-term potential and actual evapotranspiration of two differentforests on the Atlantic Coastal Plain. Trans. ASABE 2016, 59, 647–660.53. Brauman, K.A.; Freyberg, D.L.; Daily, G.C. Potential evapotranspiration from forest and pasture in thetropics: A case study in Kona, Hawai‘i. J. Hydrol. 2012, 440–441, 52–61. [CrossRef]54. NREL 2013 US Department of Energy, National Renewable Energy Laboratory. National Solar RadiationData Base. Available online: http://rredc.nrel.gov/solar/old_data/nsrdb/ (accessed on 26 April 2016).55. GSOD 2013 National Oceanic and Atmospheric Administration, Global Summary of the Day.Available online: https://data.noaa.gov/dataset/global-surface-summary-of-the-day-gsod (accessed on31 October 2013).56. Brouwer, C.; Heibloem, M. Irrigation Water Management: Irrigation Water Needs; Food and AgricultureOrgainzation of the United Nations: Rome, Italy, 1985.57. Lu, J.; Sun, G.; McNulty, S.G.; Amatya, D. A comparison of six potential evapotranspiration methods forregional use in the Southeastern United States. J. Am. Water Resour. Assoc. 2005, 41, 621–633. [CrossRef]58. Kelliher, F.; Leuning, R.; Schulze, E. Evaporation and canopy characteristics of coniferous forests andgrasslands. Oecologia 1993, 95, 153–163. [CrossRef]59. Landsberg, J.; Sands, P. Physiological Ecology of Forest Production: Principles, Processes and Models, 1st ed.;Elsevier Inc.: Boston, MA, USA, 2011.60. Pereira, A.R. The Priestley–Taylor parameter and the decoupling factor for estimating referenceevapotranspiration. Agric. For. Meteorol. 2004, 125, 305–313. [CrossRef]61. Bladon, K.D.; Silins, U.; Landhäusser, S.M.; Lieffers, V.J. Differential transpiration by three boreal tree speciesin response to increased evaporative demand after variable retention harvesting. Agric. For. Meteorol. 2006,138, 104–119. [CrossRef]62. Humphreys, E.; Black, T.; Ethier, G.; Drewitt, G.; Spittlehouse, D.; Jork, E.; Nesic, Z.; Livingston, N. Annualand seasonal variability of sensible and latent heat fluxes above a coastal Douglas-fir forest, British Columbia,Canada. Agric. For. Meteorol. 2003, 115, 109–125. [CrossRef]63. Martin, T.A.; Hinckley, T.M.; Meinzer, F.C.; Sprugel, D.G. Boundary layer conductance, leaf temperature andtranspiration of Abies amabilis branches. Tree Physiol. 1999, 19, 435–443. [CrossRef] [PubMed]64. Jarvis, P.G.; McNaughton, K.G. Stomatal Control of Transpiration: Scaling Up from Leaf to Region.In Advances in Ecological Research; Elsevier: Amsterdam, The Netherlands, 1986; Volume 15, pp. 1–49.65. Running, S.; Coughlan, J. A General-Model of Forest Ecosystem Processes for Regional Applications. 1.Hydrologic Balance, Canopy Gas-Exchange and Primary Production Processes. Ecol. Model. 1988, 42,125–154. [CrossRef]66. Bates, D.M. Lme4: Mixed-Effects Modeling with R; Springer: New York, NY, USA, 2010.67. Kalnay, E.; Kanamitsu, M.; Kistler, R.; Collins, W.; Deaven, D.; Gandin, L.; Iredell, M.; Saha, S.; White, G.;Woollen, J.; et al. The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteorol. Soc. 1996, 77, 437–471.[CrossRef]68. Schroeder, T.; Hember, R.; Coops, N.; Liang, S. Validation of solar radiation surfaces from modis andreanalysis data over topographically complex terrain. J. Appl. Meteorol. Climatol. 2009, 48, 2441–2458.[CrossRef]69. NARR 2013 National Oceanic and Atmospheric Administration, Earth System Research Laboratory. NorthAmerican Regional Reanalysis. Available online: www.esrl.noaa.gov/psd/data/gridded (accessed on18 December 2013).70. ESRI. ArcGIS 2013 ArcGIS 9.2 Desktop Help. Available online: webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Area_Solar_Radiation (accessed on 23 July 2013).71. Cogley, J.G. Greenland accumulation: An error model. J. Geophys. Res. Atmospheres 2004, 109. [CrossRef]72. Press, W.; Teukolsky, S.; Vetterling, W.; Flannery, B. Numerical Recipes: The Art of Scientific Computing, 3rd ed.;Cambridge Universiy Press: New York, NY, USA, 2007.73. Wang, T.; Hamann, A.; Spittlehouse, D.; Aitken, S. Development of scale-free climate data for westernCanada for use in resource management. Int. J. Climatol. 2006, 26, 383–397. [CrossRef]74. Wang, T.; Hamann, A.; Spittlehouse, D.; Carroll, C. Locally downscaled and spatially customizable climatedata for historical and future periods for north america. PLoS ONE 2016, 11, e0156720. [CrossRef] [PubMed]75. EC-NCDIA 2013 Environment Canada, National Climate Data and Information Archive. 2013. Availableonline: http://climate.weather.gc.ca/index_e.html (accessed on 5 January 2013).Hydrology 2017, 4, 5 24 of 2476. Menne, M.J.; Williams, C.N.; Vose, R.S. United States Historical Climatology Network (USHCN); Version 2.5Serial Monthly Dataset; Oak Ridge National Laboratory, Carbon Dioxide Information Analysis Center:Oak Ridge, TN, USA.77. EC-AHCCD 2013 Environment Canada, Adjusted and Homogenized Canadian Climate Data. 2013. Availableonline: http://www.ec.gc.ca/dccha-ahccd/ (accessed on 15 July 2013).78. CRU31 2013 UK Met Office, Climate Research Unit. Climate Database, Version 3.1; Available online: http://badc.nerc.ac.uk/view/badc.nerc.ac.uk__ATOM__dataent_1256223773328276 (accessed on 25 July 2013).79. Kimball, J.; Running, S.; Nemani, R. An improved method for estimating surface humidity from dailyminimum temperature. Agric. For. Meteorol. 1997, 85, 87–98. [CrossRef]80. Thornton, P.; Running, S.; White, M. Generating surfaces of daily meteorological variables over large regionsof complex terrain. J. Hydrol. 1997, 190, 214–251. [CrossRef]81. USGS National Atlas of the United States. Available online: https://nationalmap.gov/small_scale/ (accessedon 13 March 2015).82. Coakley, J.A. Reflectance and albedo, surface. In Encyclopedia of the Atmospheric Sciences; Holton, J.R.,Curry, J.A., Pyle, J.A., Eds.; Academic Press: Cambridge, MA, USA, 2002; pp. 1914–1923.83. Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC). FLUXNET Web Page;ORNL DAAC: Oak Ridge, TN, USA, 2015.84. CEC 1997 Commission for Environmental Cooperation (CEC): Ecological Regions of NorthAmerica—Towards a Common Perspective. GIS Data (shapefiles, metadata and symbology). Availableonline: http://www.cec.org/tools-and-resources/map-files/terrestrial-ecoregions-level-i (accessed on29 November 2009).85. Isaac, V.; van Wijngaarden, W.A. Surface water vapor pressure and temperature trends in North Americaduring 1948–2010. J. Clim. 2012, 25, 3599–3609. [CrossRef]86. Vincent, L.A.; Zhang, X.; Brown, R.D.; Feng, Y.; Mekis, E.; Milewska, E.J.; Wan, H.; Wang, X.L. ObservedTrends in Canada’s Climate and Influence of Low-Frequency Variability Modes. J. Clim. 2015, 28, 4545–4560.[CrossRef]87. Vörösmarty, C.J.; Federer, C.A.; Schloss, A.L. Potential evaporation functions compared on US watersheds:Possible implications for global-scale water balance and terrestrial ecosystem modeling. J. Hydrol. 1998, 207,147–169. [CrossRef]88. Yao, H. Long-term study of lake evaporation and evaluation of seven estimation methods: Results fromDickie Lake, south-central Ontario, Canada. J. Water Resour. Prot. 2009, 1, 59–77. [CrossRef]© 2017 by the authors; licensee MDPI, Basel, Switzerland. This article is an open accessarticle distributed under the terms and conditions of the Creative Commons Attribution(CC-BY) license (http://creativecommons.org/licenses/by/4.0/).
Article
eng
Reviewed
Vancouver : University of British Columbia Library
Multidisciplinary Digital Publishing Institute
CC BY 4.0
Faculty
Potential evaporation
Evaporative demand
Pan evaporation
Climate trends
Spatial and Temporal Variability of Potential Evaporation across North American Forests
Text