Spatial and Temporal Variability of Potential Evaporation across North American Forests

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.


Introduction
The evaporating power of the atmosphere, or simply evaporative demand, exerts a primary control on the distribution, abundance, and diversity of life on land [1][2][3].On one hand, the atmospheric demand for water allows plants to draw water above the surface without exerting vast amounts of energy.On the other, fluctuations in evaporative demand simultaneously pose a risk to plants if soil water supply and transport efficiency cannot meet peak evaporative demand [4][5][6].These compensating factors appear to put forest biomes at risk under anticipated climate change [7][8][9][10][11].While declining precipitation must impose major constraints on productivity of trees across the forest-grassland ecotone, there is a growing appreciation that evaporative demand also strongly constrains productivity of humid forest biomes [12][13][14].For these reasons, evaporative demand is an important 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 the rate of evaporation that would occur from a surface with unlimited water supply and no resistance to transfer of water.It is, therefore, determined by the energy available to change the phase of water from liquid to gas [17][18][19][20].Although PE has a conceptual definition, it is not easy to measure directly above land 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 at many meteorological stations across Canada and the U.S. and has been used to assess climatological variability of PE [22,23].
Despite extensive investigation of evaporation from pans, lakes, and lysimeters, empirical studies have not resulted in a universal approach to predict PE.Some of the earliest models were based on temperature and were developed with the limited availability of radiation, humidity, and wind measurements in mind [16,24,25].Other models were simultaneously developed to predict PE based on a broader number of driving variables, including net radiation (R n ), temperature (T a ), vapour pressure deficit (D a ), and wind speed (U) [26,27].Later, Penman's model was modified to produce the Penman-Monteith equation, which accounts for varying surface conductance to water transfer [28].While the Penman-Monteith model incorporates an aerodynamic component, other models were developed that omit an aerodynamic component.For example, it was later argued that at regional scales PE could be approximated by multiplying the equilibrium rate of evaporation by a constant value, α = 1.26 [29].Although empirical studies suggest that the overall magnitude of PE using this Priestley-Taylor coefficient of α = 1.26 will bring estimates within closer proximity to models that include an aerodynamic component, spatial and temporal variability in estimates will only reflect variation in T a and R n .In forests, α must be modified [30,31], yet it is still a useful concept if it is recognized that there is no universal value of α.Lastly, other studies express PE based solely on a function of R n [32], or solely based on D a [33].
The above modelling approaches all continue to be used broadly in hydrological and ecological research.For example, temperature-based models are applied in relatively basic surface water balance models [34].Use of Penman's model has also been applied in ecological studies [2].The Penman-Monteith model is widely applied [35][36][37][38].Other studies represent potential evapotranspiraiton based on the Priestley-Taylor model [30,39], or close equivalents based on R n and T a [25,40], or R n alone [32], or D a alone [41,42].Despite its importance in modelling the surface water balance and subsequent application in ecological research, differences in the variation of evaporative demand resulting from these approaches, including long-term trends at the regional scale, are not well documented.
Historical and projected changes in PE may arise from natural climate variability, and changing atmospheric greenhouse gas and aerosol concentrations.However, the response of PE to regional climate change involves scale-and surface-dependent relationships with multiple meteorological drivers that are not always directly measured.As such, continuous high temporal and spatial resolution maps of PE and projections of future change required to study many ecological processes remain relatively scarce.For regions with long-term continuous measurements of pan evaporation, some studies report a decrease in pan evaporation [43][44][45].Other studies resort to use of models to understand trends.For example, Penman's model suggested a decrease in historical PE across some regions of India [43], but also projected long-term future increases of PE.Across Australia, trends in PE inferred 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 and the Penman-Monteith model, concluded that differences among models partially led to differences in the direction of future trends in aridity index [47].Trends in PE based on the Penman-Monteith equation generally indicate increasing trends in evaporative demand [7,[48][49][50].Positive future trends in PE are largely attributed to increasing saturation vapour pressure that accompanies global warming [48].Comparison of historical trends in drought using Penman's model with previous studies that applied Thornthwaite's temperature-based model of PE emphasized that the latter might exaggerate increasing trends in aridity [51].Conversely, time trends in global land aridity of Coupled Model Intercomparison Project Phase 5 (CMIP5) model simulations based on estimates of PE inferred from just net radiation were commonly not significant [32].
There is considerable demand for a spatial and temporal index of change in evaporative demand to 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 America using six different model-based estimates of PE.To put the magnitude and underlying mechanics of model-based estimates of PE into context, we compared predictions against measurements of pan evaporation (E Pan ) for 28 sites across Canada and the conterminous U.S..One of the six estimates of PE was defined by a statistical model fit to E Pan measurements, while the other five estimates reflected approaches that are widely used in applied research, but that require no calibration.First, we evaluated the agreement between models and E Pan .Second, we evaluated the variance in PE that is attributed to each meteorological driver in each of the models.Lastly, to understand spatial and temporal variation in the model estimates, we developed spatially-continuous monthly mean values of the driving variables covering North America from 1951-2014.Maps were then generated to document long-term averages ('climatologies') of PE, and historical time trends in PE from 1951-2014.

Pan Evaporation Measurements
Daily total pan evaporation (E Pan ) measurements from Class A pans were collected across North America.Data collected at stations in Canada were requested and received from Environment Canada during March 2014.Data collected at stations in the U.S. were downloaded from the National Oceanic and Atmospheric Administration during March 2016.Any stations with available daily measurements between 1950 and 2014 were considered.
In Canada, large periods of the pan evaporation measurements were accompanied by simultaneous measurements of pan water temperature, air temperature (T a ), and wind speed (U).Separate databases of daily mean downward global shortwave radiation (R sgd ) and dew-point temperature were also available.Daily mean vapour pressure deficit of the air (D a ) was calculated from air temperature and dew-point temperature.From the initial databases, ten stations in Canada had multiple years with overlapping measurements of E Pan and meteorological variables (all of R sgd , T a , D a , and U).For all U.S. stations, daily mean values of R sgd were downloaded from two National Renewable Energy Laboratory (NREL) databases [54].All other meteorological variables were downloaded from the Global Summary of the Day (GSOD) database [55].Relatively few U.S. stations simultaneously measured pan evaporation and meteorological drivers.By manually cross-referencing between databases, we identified eighteen stations that collected pan evaporation measurements and were also within close proximity to stations with the necessary meteorological measurements (a distance of less than 50 km in most cases).Daily values of R sgd , T a , D a , and U for each station with pan evaporation measurements were then estimated from inverse-distance weighting interpolation of measurements from the ten nearest stations.As a quality control measure, multivariate regression was used to relate pan evaporation against meteorological variables at each candidate station.Despite reliance on off-site meteorological variables in the U.S., we did not observe a significant difference in the skill or behaviour of the models between Canada and the U.S. (not shown).Combining the stations from Canada and U.S. produced a total sample of 28 stations, each consisting of multiple years of measurements for analysis (Table 1).Stations were relatively evenly distributed across North America, with the exception of Mexico (Figure 1).Although the emphasis was on documenting variability of PE for forest biomes of North America, we chose to analyze all 28 available stations.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.

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: where s is the slope of the saturation vapour pressure-over-temperature relationship for air (kPa , 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]: A pan coefficient, k pan , 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 k pan = 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 k pan = 0.77 [45], while irrigation manuals report a value of k pan = 0.80 for moderate wind and moderate humidity [56].In this study, we applied the mid-range value, k pan = 0.77.

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: where s is the slope of the saturation vapour pressure-over-temperature relationship for air (kPa ), e s is saturation vapour pressure (kPa), e a 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]: where G, s, γ, and λ L are as defined above.The Priestley-Taylor model of potential evaporation was then calculated from [29]: Above wet surfaces, they argued that α = 1.26 yields a strong predictor of potential evaporation at the regional scale.
A 'reference' rate of evaporation (E Ref ) is commonly defined for a specific type of surface vegetation.
It is generally defined as an intermediate variable between PE and actual evapotranspiration.In essence, it indicates the maximum rate of actual evapotranspiration [19], or the expected rate of actual evapotranspiration at maximum hydraulic conductivity of plants and soil.For example, the reference evapotranspiration for a hypothetical surface of vegetation with a characteristic height and surface conductance can be calculated from the Penman-Monteith equation as [35]: where ρ a is the density of air (kg•m −3 ), and c p is the specific heat of the air (J•kg −1 •K −1 ), g a is aerodynamic conductance (m•s −1 ), and g s is surface conductance (m•s −1 ).In reality, surface conductance is highly variable across forest types and environments, yet here we demonstrate the behaviour of predictions from the Penman-Monteith equation when g a and g s are set as constants.
For example, for a reference grassland, conductance values are commonly set to g a = 0.010 m•s −1 and g s = 0. 014 m•s −1 .In so doing, the reference evapotranspiration rate specifically for grasslands (E RG ) 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 representative of forest canopies [58,59].We, therefore, considered application of the Penman-Monteith model to simulate reference values for forest canopies.Variation in forest canopy structure affects the conversion of sensible heat from the surrounding air into latent heat [60][61][62][63].As surface roughness increases, so does the potential extraction of sensible heat and the surface is said to be increasingly 'coupled' to the atmosphere.Conversely, as surface roughness decreases, less sensible heat can be extracted from the 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 to the atmosphere (E RFu ) and those that are well coupled to the atmosphere (E RFc ).For canopies that are moderately coupled to the atmosphere, we set aerodynamic conductance as g a = 0.058 m•s −1 , consistent with average conditions recently reported for coniferous, deciduous, and mixed-wood forests using eddy-covariance systems [31].For canopies that are well coupled to the atmosphere, we set aerodynamic conductance as g a = 0.150 m•s −1 , more consistent with the value commonly assumed in models [59,65].In choosing to focus on these two levels of aerodynamic conductance for forest canopies, it was not intended to advocate either one, but rather to understand the implications of each level.In both estimates for reference forests, we applied an equivalent value of surface conductance, g s = 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 fit against daily measurements of E Pan .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 assumed that daily values of E Pan would be related to four variables, including downward global shortwave radiation, air temperature, vapour pressure deficit, and wind speed, yielding: where β 0 ...β 4 are fixed effects and b 0,i ...b 4,i , are random effects applied to i = 1...n stations.To the extent that the assumed drivers and the assumption of additive linear effects can explain PE, this estimate also reflects the upper limit of how well models should perform.In all cases, the response variable, E Pan , and the predictor variables, R sgd , T a , D a , and U, were standardized (i.e., 'z-scored') by subtracting the sample mean and dividing by the sample standard deviation.The resulting coefficients express the magnitude of change in the response variable (in terms of number of standard deviations) that is imposed by an increase of one standard deviation in each predictor variable.All models were fitted in R (R Development Core Team, 2015) using the lmer function from lme4 package [66].As interest in E PanFit was to describe the behavior of the calibrated model rather than as a predictive tool, we evaluated 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 .Unless otherwise specified, we focused on reporting mean daily PE over a consistent period during the warm season when most forest regions of North America are above freezing, defined as the months from May to September.

Comparison between Models and Pan Evaporation
Comparison between models and pan evaporation was first performed on a daily time scale by regressing model predictions of daily total PE with daily total E Pan .For each comparison, we produced scatterplots and reported the coefficient of determination (R 2 ), 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 predict monthly mean PE has broader practical value in applied research.Monthly mean values of daily total pan evaporation were calculated from the daily values, only considering months with no missing daily 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 evaporation were plotted against individual monthly mean values for each of the driving variables.Scatterplots and performance statistics were reported for predictions of monthly mean potential evaporation based on E Pen , E PT , E RG , E RFu , E RFc , and E PanFit .

Analysis of Variance Attribution
The relative importance of the drivers was inferred from estimates of sensitivity assessed first on a daily time scale.For the fitted model, the sensitivity of drivers was simply inferred from the fixed effects coefficients of the LME model.For example, the sensitivity of daily E Pan to daily R sgd was 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 E Pen and E PT , E RG , E RFu , and E RFc .By comparing coefficients, we were able to infer whether attribution in the models was consistent with that 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 E Pan to produce a second (monthly) version of E PanFit , as well as sets of coefficients for monthly mean E Pen , E PT , E RG , E RFu , and E RFc .The coefficients were compared to understand any differences in attribution between methods.As the coefficients were developed at daily and monthly time scales, the analysis was able to assess whether attribution was scale-dependent.

Model Application
To apply the models across North America, monthly climate data were assembled for a 1 km grid covering North America between 1951 and 2014.Where possible, the database was developed by compiling existing datasets.Other variables, developed from raw data sources, are described in more detail below.For each variable, processing and storage of grids was performed separately for long-term monthly normals and time series of monthly anomalies for 1951-2014.The former reflect twelve grids indicating the mean monthly values for January to December for the base period 1971-2000 for each variable.The latter (anomalies) indicate the monthly deviations from the 1971-2000 normal for all months during 1951-2014.
Transient monthly variation in R sgd was taken from the National Centers for Environmental Predition/National Center for Atmosphere Research (NCEP/NCAR) global reanalysis [67].Similar to the methodology we applied previously [68], a first approximation of monthly normal R sgd was developed by combining the synoptic and zonal variation captured by values from the North American Regional Reanalysis [69] with topographically-driven variation captured by values from the ArcGIS Solar Radiation Tool [70].In the first step, estimates of R sgd from North American Regional Reanalysis (NARR) were downloaded and resampled from the original resolution of 32 km down to a resolution of 1 km using cubic interpolation.In the second step, potential R sgd was calculated using the ArcGIS Solar Radiation Tool with default settings (eight zenith divisions, eight azimuth divisions, uniform sky diffuse model type, a diffuse proportion of 0.3 and a transmittivity of 0.5) for each month at 1 km resolution.The tool was run with a digital elevation model developed by the U.S. Geological Survey the HYDRO1k project (https://lta.cr.usgs.gov/HYDRO1K).The ArcGIS Solar Radiation Tool does not take latitude as an input argument so the function was executed for 0.5 • -wide zonal bands with 0.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 R sgd by correcting the ArcGIS Solar Radiation Tool estimates so that they were equivalent to the NARR estimates at the original scale of the NARR dataset (i.e., 32 km).This was achieved by applying a two-dimensional 32-km average filter ("fspecial" and "filter2" functions, The Mathworks Matlab 7.9.0) to the difference between NARR and ArcGIS estimates and adding it to the estimates of potential radiation from ArcGIS.The intended outcome was an unbiased estimate of the initial NARR dataset with added subgrid-scale variation caused by topographic effects.That is, the ArcGIS Solar Radiation Tool accounts for expected variation in irradiance caused by earth-sun geometry and shading by neighbourhood terrain, while NARR estimates express the synoptic-scale variation.
To expand on previous accuracy assessment [68], we compared the first approximation with measurements of monthly R sgd across North America provided by Environment Canada and from the National Solar Radiation Data Base for the United States produced by the National Renewable Energy Laboratory [54].The comparison indicated that initial NARR estimates systematically overestimated R sgd by 15%-30% depending on month of the year.We made the assumption that the estimates from station networks were closer to the true values and developed a second approximation of R sgd that corrected the first approximation to better match the station values.This was achieved in two steps: First, the initial monthly biases between NARR and the station observations (an average offset of 1.73 MJ•day −1 ) was subtracted from the first approximation.Second, the remaining residual errors at each station were interpolated to the grid using a polynomial surface fit with spatial variation defined by the errorgram and solved using singular value decomposition [71,72].The resulting continuous grid of differences suggested that NARR dataset overestimated R sgd in the southwest U.S. and large regions of central Canada and underestimated it in the northeast conterminous U.S. The second (and final) approximation of R sgd was then developed by adding the residual error field to the first approximation.
To validate the second approximation of R sgd , nearest grid cell estimates were compared with 25 independent Fluxnet tower observations.Tower observations spanned different years between 1998 and 2010 so a proper validation of the base period was not possible.Deviations from the long-term mean from the NARR dataset, corresponding with the period of comparison with each tower were, therefore, added to the second approximation grid to ensure the comparison was for estimates representative of the same period.Removing regional biases between the first approximation and Environment Canada/NREL station networks increased the level of explained variance from 75% 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 consistently below 10% for each month of year.
A first approximation of monthly normal T a was compiled from the ClimateNA dataset [73,74].Transient variability in monthly temperature was derived from three sources, including Environment Canada [75], the U.S. Historical Climate Station Network (USHCN) [76] for the conterminous U.S. and from the Global Summary of the Day (GSOD) product [55] for Alaska.For Canada, daily climate records were extracted for 697 stations that included all stations listed in Environment Canada's reference climate station network, all stations in the Adjusted and Homogenized Canadian Climate Data archive [77], and some additional remote stations with long records.Canadian station time series 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 the difference.Seven stations were removed in the process due to poor record quality.Gaps in the original records were filled using regression analysis to relate each station record with that of four nearest neighbours.Monthly mean air temperature was then calculated from the gap-filled records of daily minimum and maximum air temperatures.
Transient monthly variation of D a was then derived from resampling of the Climate Research Unit estimates of vapour pressure [78] and saturation vapour pressure (e a ) derived from T a .A first approximation of long-term mean monthly D a was developed based on a function of minimum and maximum temperature [79,80], with the temperature data coming from ClimateNA [74].As with radiation, the first approximation was compared with measurements compiled from the Global Summary of the Day dataset [55].Environment Canada measurements were also added after noticing that not all Environment Canada data appears to have been entered into the GSOD database.The difference between nearest grid cells of the first approximation and stations was then calculated and several outliers were censored from the station record.The differences were then interpolated to produce continuous maps of the error, showing considerable differences over the southern conterminous U.S. Adding the difference field to the first approximation to produce the second approximation led to improved agreement with independent observations of vapour pressure deficit at Fluxnet towers.The level of explained variance in warm-season mean vapour pressure deficit increased 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 preparation of climate input grids included R sgd , yet models of PE required input grids of R n .Net radiation was derived from incident shortwave radiation in two steps.First, we converted incident shortwave radiation to net shortwave radiation to account for differences in mean albedo of various forest surface types.This was achieved by producing a lookup table linking surface albedo values for land cover classes.Land cover classes were derived from the 2002 version of the North American Land Cover Characteristics 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 were drawn from a previous synthesis of values [82].Net shortwave radiation (R sgn ) was then calculated by multiplying R sgd by (1-albedo).In a second step, we developed linear mixed effects models that related monthly mean observations of net shortwave radiation (W•m −2 ) to observations of monthly mean 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 Ridge National 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 of 1054 monthly observations that were representative of deciduous and evergreen forest stands from boreal and temperate forest biomes.The global model yielded a relationship, R n = 0.830R sgn − 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 R 2 = 0.95, and the root mean squared error was 13.5 W•m −2 .We additionally tested whether differences existed between deciduous and evergreen stands, or between boreal and temperate forest biomes.We found that parameters were not statistically different between boreal and temperate stands (p = 0.23), but we did find that parameters differed between deciduous and evergreen stands (p = 0.009).The model for deciduous stands yielded a relationship, R n = 0.882R sgn − 39.500.The slope and intercept had standard errors of 0.011 and 2.300, respectively, and were both significant (p < 0.001).The model for evergreen stands yielded a relationship, R n = 0.817R sgn − 29.512.The slope and intercept had standard errors of 0.031 and 5.894, respectively, and were both significant (p < 0.001).Parameters for deciduous and evergreen forest were linked to those corresponding land cover classes in the USGS land cover classification 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 seasonal cycles in heat capacity, for simplicity, we assumed that the ground heat flux was zero at a decadal time scale.Evaluation of the long-term mean ground heat flux from 12 of the aforementioned Fluxnet sites located in forest ecosystems suggested that this assumption may lead to an overestimate in net radiation of 0.11 percent.
We summarized long-term warm-season mean values of PE from each method by forest region 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 summarized the 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 forest ecozone, 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 stations in the eastern temperate forest ecozone.There is a high degree of uncertainty in values of g a and g s at this scale.We therefore investigated the sensitivity of trends to the parameterization of the Penman-Monteith model by calculating time trends in the resulting values of PE for a broad range in values of g a and g s and then applied multivariate regression analysis to derive a general relationship between time trends in PE and values of g a and g s used in the Penman-Monteith model.

Daily Model Comparison
All candidate predictor variables in the LME model fit to daily E Pan were significant (p < 0.001) and fixed effects coefficients ranged from 0.402 for R sgd to 0.126 for U (Table 2).Variance inflation factors for the predictor variables were all well below 10.0 so we assumed multicollinearity did not confound interpretation of the coefficients.Large systematic differences, perhaps reflecting measurement error and 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 E Pan exhibited low precision (Figure 2).Coefficients of determination were similar across models.E Pen exhibited a model efficiency coefficient of 0.27, E PT exhibited a model efficiency coefficient of 0.33, E RFu and E RFc both exhibited model efficiency coefficients close to zero (i.e., no better than applying the mean as a predictor), while E PanFit exhibited a model efficiency coefficient of 0.40.In terms of overall bias, E Pen exhibited a mean error of 0.5 mm•day −1 (10%) (i.e., E Pen overestimated E Pan by 10 percent), E RG , E RFu , and E RFc exhibited larger biases (i.e., mean errors of 0.3, −0.1, and 0.1 mm•day −1 (7%, −3% and 2%), respectively, while E PT and E PanFit 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 pan evaporation measurements.The marginal coefficient of determination (i.e., without application of the random effects) was R 2 = 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 .

Monthly Model Comparison
Monthly variability of E Pan 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 E Pan and R sgd and this relationship showed some indication of curvature (Figure 3a).Relationships with T a and D a were weaker (Figure 3b,c), while there was no relationship between monthly E Pan and U (Figure 3d).
Consistent with fits previously reported at the daily scale, all candidate predictor variables in the LME model fit to monthly E Pan were significant (p < 0.001) (monthly LME model coefficients not shown).The marginal coefficient of determination for the monthly fit to pan evaporation was R 2 = 0.68, with a mean (observation minus prediction) error of −0.16 ± 0.06 mm•day −1 , and a root mean squared error of 0.93 mm•day −1 .
Models listed in decreasing order of performance (i.e., agreement with monthly E Pan measurements based on model efficiency coefficients), were E PanFit (MEC = 0.67), E PT (MEC = 0.64), E Pen (MEC = 0.53), E RG (MEC = 0.31), E RFu (MEC = 0.10), and E RFc (MEC = −0.20)(Figure 4).The rate of evaporation for a reference forests with canopies that are moderately coupled to the atmosphere substantially underestimated E Pan (Figure 4d).By design, these values were not expected to be unbiased expressions of PE, but we plotted them here to compare precision and explained variance for context.E PanFit slightly outperformed E PT , explaining 68% of the variance, with a mean error of −3% (Figure 4f).It is interesting to note that, despite slightly greater overall performance of E PanFit and E PT , these models still had limited capacity to reproduce the highest monthly values of E Pan above 7.0 mm•day −1 (Figure 4b,f).

Attribution of Variance in PE to Meteorological Drivers
In decreasing order of sensitivity, daily variability of E Pan was controlled by R sgd , T a , D a , and U (Figure 5).On a monthly time scale, the controls of T a and D a on E Pan dissipated and were compensated by increases in sensitivity to R sgd , while control by U was no longer significant.

Attribution of Variance in PE to Meteorological Drivers
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.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.

Spatial and Temporal Variability
Differences in the spatial variation of PE were evident from visual inspection of the mean warmseason 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).On a daily time scale, E Pen was most strongly driven by R sgd and secondarily by D a and U, while T a exerted a low positive effect.There were no significant differences in sensitivity between daily and monthly time scale for E Pen .E PT exhibited the strongest sensitivity to R sgd , with secondary dependence on T a .On a monthly time scale, control on E PT shifted slightly from R sgd to T a .Despite differences in aerodyne mic, the behaviour of E RFu and E RFc only showed modest differences in sensitivity to D a .Both exhibited strong dependence on D a , followed by R sgd and weak positive dependence on T a .Hence, none of the models perfectly matched the variance attribution of monthly E PanFit ; sensitivity of monthly E PT to T a was too high, while sensitivities of monthly E Pen , E RFu , and E RFc to D a were too high.

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 E Pen 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 E RG (Figure 6b).Reference forest values showed more prominent zonal gradients and the differences between spatial patterns of PE for moderatelyand 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 (E RG ) to 3.8 mm•day −1 (E PanFit ) (Table 3).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  On average over forest ecozones, there was little change in R sgd from 1951-2014 (Table 4).According to NCEP/NCAR global reanalysis, R sgd 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 E Pen , E PT , E RG , E RFu , E RFc , and E PanFit 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 Penman's model and the LME fit to pan evaporation measurements at 0.16 and 0.14 mm•day −1 respectively, 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 of aerodynamic and surface conductance, a general relationship was produced from regression analysis (Figure 9).The analysis confirmed that time trends increased strongly with increasing aerodynamic conductance and weakly with increasing surface conductance.

Discussion
Ecological studies commonly make use of basic expressions of evaporative demand.Here, we explored the behavior of several common approaches to calculating PE.While such basic expressions of PE may have strong practical value in studies that attempt to predict climate change impacts on terrestrial ecosystems across large spatial and temporal domains, further work is needed to improve the representation of spatial variation in meteorological conditions above forests, canopy structure and physiology, and ultimately the prediction of actual evapotranspiration rather than PE or reference rates [52,53].
We considered measurements of pan evaporation as one independent indicator of PE.While there are likely key differences between PE above forests and evaporation pans, the comparison helped place 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 were only rarely accompanied by on-site meteorological variables, so considerable effort was placed on deriving best available values from nearby climate stations, which constitutes an important source of 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 by applying a constant pan coefficient and including a random effect on the intercept coefficient at the scale of individual pan evaporation stations in the LME models.Applying a mid-range value of the pan coefficient, k = 0.77, led to a mean E PanFit of 3.8 mm•day −1 across North American forests.As this value was 0.5 mm•day −1 greater than the next highest value (E Pen ), any future consideration of the dataset should likely consider reducing k.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.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 forest ecosystem modelling.Analysis of variance attribution suggested that controls on E Pan were scale dependent.On a daily time scale, E Pan was significantly related to R sgd , T a , D a , and U, whereas on a monthly time scale E Pan was only significantly related to R sgd , T a , and D a .We interpreted this difference to mean that systematic effects of advection (as an energy source over pans) and ground heat flux (as an energy sink) diminish as time scale increases.That is, on any given day, advection and 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 of spatial and temporal variation (including trend analysis) may not suffer from failure to represent wind.
Varying sensitivity to D a 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 D a -sensitivity and strong D a -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 D a across North America [85].Interestingly, the different dependence on D a 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 E RFu and E RFc start to exceed E Pen when D a > 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 E PT on the basis that it does not include dependence on D a , 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 E RFu and E RFc , owing to strong dependence on increasing D a (and the general pattern of increasing trends with increasing g a and g s ) 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 E Pan .That is not to say, however, that the more modest changes in evaporative demand implied by E PT cannot impose severe water stress, as the sensitivity of plant health and productivity to changes in E PT , specifically, remains largely unexplored.At the regional scale, all trends in PE were positive with the exception of trends in the eastern temperate ecozone, where E PT was negative and E Pen was zero.Interestingly, this could be the signature of high tropospheric aerosol forcing on surface climates over the east-central conterminous U.S. (Figure 8).Similar patterns exist for pockets of industrial activity across the continent.Further research 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 onset of aerosol emission regulations that began late in our study period (i.e., 1990's), and conversely by any new or growing emission sources.

Conclusions
We investigated spatial and temporal variation of PE across forest regions of North America.While records of E Pan were not long enough to perform direct trend analysis, we instead chose to produce a statistical model, E PanFit , that helped to place the magnitude of model predictions in context and evaluate differences in the relative effect of the driving meteorological factors.Consistent with previous model inter-comparison studies [1,87,88], we found substantial differences in PE among 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 whether it was calculated with the Priestley-Taylor model, or with the Penman-Monteith model set for a reference 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 reasons for divergent behavior in quantitative terms.We conclude that until a consensus emerges on which expression is more representative of reality in forests, studies that adopt such basic expressions of PE should 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 outbreak risk, drought-related mortality, and health decline syndromes, indicates a need to improve monitoring of the key drivers of PE, including R sgd and D a .This source of high resolution, continuous monthly variation of PE from various models can facilitate research on a wide range of applications that require knowledge of the surface water balance.

Figure 1 .
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.

Figure 1 .
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.

Figure 2 .
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 mixedeffects model fit to daily pan evaporation measurements (EPanFit).Linear curves indicate the ordinary least squares fits.Statistics include: coefficient of determination (R 2 ), 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 .
Figure 2. Comparison between observations of daily total pan evaporation (E Pan ) (after applying the pan coefficient of 0.77) and models of potential evaporation based on: (a) Penman's model (E Pen ); (b) the Priestley-Taylor model (E PT ); (c) reference rate for grassland based on the Penman-Monteith model (d) reference rate for forests canopies that are moderately coupled to the atmosphere (E RFu ); (e) reference rates for forest canopies that are well coupled to the atmosphere (E RFc ); (f) linear mixed-effects model fit to daily pan evaporation measurements (E PanFit ).Linear curves indicate the ordinary least squares fits.Statistics include: coefficient of determination (R 2 ), 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 4 .
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 .
Figure 3. Relationships between monthly mean pan evaporation (E Pan ) and climate variables, including incident global shortwave radiation (R sgd ), air temperature (T a ), vapour pressure deficit (D a ), and wind speed (U).

Figure 4 .
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 4 .
Figure 4. Comparison between observations of monthly mean pan evaporation (E Pan ) 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 (E PanFit ).

Figure 5 .
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.

Figure 5 .
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 (E PanFit ); 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.

Figure 6 .
Figure 6.Comparison between long-term (1971-2000) mean daily potential evaporation (PE) over the warm season (May-September).Units of PE are expressed as mm•day −1 : (a) Penman's model (E Pen ); (b) Priestley-Taylor model (E PT ); (c) reference grassland (E RG ) based on the Penman-Monteith model; (d) reference forest with moderately coupled canopy (E RFu ) based on the Penman-Monteith model; (e) reference forest with coupled canopy (E RFc ) based on the Penman-Monteith model; (f) linear mixed-effects model fit to monthly mean pan evaporation (E PanFit ).

Figure 7 .Table 3 .
Figure 7. Zonal mean warm-season (May-September) potential evaporation from various models over the land surface of North America.

Figure 7 .Table 3 .
Figure 7. Zonal mean warm-season (May-September) potential evaporation from various models over the land surface of North America.

Figure 8 .
Figure 8.Long-term change in warm-season (May-September) mean potential evaporation (PE) from 1951-2014 (mm•day −1 ): (a) Penman's model (E Pen ); (b) Priestley-Taylor model (E PT ); (c) reference grassland (E RG ); (d) reference forest with moderately coupled canopy (E RFu ); (e) reference forest with well coupled canopy (E RFc ); (f) linear mixed-effects model fit to monthly mean pan evaporation (E PanFit ).Values were calculated by multiplying the slope coefficient from regression relationship between PE and time by the length of the study period (64 years).

Figure 9 .
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).

Figure 9 .
Figure 9. Dependence of long-term changes in reference rates of evaporation on levels of aerodynamic conductance (g a ) and surface conductance (g s ) 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).

Table 1 .
Station information: Identifiers (ID) from Environment Canada or NOAA, geographic coordinates, mean annual temperature (MAT) and mean annual precipitation (MAP) for the 1971-2000 base period, number of days with measurements of pan evaporation (E Pan ) used in the analysis (n), and time period of containing E Pan measurements.