Comparison of Latent Heat Flux Using Aerodynamic Methods and Using the Penman – Monteith Method with Satellite-Based Surface Energy Balance

A surface energy balance was conducted to calculate the latent heat flux (λE) using aerodynamic methods and the Penman–Monteith (PM) method. Computations were based on gridded weather and Landsat satellite reflected and thermal data. The surface energy balance facilitated a comparison of impacts of different parameterizations and assumptions, while calculating λE over large areas through the use of remote sensing. The first part of the study compares the full aerodynamic method for estimating latent heat flux against the appropriately parameterized PM method with calculation of bulk surface resistance (rs). The second part of the study compares the appropriately parameterized PM method against the PM method, with various relaxations on parameters. This study emphasizes the use of separate aerodynamic equations (latent heat flux and sensible heat flux) against the combined Penman–Monteith equation to calculate λE when surface temperature (Ts) is much warmer than air temperature (Ta), as will occur under water stressed conditions. The study was conducted in southern Idaho for a 1000-km area over a range of land use classes and for two Landsat satellite overpass dates. The results show discrepancies in latent heat flux (λE) values when the PM method is used with simplifications and relaxations, compared to the appropriately parameterized PM method and full aerodynamic method. Errors were particularly significant in areas of sparse vegetation where differences between Ts and Ta were high. The maximum RMSD between the correct PM method and simplified PM methods was about 56 W/m in sparsely vegetated sagebrush desert where the same surface resistance was applied. OPEN ACCESS Remote Sens. 2014, 6 8845


Introduction
The Penman-Monteith (PM) method for estimating latent heat flux (λE) calculation has long been extensively used in hydrological, atmospheric and environmental modeling.In 1948, Penman [1] combined aerodynamic bases and energy bases for calculating evaporation, which eliminating the use of the most complex meteorological parameter to be measured in the field, i.e., surface temperature.The development of satellite-based techniques has facilitated the computation of surface temperature, though data are not always available in the required temporal and spatial resolution.Mallick et al. [2] demonstrated a method that physically integrates the radiometric surface temperature into the PM equation for estimating terrestrial surface energy balance fluxes.The PM method computes relatively accurate latent heat flux (λE) in different meteorological and surface conditions when the correct parameterizations are used.However, it is commonly known to the hydrological community that the PM method can produce errors in latent heat flux calculations in sparsely vegetated areas due, in part, to: (1) the assumption of neutral atmospheric stability common to many applications for reference ET; (2) inaccurate calculation of the slope of the saturation vapor pressure curve; (3) inaccurate calculation of long-wave radiation emission from the surface; (4) underestimation of soil heat flux; and (5) mismatches between sources and sinks of momentum, heat and vapor fluxes; see, for example, Tanner and Pelton [3], Wallace et al. [4], Stannard [5], Moran et al. [6], etc.This study reinforces previous work on the difficulty in calculating λE using the PM method in sparse vegetation using satellite and weather-based data.It is essential to appropriately parameterize the PM equation when it is used to estimate λE in areas of sparse vegetation.While calculating λE using the PM method on a point scale, air temperature (T a ) is generally collected from meteorological stations at designated heights above the surface (generally ~2-3 m).These meteorological stations generally do not record surface temperature (T s ), and many studies have based important parameters on the assumption of similarity between T s and T a when computing λE using the PM method.Even though T s is not required to compute λE using the PM method, assumptions implicit to the PM can lead to error if differences between T s and T a are large.
Many studies conducted sensitivity analyses on parameterization and appropriate adaptation of the PM method (Thom and Oliver [7], De Bruin and Holtslag [8], Shuttleworth and Wallace [9], Allen [10], Kustas et al. [11], Moran et al. [6], Alves and Pereira [12], Gavin and Agnew [13], Nandagiri and Kovoor [14], Allen et al. [15], Gong et al. [16], Zhao et al. [17], etc.).Otles and Gutowski [18] studied atmospheric stability effects on the PM method and compared the evapotranspiration from PM to the Bowen-ratio and lysimeter data.Mukammal et al. [19] compared aerodynamic and energy balance techniques to estimate evapotranspiration from a cornfield where they discussed difficulties in using the aerodynamic method in tall crops.Results from Otles and Gutowski [18] showed that under extremely unstable boundary layer conditions, actual evapotranspiration values become smaller than estimated by the original PM formula, because under the same net input energy, the sensible heat flux is considerably larger than under neutral conditions.Thom and Oliver [7] extensively discussed the variation in estimated λE with the stability of the ratio of actual to neutral aerodynamic resistance.Allen et al. [15] discussed representative values for the bulk stomatal leaf resistance (r l ) for neutral conditions common to the so-called reference crop surfaces of alfalfa and clipped, cool season grass.Mahrt and Ek [20] discussed the influence of atmospheric stability on potential evaporation rates.Lhomme et al. [21] compared effective parameterization of the surface energy balance using various methods and conditions in heterogeneous landscape.The main objective of this study is to emphasize the advantages of using surface energy balance via separate aerodynamic equations versus using the PM combined equations while calculating λE when surface T s is significantly warmer than T a .One specific objective of this study is to explore the applicability of the fully-populated PM method and the neutral boundary-layer PM method (PM nbl ) in dry and water-limited environments.These objectives were implemented through the assistance of remote sensing techniques and data.
The availability of high-resolution gridded weather and remote sensing-based data has opened up important opportunities to understand and compute fluxes for different surface and atmospheric conditions.In this study, bulk surface resistance (r s ) was computed from inverting the fully-populated PM and aerodynamic (AERO) methods for latent heat flux (λE) using instantaneous ET (ET ins ) obtained from the Mapping Evapotranspiration at High Resolution using Internalized Calibration (METRIC) Landsat-image based surface energy balance.The analyses either use the same latent heat flux (λE) from the METRIC model when Landsat-based T s is used or modifies the METRIC-generated λE based on iteratively calculated T s .The impact of using iteratively calculated surface temperature (T s ) was examined to understand its influence on the λE calculation in both the λE PM and λE aero approaches.The results of the comparisons between the λE PM and λE aero methods show that both methods produce identical results when Δ (f = T s , T a )) is used in the PM method and other parameters, including net radiation (R n ), soil heat flux (G) and buoyancy, are based on full, physically-based computations.The use of the surface energy balance combined with remote sensing techniques facilitated sensitivity analysis for the computation of these fluxes, including λE.A backward-averaged accelerated numerical solution (BAANS, Dhungel et al. [22]) was used to simulate surface energy balance flux parameters based on gridded weather and satellite-based data for two satellite overpass dates under clear sky conditions, although the analysis can be extended to any atmospheric and surface condition.The description and computation procedure of the surface energy balance flux parameters using BAANS is shown in Figure 1 and detailed in Appendix A. BAANS utilizes instantaneous evapotranspiration (ET ins ) and other vegetation indices from the METRIC (Allen et al. [23]) model and gridded weather data from North American Regional Reanalysis (NARR) (Mesinger et al. [24]) (Figure 1).
The instantaneous ET estimates, ET ins , acquired from METRIC (Allen et al. [23]) for Landsat overpass times are based on a measured radiometrically-based T s ; the latent heat of vaporization (λ) calculated in BAANS is based on an iteratively calculated T s in order to apply the model for periods between satellite overpasses.BAANS simulates T s iteratively inside the surface energy balance based on meteorological conditions and surface roughness by applying the Monin-Obukhov (MO) similarity theory and stability corrections.The major objective of BAANS is to estimate the surface energy flux parameters, including T s , using a single, near-surface layer from available gridded weather sets (NARR, North American Land Data Assimilation System (NLDAS), etc.) and satellite-based roughness data (NDVI, LAI, α, ε o , etc.), where the near-surface layer, at 30 m, is used as a type of "blending layer" for air temperature, specific humidity and wind speed that is extrapolated to specific surface features at the 30-m scale.This procedure gives an opportunity to compute λE over large areas when no thermal-based T s data are available.In this study, BAANS utilizes ET ins from the METRIC model during parameterization of surface resistances to reduce the number of iterations in the convergence process of the surface energy balance, as well as to serve as a boundary condition of latent heat flux through ET ins .It further evaluates the accuracy of simulated T s based on the METRIC-based ET ins .When BAANS is implemented for non-satellite overpass times, ET ins is computed iteratively inside the surface energy balance along with T s .In this study, BAANS provides an opportunity to compare the surface resistance (r s ) and λE obtained from two independent and different approaches.A central question still remains regarding the accuracy of the simulated T s from BAANS and its possible implications for the final results and conclusions.To address this, Section 4 provides an analysis based on both thermal-based T s and simulated T s .Both methods (AERO and PM) produced identical results, irrespective of the T s simulated from BAANS.

Methodology
The procedures for calculating latent heat flux (λE) from the aerodynamic method and PM method are discussed in this section.

Aerodynamic Method of Calculating Latent Heat Flux (λE aero )
The aerodynamic method for computing latent heat flux (λE aero ) can be expressed in a resistance form as Equation (1) (for example, Van de Griend and Owe, [25], Lhomme et al. [21], etc.).
where λE aero is latent heat flux (W/m 2 ) from the aerodynamic method, ρ a is atmospheric density (kg/m 3 ), C p is the specific heat capacity of moist air (J/kg/K), e°s ur is the saturation vapor pressure of the evaporating surface (kPa), e a is the actual vapor pressure of the air at the same height as the measured T a (kPa), r s is the bulk surface resistance (s/m), r ah is the bulk aerodynamic resistance (s/m) between the surface and height of e a and T a measurements and γ is the psychrometric constant (kPa•°C −1 ).The saturation vapor pressure at the surface (e°s ur ) in Equation ( 1) is computed using Equation (2) (Tetens [26]).In this study, surface temperature (T s ) was computed iteratively inside the surface energy balance by inverting the aerodynamic equation of sensible heat flux (H) using Equation (A2), coupled with equations for R n , G, where H = R n − G − λE.
where T s is in Kelvin (K).In the developed surface energy balance, outgoing longwave radiation (R L↑ ) is computed from T s using Equation (3).Most of the currently operational remote sensing-based evapotranspiration models, such as surface energy balance algorithm for land (SEBAL) (Bastiaanssen et al. [27]), METRIC (Allen et al. [23]) and Surface Energy Balance System (SEBS) (Su [28]) utilize satellite-radiatively-based T s in the surface energy balance.In Equation (3) (BAANS), T s is computed iteratively by inverting the aerodynamic equation of sensible heat flux (Equation (A2)).
where R L↑ is outgoing longwave radiation (W/m 2 ), ε o is broadband emissivity and σ is the Stefan-Boltzmann constant (W/(m 2 •K 4 ).In one of the comparative scenarios, the impact of using T a to calculate R L↑ is also investigated (Equation ( 4)).
where T a is air temperature (K).Equation ( 4) is frequently used in hydrological applications as a relaxation of data requirements.As shown later, using T a for outgoing longwave radiation can cause substantial error in the R n and λE estimates for sparse and/or water-stressed vegetation.In BAANS, sensible heat flux in the surface energy balance is computed as a residual, as shown in Equation (5).
where H is sensible heat flux (W/m 2 ), R n is net radiation (W/m 2 ), G is ground heat flux (W/m 2 ) and λE ins (separately λ ET ins ) is latent heat flux from METRIC (W/m 2 ).R L↑ of Equation ( 3) is used to compute net radiation (R n ) (Equation ( 6)).
where R s↓ is incoming global solar radiation (W/m 2 ), R L↓ is downwelling longwave radiation (W/m 2 ), R L↑ is outgoing longwave radiation (W/m 2 ) and α is surface albedo.Bulk surface resistance (r s_aero ) was computed on a 30-m pixel basis for the Landsat image domain by inverting the aerodynamic equation of latent heat flux (Equation ( 1)) using λE from METRIC.Shuttleworth and Wallace [9] used a similar inversion approach to compute the surface resistance at the substrate surface.

The Penman-Monteith Method for Calculating Latent Heat Flux (λE PM )
The Penman-Monteith (PM) equation (Monteith [29]) for computing latent heat flux (λE PM ) is shown in its general form in Equation (8), as given, for example, by Allen et al. [30].The PM equation, Equation (8), reverts to a suite of aerodynamic and radiation-based equations describing the full surface energy balance (including Equations (1), (3) and (A15)) when T s is used along with T a when calculating the slope of saturation vapor pressure curve (Δ (f = T s , T a )).This reversion is shown in Appendix B.
where Δ is the slope of the saturation pressure (kPa•°C −1 ), R n is the net radiation (W/m 2 ) and G is ground heat flux (W/m 2 ).Valiantzas [31] discussed the difficulties of calculating evapotranspiration using the PM method, largely because of the difficulties of measuring or estimating weather variables, including net radiation, vapor pressure deficit, etc.The saturation vapor pressure of air (e°a ir ) was computed in this study using T a at 30 m from the NARR reanalysis data set as in Equation (9) (Tetens [26]).
where e°a ir is in kPa.As mentioned earlier, the ease of use and popularity of the PM method is high, since, in many applications, it does not directly depend on T s, to calculate evapotranspiration.However, as previously discussed, applying the PM method under the PM nbl assumptions that essentially presume a reference evapotranspiration (ET r ) condition (Allen et al. [15]; American Society of Civil Engineers-Environmental and Water Resources Institute (ASCE-EWRI) (Walter et al. [32]), where the majority of available energy is partitioned into λE, i.e., λE~(R n − G), so that sensible heat flux is nearly equal to zero, can lead to large errors in estimated λE for dry surfaces, due to the assumption that T a ~Ts .It is only under conditions of a large λE rate that the assumption that T a ~Ts is valid and the PM nbl does not produce large errors in λE.While calculating evapotranspiration in sparse vegetation using the PM method, the calculation of ∆ varies largely compared to fully-vegetated and well-watered agricultural land; this is because of the large temperature difference between surface and air in sparse vegetation.In the full PM application, the slope of saturation vapor pressure (Δ) is calculated through a linear relationship of the saturation curve between T s and T a (Equation ( 10)).
This slope of saturation vapor pressure value (∆ (f (= T s , T a )) is appropriate for sparse vegetation and/or drying surfaces, as it incorporates the influence of surface temperature (T s ) on the extrapolation of air temperature to the surface, implicit to the Penman combination equation derivation (Penman [1] Allen et al. [15]); however, Equation (10) requires T s , which is generally not measured.The estimate for Δ using T a only is shown in Equation (11), which has been widely used to calculate λE in the PM method (Allen et al. [30], Xu and Singh [33], Nandagiri and Kovoor [14], Donatelli et al. [34], etc.).
In this study, the impact of iteratively determining T s was explored, where λE is calculated using the PM method with Δ from Equation (10) (∆ (f = T s , T a )) and from Equation (11) by Murray [35] (∆ (f = T a )) was tested for use as an appropriate procedure.Most of the previously discussed simplifications do not account for Equation (10), while computing ∆ in the PM method and using ∆ approx (f = T a ) as the appropriate procedure for calculating ∆.
Apparent values for bulk surface resistance for use in the PM equation (r s_PM ) were computed by inverting the PM equation (Equation ( 12)) using λE from METRIC.These values for r s_PM are compared to r s_aero (Equation ( 7)) below in Section 4.1 to evaluate the consistency in estimates.

Aerodynamic Resistance (r ah )
Aerodynamic resistance for the neutral atmospheric condition (r ah_neu ) was computed from Equation (13) (Thom [36], and Monteith and Unsworth [37]).The impact of using r ah_neu in the PM method is discussed in one of the scenarios.
Equation ( 14) is used to calculate r ah under conditions requiring stability correction (Thom [36]).
where ψ m is a stability correction term for momentum, ψ h is a stability correction term for heat, L is MO length (m), z is the measurement of the height of wind and temperature (m), Z om is the roughness length of momentum (m), Z oh is the roughness length of heat (m), u z is wind speed (m/s), d is zero plane displacement (m) and k is the von Kármán constant, i.e., 0.41.For simplicity, the roughness length of heat transfer (Zoh) is computed as 0.1 times Zom for agricultural land use classes, i.e., Zoh~0.1 × Zom (m).For sparse vegetation areas, like sagebrush desert and grasslands, the relationship is not kept constant.BAANS computes r ah using Equation ( 13) and ∆ using Equation ( 11) for numerical stability when differences between T s and T a are less than 0.01 K (near neutral atmospheric conditions).
In an ideal scenario, the λE obtained from the aerodynamic method (λE aero ) and the Penman-Monteith (λE PM ) method should be identical.For example, if identical latent heat flux (λE), aerodynamic resistance (r ah ) and surface energy flux parameters are provided for Equation (1) and Equation ( 8), the inverted bulk surface resistance (r s ) should be identical for both methods.Because of the various assumptions in the surface energy balance and surface roughness, as well as the closure problem in energy balance, it is not always possible to produce identical results from these methods.However, when utilizing surface energy balance, BAANS facilitated simulation of identical r s values through the use of the same equations for aerodynamic roughness, R n and G. Section 4.1 shows the comparison of r s between the two models resulting in identical values from both methods.

Relaxations in the Parameterization of the Penman-Monteith Method
In this section, evaluations are carried out for the most frequently used relaxations while calculating λE using the PM method.These relaxations concern the calculation of ∆, assumptions in stability correction, the calculation procedure of R L↑ and a combination of these.Table 1 lists different parameterizations and assumptions that represent a combination of relaxations commonly.Table 1 also includes associated parameter equations used while calculating λE from the aerodynamic method that can change during relaxations in the PM method.λE aero was produced by the aerodynamic method using Equation ( 1) with stability correction.The application of λE aero followed the complete surface energy balance procedure outlined in Appendices A and B and in Figure 1.The results of relaxations in Table 1 are labeled (λE PM_Δra , λE PM_Δ and λE PM_ΔRL ), representing when Δ and r a are estimated assuming neutral stability (T s = T a ), when Δ is estimated assuming neutral stability, but r a is estimated using iteratively estimated T s , and when Δ, r a and R L↑ are estimated assuming neutral stability, respectively.The λE PM_ΔRL product is commonly (and improperly) applied in practice.
The fully-parameterized PM, labeled λE PM , is produced by the PM method using Equation ( 8) with stability correction, and ∆ and R L↑ are calculated with full use of iteratively derived T s (Equation ( 10)).λE PM is considered to be the correct estimate of latent heat flux from the PM method.As both λE aero and the fully-parameterized λE PM are calculated using converged surface energy flux parameters based on the iterative solution of T s , they are expected to be identical, because the PM method reverts to the surface energy balance when ∆ (f = T s , T a ) (Equation ( 10)) is used.In the fully-populated application of the PM and aerodynamic procedures, ground heat flux (G) estimation considered the influence of T s when T s and T a were different.The estimation of G was based on recent algorithms from METRIC (Allen et al. [38]).7) Equation ( 1) Equation ( 3) Equation ( 14) λE PM Equation ( 8) Equation ( 10) Equation ( 3) Equation ( 14) Equation ( 12) ---λE PM_Δra Equation ( 8) Equation ( 11) Equation ( 3) Equation ( 13) Equation ( 12) ---λE PM_Δ Equation ( 8) Equation ( 11) Equation ( 3) Equation ( 14) Equation ( 12) ---λE PM_ΔRL Equation ( 8) Equation ( 11) Equation ( 4) Equation ( 14) Equation ( 12) --- The term "error" refers to an underestimation or overestimation of λE compared to the appropriately (i.e., fully) parameterized PM method.Both overestimation and underestimation of λE with the relaxations create problems in the closure of the surface energy balance, even if accurately measured R n and G are used in the PM method.
In the first simplification, λE PM_Δra was computed assuming a neutral atmospheric condition where stability correction factors Ψ m and Ψ h are assumed to be zero (Equation ( 13)).In this simplification, r ah is replaced by r ah_neu and ∆ is replaced by ∆ approx in the PM equation, while the rest of the parameters and fluxes (i.e., R n , G and r s ) come from the surface energy balance.λE PM_Δra is computed to understand the consequences of not applying stability correction while computing λE using the PM method.This PM application includes some bias from common practice, because R n , G and r s may, in practice, also be affected by the simplifications in the surface energy balance, which are not accounted for.
In the second relaxation, while calculating λE PM_Δ from the PM method, the slope of the saturation vapor pressure curve ((∆) (f = T a , T s )) of λE PM was replaced by ∆ approx (f = T a ), and the rest of the parameters (i.e., R n , G, r ah and r s ) are utilized from the fully-parameterized λE PM .In this simplification, the consequence of using ∆ approx (f = T a ) with converged surface energy balance flux parameters and stability corrected r ah is explored.This is one of the frequently adopted relaxations for calculating λE from the PM method; however, knowledge of surface temperature is required to estimate buoyancy correction and outgoing long-wave radiation.
In the third relaxation, while calculating λE PM_ΔRL , stability correction is based on iteratively determined T s , where the iteration is carried out for the rest of the surface energy balance flux parameters until convergence is achieved.In the λE PM_ΔRL relaxation, R n and G are based on R L↑approx (f = T a ), and ∆ is computed using ∆ approx (f = T a ).The result of this simplification is the combined effect of R L↑approx (f = T a ) and ∆ approx .This relaxation application faces bias from common applications because λ and G are calculated based on T s in the iteration, but the rest of the surface energy balance flux parameters are based on T a .In this relaxation, a separate surface energy balance needs to be conducted, and the flux parameters are different than the rest of the λEs (λE aero , λE PM , λE PM_Δra and λE PM_Δ ), because of the use of T a in the surface energy balance.
The fourth relaxation is where T s is neglected in all parameters in the surface energy balance.There was little variations in results between this and the third relaxation (λE PM_ΔRL ), so the discussion of this relaxation is not carried out in detail.

Study Area and Data
This study was conducted in a semi-arid environment in southern Idaho for a study AOI of about 1000 km 2 , consisting of a variety of land use classes (Figure 2).The representative study area consists of about one million Landsat pixels (each one a single NARR reanalysis pixel), including irrigated agricultural, grassland, sagebrush desert areas and open water bodies.NARR reanalysis (NARR, Mesinger et al. [24]), gridded weather data and satellite-based METRIC (Allen et al. [23]) data were used to compute and close the surface energy balance.Landsat 5 images from 17 May 2008 (11:02 am), and 18 June 2008 (11:01 am), were processed using the METRIC model, where products included leaf area index (LAI), surface albedo (α), broadband emissivity (ε o ) and instantaneous evapotranspiration (ET ins ).There were generally crops growing during this period in agricultural parts of the image, and the application of irrigation was common.Desert soils may still have had some moisture because of winter precipitation.
A short description of the parameters used in the METRIC model is given in this section.The METRIC model computes ET ins based on the surface energy balance as a residual of the energy balance.While calculating ET ins , METRIC uses a water balance model to infer the impact of antecedent soil moisture on drier (and hotter) parts of the image.BAANS uses ET ins from METRIC for calibration and initialization, so that it is impacted by the antecedent soil moisture condition.When BAANS is used for time periods in between the satellite overpass dates, a soil water balance tracks soil moisture.METRIC and the similar SEBAL model (Bastiaanssen et al. [27]) calculate reasonably accurate ET ins , but can be impacted by the behavior of population end-members during calibration: the surface temperatures of hot and cold pixels and the available energy for the hot pixel.Although meteorological variables are necessary components of SEBAL-type models, the most critical variables that determine the performance of SEBAL/METRIC are surface temperatures and net radiation at end-members as selected by the operator.In the METRIC model, the aerodynamic function of sensible heat flux is based on the near-surface temperature difference (dT), i.e., the temperature difference between two near-surface heights.METRIC computes aerodynamic roughness Z om based on the leaf area index (LAI) for agricultural land and uses the Perrier [39] Z om function for tall grass and desert areas (Allen et al. [23,40]).Broadband emissivity (ε o ) is computed from LAI. Surface albedo (α) is calculated by integrating reflectivities from six shortwave bands of Landsat.Bastiaanssen et al. [27] and Allen et al. [23] provide details on computation of ET ins from SEBAL and METRIC, respectively.
The central idea of this study was to utilize commonly available gridded weather data (NARR reanalysis, NLDAS, etc.) and remote sensing data to calculate λE based on the full form of the surface energy balance; a similar method has been demonstrated in Dhungel et al. [22], Dhungel and Allen [41].Table 2 shows weather data from NAAR reanalysis for the two satellite dates, 17 May and 18 June 2008, for one 32-km NARR cell.Instantaneous weather data at 11:00 am from NARR reanalysis were used for the surface energy balance near satellite overpass time.NARR reanalysis data used from the 30-m height of NARR included air temperature (T a ), wind speed (u z ) and specific humidity (q a ).The 30-m height functioned as an effective blending height for extrapolating these data to surfaces of individual 30-m Landsat pixels.NARR reanalysis data for the surface included incoming shortwave radiation (R S↓ ) and incoming longwave radiation (R L↓ ).The incoming solar radiation (R S↓ ) was relatively high at the satellite overpass times for both dates due to clear sky conditions.Those clear sky values created a large difference between the air and surface temperature for dry Landsat pixels.NARR reanalysis has a large amount of additional gridded weather data and layers other than those utilized in this study, including surface temperature, outgoing longwave radiation and outgoing shortwave radiation.However, incoming surface energy fluxes input from NARR reanalysis data have much coarser resolution than Landsat, and the variation of these meteorological parameters within the 32-km resolution pixel is assumed and precludes the use of T s and the emitted longwave radiation at the 30-m scale, where large ranges in water availability and, thus, T s are common.For example, localized precipitation or irrigation events, as well as surface roughness, albedo and emissivity create variation in soil moisture in the surface and root zone and, thus, outgoing energy in each Landsat pixel.As this study aimed to compute the surface energy flux parameters at the higher spatial resolution of Landsat, these variations needed to be captured using a surface energy balance at the 30-m scale.The use of the coarser resolution NARR reanalysis data to represent surface conditions and outgoing fluxes would fail to capture these variations among Landsat pixels.The question still remains whether 32-km coarse spatial resolution weather data at 30 m from NARR reanalysis can represent the incoming meteorological inputs of all approximately one million 30-m resident pixels in the surface energy balance.To assess the impact of using large-scale weather data from NARR reanalysis, three separate analyses were conducted.In the first, NARR data from neighboring cells to the AOI were compared with that for the study AOI to evaluate spatial variation.An example is shown in Figure 2. The results suggest that the NARR reanalysis data contain negligible variation in neighboring pixels in the relatively dry Idaho climate.Secondly, NARR reanalysis data were extensively compared to various AgriMet and remote automatic weather stations (RAWS) for agricultural and desert environments, respectively, every 3 h of temporal resolution throughout the year 2008.The meteorological data compared included air temperature, wind speed, relative humidity, vapor pressure and solar radiation.Finally, theoretical clear sky radiation (R so ) from the Reference Evapotranspiration Calculator (RefET) (Allen [42]) model was compared to solar radiation from NARR reanalysis data sets.The analyses indicate that NARR reanalysis approximated actual surface meteorological measurements relatively closely.Figure 2 shows the study area using a Landsat image from 18 June 2008, and NARR reanalysis image for the same date.The bottom left side of Figure 2 shows NARR reanalysis incoming longwave data (W/m 2 ) for 18 June 2008, which is overlain by an Idaho county border map.The upper left figure is an expanded view of the NARR reanalysis data used in the study AOI.Incoming longwave radiation from NARR reanalysis varied between 273 and 322 W/m 2 in the expanded AOI with about 50 W/m 2 difference.Two sub-AOIs, one in agricultural land and another in sagebrush desert, were selected for statistical analysis in the 32-km study area (Figure 2).
The main objective of this paper is to emphasize the advantage of using separate aerodynamic equations versus the combined PM equation in the surface energy balance studies, when T s is much higher than T a .Therefore, the size of the study area and time periods of the study were constrained for this study to focus on the behavior of the methodology (e.g., Penman-Monteith), rather than on large spatial data coverage.The representative area and time periods have sufficiently explained the behaviors of the both methods.It is assumed that the meteorological data acquired from NARR reanalysis will vary little within 32-km resolution.In addition, any bias in METRIC-generated parameters and ET ins will also be carried out in this study, but these biases do not directly influence the objectives of the study.

Apparent Bulk Surface Resistance (r s )
One question related to the relaxation of various environmental parameters in the PM equation is the effect of the relaxation on the value for effective bulk surface resistance (r s ) required to reproduce the same λE as produced by the full PM and full aerodynamic methods.Required values for r s depend on the degree of relaxation, for example neglecting the effects of T s on outgoing long-wave radiation or on Δ, or neglecting buoyancy correction.A large change in effective r s may point to increased risk of error in estimated λE, as differences between T s and T a increase.In this study, apparent, effective values for r s were computed by inverting the aerodynamic (λE aero , Equation ( 7)) method and PM method (λE PM , Equation ( 12)) and compared.Computations were made at Landsat over pass times using ET ins from METRIC and simulated T s derived through iterative calculation of the surface energy balance where λE from METRIC was an energy balance target.Satellite-based surface energy balance models, like METRIC and SEBAL, do not compute nor use r s to calculate λE, but rather calculate λE as a residual of a thermally-driven surface energy balance.Figure 3 shows the variation in effective r s for the study area on 18 June 2008, as determined by inverting the full PM method.Values for r s_PM ranged from about 20 to 300 s•m −1 for irrigated fields and from 500 to 2000 s•m −1 for rangeland areas.
As discussed earlier, when the PM model is applied with outgoing longwave radiation, Δ, G and buoyancy correction, all estimated using T s , then r s derived from inversion of the aerodynamic and PM methods should be identical.Tables 3 and 4 show the results of the surface energy balance for single pixels sampled from a range of land use classes at the two satellite overpass times.The computed r s from the full PM method shows insignificant variation from that computed using the full aerodynamic method.This comparison confirms that effective r s derived from both methods will be equal despite the accuracy of other parameters computed from BAANS if identical surface energy flux parameters are input to both methods.In the calculations, both methods utilized identical values for ρ a , C p , e a , λE, γ, r ah, R s , R L↓, G, T a and e°a ir .The full PM method utilized ∆ from T s (Equation ( 10)), whereas, the aerodynamic method used an iteratively determined e°s ur .As an example, the r s of an agricultural pixel (5,103,099 m, 3,807,774 m) having low LAI (0.063) and with an apparently relatively dry soil surface was about 820 s•m −1 when using the PM method and 823 s•m −1 when using the aerodynamic method on 18 June 2008 (Table 3).The bulk surface resistance (r s ) of another agricultural pixel (5,089,869 m, 3,826,134 m) having high LAI (5.65) was about 33 s•m −1 for both methods on 18 June 2008 (Table 3).That pixel had about 298 K simulated T s , with the NARR-based air temperature ~296 K and a sensible heat flux (H) of about 38 W•m −2 .When the ET ins specified from METRIC was high, for example, 0.88 mm•h −1 on 18 June 2008, for an agricultural pixel (5,085,190 m, 3,837,433 m) having an LAI of about 3.5, the value of r s became very small (near zero) for both models.Of the compared pixels (Tables 3 and 4), that pixel had the smallest value for H, as the difference between T a and T s was small.For an agricultural pixel with low LAI (about 0.08), but an apparently wet soil surface, so that ET ins from METRIC was about 0.67 mm•h −1 , the r s was about 93 s•m −1 .As a comparison, Food and Agriculture Organization (FAO-56) recommends an r s value of about 70 s/m for the well-watered grass reference evapotranspiration (ET o ) definition, based on the bulk stomatal resistance of the vegetation per unit LAI (r l ) of 100 s•m −1 for 24-h calculation time steps.The r s value will generally be smaller for hourly time steps than for daily time steps (Allen et al. [15]).Values for r s from the inversion of the aerodynamic and full PM method are in relatively good agreement with the standardized reference r s value.
In terms of iteratively determined T s , BAANS simulated a lower T s in desert and grassland areas as compared to the thermally-measured T s from Landsat using initially-specified parameterizations for Z om and Z oh , while estimated T s for agricultural land was in good agreement with Landsat-measured values (Tables 3 and 4).One of the uncertain parameters that has to be specified for the application of the surface energy balance and MO stability correction is surface roughness (Z om , Z oh and their interrelationships).The values for these parameters affect the aerodynamics of the energy balance (Beljaars and Viterbo [43], Mascart et al. [44], etc.) and ultimately the simulated T s from BAANS.After modifying surface roughness values (Z om , Z oh and their interrelationship) in desert and grassland, BAANS was able to simulate T s that compared better to the thermal-based T s for sagebrush desert and grassland.With this modification, the value estimated for r s was altered, even though both methods generated identical estimates before and after modifications to roughness (Tables 3 and 4).Table 3 shows values for r s for land use Class 52 (sagebrush) on 18 June 2008, where LAI was 0.043 and when Z om and Z oh were set at initial values of 0.3 and 0.03 following Allen et al. [40].r s averaged about 599 s/m from both methods, with simulated T s ~305 K as compared to a mean T s for the desert AOI of about 320 K.After modification of values for Z om and Z oh to values of 0.045 and 0.00045, r s determined from the inversion of the PM method took on values of about 1400 s/m for both methods, with simulated T s ~319 K close to the measured T s of 320 K. Tables 3 and 4 show insignificant variation in computed latent heat flux, even though the T s varied significantly after modified parameterizations.This is because ET ins from METRIC was specified as a boundary condition in both initial and corrected parameterizations of the surface energy balance.Other than ET ins , the rest of the surface energy fluxes varied (i.e., R n , G, H) based on the new simulated T s .In this study, the final values from BAANS, after modifications, were used to compare λE computed using the relaxed parameterizations in the PM method, including the r s derived by inverting the full PM application.A detailed discussion of the results of BAANS is presented in Dhungel et al. [22].It is worth noting in Table 3 the differences in values for R n that reach nearly 100 W•m −2 , even for similar albedos, due to the impact of differences in T s on outgoing long-wave radiation.In a similar manner, values for Δ differ by up to 80% due to differences in T s between cooler agricultural pixels and hotter rangeland pixels.
Figure 4 shows a scatter plot of r s derived by inverting the full PM versus r s derived by inverting the aerodynamic method for about 27,000 pixels in sagebrush desert and agricultural land AOI's for two different dates (18 June and 17 May 2008).The results show higher values for r s in the sparsely vegetated areas (Figure 4a), as expected, and lower values for r s in the agricultural area, where LAI values were higher (Figure 4b).Additionally, when agricultural areas have low LAI due to bare soil conditions, r s varies widely, exceeding the r s of desert (Figure 4c,d).The majority of points lay on the 1:1 line, showing that both methods result in nearly identical values of r s (Figure 4, Tables 3 and 4).The results presented in Figure 4 represent results after the modification of Z om parameterizations in BAANS with simulated T s .

Latent Heat Flux (λE)
Comparisons of λE are carried out in this section using the applied relaxations to the PM method to assess the impacts of ∆ (f = (T s , T a ) and ∆ (f = (T a ) on ET estimates.Limited studies have made rigorous quantitative assessments of the impact of relaxations, for example Paw U [45].A useful way to make the assessment is to cross compare estimates produced by relaxations to the PM method with calculations from the fully and properly implemented PM method or using a full surface energy balance-aerodynamic procedure, both of which are utilized in this study.Identical r s calculated from the converged surface energy balance flux parameters (λE areo and λE PM , Figures 1 and 4) is used to calculate λE (λE PM_Δra, λE PM_Δ, λE PM_ΔRL ) assuming that those values for r s represent the best and most true current values for r s based on current atmospheric and soil water content conditions at the satellite overpass time.In this section, the comparison of λE is carried out using both simulated and Landsat-based T s values.Thermal-based T s from Landsat is used to verify the accuracy of simulated T s and to explore the sensitivity of λE on error in simulated T s .The results from Section 4.1 confirmed that a robust and workable parameterization is obtained from the full PM method, where T s is used to calculate outgoing longwave radiation, the slope of the saturation vapor pressure curve, the soil heat flux and the boundary layer stability.Estimates by the full PM application are denoted as λE PM .The primary surface energy balance fluxes of this study are driven by NARR reanalysis gridded weather data using vegetation indices from the METRIC model (Section 4.1).This process utilizes ET ins from METRIC in the surface energy balance to assess values for r s determined by inversion of the PM model when applied with various relaxations.Figure 5 shows the λE PM based on the simulated T s for 17 May and 18 June 2008, which is used to compare to results of relaxations to the PM over an extended area that covers both agricultural and desert AOIs.The use of ∆ approx and the assumption of neutral atmospheric conditions are commonly practiced simplifications or relaxations when computing λE from the PM method; see, for example, Lankreijer et al. [46], Penman [1], etc. Figure 6 shows Landsat-based T s , simulated T s from BAANS and differences between Landsat-based T s and NARR-based air temperature, T a , respectively, on 18 June 2008.Simulated T s from BAANS followed Landsat-based T s relatively closely for all land use classes (Figure 6b). Figure 6c shows differences between Landsat-based T s and T a that was very large, up to 30 K, in desert and grassland areas, as well as in the sparsely vegetated agricultural land.NARR-based T a at 30 m was about 296 K at the Landsat overpass time on 18 June 2008.
The neutral stability condition basically assumes sensible heat flux (H) to be zero, partitioning the available energy into latent heat flux (λE).In a neutral atmospheric condition, aerodynamic resistance (r ah_neu ) is larger than the stability-corrected r ah in an unstable condition.Lower u z and smaller Z om values also increase the value of r ah for a neutral atmospheric condition.The slope of saturation vapor pressure (∆ (f = (T s , T a )) is equal to ∆ approx (f = T a ) for the neutral atmospheric condition, because T s and T a , by definition, are nearly identical.The slope of saturation vapor pressure (∆ approx (f = T a )) is smaller than ∆ (f = (T s , T a ) for the unstable condition where values for T s are greater than for T a .The slope of saturation vapor pressure (∆) exponentially increases with increasing T s in Equation ( 10) (Monteith [29]), while ∆ in Equation ( 11) remains constant based on T a .As both r ah and ∆ are in the numerator and denominator of the PM equation, the impact of these relaxations (larger value for r ah for the assumed neutral atmospheric condition and smaller value for ∆ using T a only) is somewhat moderated in the PM method.However, the impacts can still be large.Figure 7a-c shows scatter plots for three estimates for λE based on different combinations of relaxations (simplifications) of the PM equation versus full application of the PM using Landsat-based T s .Figure 7d-f shows the results for the same relaxations, but based on the simulated T s produced in BAANS.Results were sampled from an agricultural AOI for 18 June 2008.Similarly, Figure 7g-i shows the results for Landsat-based T s , and Figure 7j-l shows the results for simulated T s for the desert AOI on the same date.NARR-based air temperature (T a ) at the 30-m blending height was about 296 K, while simulated T s varied from 295 to 330 K for different land use classes, except water bodies (Figure 6).As discussed earlier, BAANS was able to simulate the relatively lower values for T s in fully covered agricultural land and the higher values in sagebrush desert and grassland areas relatively well, producing values similar to Landsat-based T s (Figure 6a,b).Figure 7 shows that application of the fully-parameterized PM method using simulated T s produced nearly identical trends in λE as compared with those estimated based on Landsat T s in both agricultural and desert AOIs.
The scatter plots of λE in the agricultural AOI (Figure 7a,d) show good agreement in λE produced by relaxations of the PM method and the fully-parameterized PM, due to the similarity between T s and T a for many of the irrigated fields, while in the sagebrush desert AOI, λE was systematically underestimated by the relaxations (Figure 7g,j).The root-mean-square deviation (RMSD) between λE PM_∆ra and λE PM for the Landsat overpass time on 18 June 2008, was about 25.6 W/m 2 and 18.5 W/m 2 for the agricultural AOI using Landsat-based and simulated T s, respectively, using λE PM as the basis.Similarly, for the desert AOI, RMSD was about 32.4 W/m 2 and 32.1 W/m 2 from Landsat-based and simulated T s, respectively.The coefficient of determination (R 2 ) was 0.99 for all four scenarios (Figure 7a,d,g,j), indicating strong correlation among estimates, but with strong biases (error) for some of the relaxations.
Figure 8a,d,g,j presents similar information as for Figure 7, but for the Landsat overpass time on 17 May 2008. Figure 8a,d shows that the impact of the λE PM_∆ra relaxation on 17 May 2008, when wind speed was smaller (1.8 m/s) than on 18 June 2008 (4.4 m/s), was to increase estimates for λE slightly in lower λE areas (having sparse or stressed vegetation) and to slightly decrease λE in higher λE areas in the agricultural AOI.Apart from the reduced wind speed, Z om values were generally smaller for agricultural lands on 17 May, as compared to 18 June, because of being earlier in the growing season.Similarly, λE estimates for the sagebrush desert AOI were less impacted by the relaxations on 17 May 2008 (Figure 8g,j), than for the June date.This is probably because of some moderating effects caused by the combination of assumed neutral aerodynamic resistance (rah _neu ), lower wind speed (u z ) and neutral ∆ approx .The RMSD between λE PM_∆ra and λE PM was 23.4 W/m 2 and 17.5 W/m 2 using Landsat-based and simulated T s , respectively, for the agricultural AOI on 17 May 2008.Similarly, for the desert AOI, RMSD was 1.62 W/m 2 and 2.6 W/m 2 for Landsat-based and simulated T s , respectively, for the same date.The coefficient of determination (R 2 ) is 0.99 for all four scenarios (Figure 8a,d,g,j) on 17 May 2008.
Table 5 compiles summaries for λE PM_∆ra , λE PM_∆ and λE PM_∆RL for different land use class pixels where the percentage change is compared to λE PM based on the simulated T s from BAANS.Characteristics of these pixels are shown in Tables 3 and 4. Simulated T s is used to compare λE, as simulated values for Ts are more widely used, since radiometrically-measured T s is generally not available.Table 5, Column 6 summarizes the impact of ∆ approx and the assumption of a neutral condition (λE PM_∆ra ) on the PM method.The maximum decrease in λE between λE PM_∆ra and λE PM was about 25 percent for desert pixels, and the maximum decrease was about six percent in the agricultural AOI on 18    Use of ∆ approx with stability correction in the PM equation is widely used for λE calculation.Figure 7b,e shows scatter plots of λE estimated using the PM with ∆ approx and stability correction (λE PM_∆ ) versus λE estimated by the fully and appropriately parameterized PM (λE PM ) for an agricultural AOI using Landsat-based and simulated T s , respectively, on 18 June 2008.Similarly, Figure 7h,k shows the scatter plots for the desert AOI for the same date.Overall, the results show that estimates for λE systematically decreased with the use of ∆ approx in the PM method (Figure 7b,e,h,k).For the compared pixels, the maximum decrease in latent heat flux (λE PM_∆ ) was about 22 percent for an agricultural pixel with a low LAI (i.e., 0.063) on 18 June 2008, while a maximum decrease in a desert pixel was about 31 percent on the same date (Table 5, Column 8).These errors are considered to be substantial and generally unacceptable in practice.Scatter plots show a larger deviation in estimated λE for the desert AOI (Figure 7h,k) as compared to the agricultural AOI (Figure 7b,e) on 18 June 2008.On that day also, the RMSD between λE PM_∆ and λE PM was about 17.3 W/m 2 and 19.6 W/m 2 for Landsat-based and simulated T s , respectively, for the agricultural AOI.For the desert AOI, the RMSD was about 7.8 W/m 2 and 11.8 W/m 2 using Landsat-based and simulated T s , respectively, on the same date.
A similar trend was observed for 17 May simulations with smaller λE PM_∆ than λE PM in all scenarios (Figure 8b,e,h,k) and larger deviation in estimated λE for desert AOI.On 17 May of the compared pixels, the maximum decrease was about 33 percent for the same agricultural pixel having a low LAI (i.e., 0.05) (Table 4, 5,103,099 m, 3,807,774 m).Simulated T s was about 323 K, and air temperature (T a ) at 30 m from NARR reanalysis was 297 K for the agricultural pixel (5,103,099 m, 3,807,774 m) at Landsat overpass time on that date.Because differences between Ta and Ts were high, even if the stability corrected rah is used in the PM method, the error in estimate λE was large for this low-vegetation agricultural pixel due to the impact on Δ.The maximum decrease in estimated λE occurred in Land Use 71 (grassland) and was about 39 percent.Again, these impacts are very significant and reflect substantial error in ET estimation.The RMSD between λE PM_∆ and λE PM was about 50.3 W/m 2 and 36.7 W/m 2 for Landsat-based and simulated T s , respectively, for the agricultural AOI on 17 May 2008.For the desert AOI, RMSD was about 41.6 W/m 2 and 41.3 W/m 2 for Landsat-based and simulated T s , respectively, on the same date.For both dates and AOIs, R 2 was about 0.99 (Figures 7b,e,h,k and 8b,e,h,k).The scatter plots between λE PM_∆ and λE PM (Figure 8b,e) show a larger deviation for the agricultural AOI on 17 May 2008, compared to 18 June 2008 (Figure 7b,e).λE was largely underestimated by λE PM_∆ in the desert AOI for both dates (Figures 7h,k and 8h,k).Table 5 also confirms that there was a systematic decrease in λE with the relaxation represented by λE PM_∆ .These results indicate that even if the stability corrected surface energy flux parameters are used in the PM method, the assumption of ∆ approx can significantly decrease λE in sparse vegetation where T s and T a are significantly different.Table 6 shows the compilation of the statistics of the simplifications of the PM method for all AOIs and dates.4.2.3.Impact of Using T a to Estimate R L↑ and ∆ approx with (ψ m , ψ h ) (λE PM_∆RL ) This relaxation of the PM method (λE PM_∆RL ) suffers two distinct types of simplifications, i.e., the use of T a in place of T s for outgoing longwave radiation, R L↑ and for ∆ approx in the PM equation.This relaxation is routinely done in practice.R L↑ is, in actuality, a function of T s and can deviate widely from T a , especially for dry surfaces.R L↑ decreases and net radiation (R n ) increases when T a < T s is used in place of T s in the surface energy balance.Scatter plots of estimated λE for the agricultural AOI show that the differences between λE PM_∆RL and λE PM are small (Figure 7c,f) when the differences between T a and T s are small; although, Figure 7c shows some variation on 18 June 2008.Out of the compared pixels, Table 5 shows the maximum decrease in estimated λE for agricultural pixels was about 14 percent; and 22 percent in the desert pixel on 18 June 2008.The RMSD between λE PM_∆RL and λE PM was about 34.6 W/m 2 and 7.3 W/m 2 using thermal-based and simulated T s , respectively, for the agricultural AOI on 18 June 2008.Similarly, for the desert AOI, RMSD was about 25.7 W/m 2 and 27.3 W/m 2 using thermal-based and simulated T s , respectively, for the same date.The most significant impact with this simplification was observed on 17 May 2008, with about a 45 percent decrease in estimated λE for both agricultural and desert AOIs within the compared pixels (Table 5, Figure 8c,f,i,l).Use of T a significantly decreased ∆, which then decreased λE PM_∆RL on this date .These kinds of relaxations on the PM equation also cause closure problems in the surface energy balance.On 17 May 2008, for the agricultural AOI, the RMSD was about 36.5 W/m 2 and 63 W/m 2 using Landsat-based and simulated T s , respectively.For the desert AOI, the RMSD was about 44.1 W/m 2 and 56.3 W/m 2 using Landsat-based and simulated T s , respectively, on the same date.The R 2 varied from 0.96 to 0.99 for different scenarios (Figure 8c,f,i,l) for agricultural and desert AOIs on 17 May 2008.Figures 7i,l and 8c,f,i,l show that estimated λE systematically decreased in agricultural, grassland and sagebrush desert AOIs with this common relaxation.

Conclusions
This study analyzed various combinations of parameterizations and relaxations on the rigor of parameter estimation in the Penman-Monteith (PM) method, primarily associated with the use of air temperature, T a , in place of surface temperature, T s .Analyses were made over a range of different land use classes and surface conditions.Gridded weather and remote sensing data and techniques were used to drive a newly developed backward averaged accelerated numerical solution (BAANS) of the surface energy balance.These applications provided an opportunity to understand and analyze various complexities in estimating latent heat flux (λE).BAANS is able to calculate and close the surface energy balance without using satellite-based and iteratively-derived surface temperature.Both the thermal values from Landsat and simulated surface temperature were used in the surface energy balance to access the accuracy of the surface energy parameters from BAANS.This study quantitatively assessed the advantage of using separate aerodynamic equations against the combined PM method for calculating λE in the surface energy balance when surface temperature (T s ) is significantly higher than air temperature (T a ).Analyses also emphasize the large errors in estimated λE that can occur from the use of common relaxations in the PM method.Only in fully covered agricultural land where T s and T a are nearly equal are errors caused by the relaxations small.Bulk surface resistance (r s ) computed using the surface energy balance was identical for the aerodynamic and PM methods, indicating the accuracy of the PM method and its equivalency to the full aerodynamic energy balance methods when the fully parameterized PM method is used.r s from both methods generated the expected small values in agricultural land (~0-120 s/m) and larger values for dry surfaces, including sagebrush desert and grasslands (up to r s ~4000 s/m).The assumption of a neutral atmospheric condition and the use of approximate slope of the saturation pressure (∆ approx ) in the PM method both increased and decreased λE in different land use classes and atmospheric conditions.For specifically compared pixels, the estimation error caused by the relaxation to parameters in the PM method was as high as 48 percent.The maximum root-mean-square deviation (RMSD) was about 32 W/m 2 for the desert area of interest (AOI) on 18 June 2008.Out of the compared pixels, the maximum increase in λE was about 22 percent, and the decrease was about 25 percent in different land use classes (17 May and 18 June 2008) with the assumption of the neutral atmospheric condition and the use of ∆ approx in the PM method.The use of ∆ approx with stability-corrected surface energy flux parameters in the PM method decreased λE in all land use classes with a maximum RMSD of about 41 W/m 2 in desert AOI.With this simplification, the maximum decrease in λE was about 39 percent in desert grassland at 17 May 2008.The errors produced by the relaxation (simplification) in outgoing longwave radiation (R L↑ ), Δ and boundary layer stability correction are large and intolerable when T s and T a are significantly different.T s should be used to estimate all of these parameters when applying the PM method, where T s can be determined iteratively using a surface energy balance.As an alternative, a full aerodynamic-surface energy balance method can be used to estimate λE, where the use of the PM method is avoided in total.

Figure 1 .
Figure 1.General inputs to the backward-accelerated numerical solution of the surface energy balance (backward-averaged accelerated numerical solution (BAANS)) model from gridded weather and satellite-based data and the strategy for convergence of surface energy fluxes.NARR, North American Regional Reanalysis; METRIC, Mapping Evapotranspiration at High Resolution using Internalized Calibration.

Figure 2 .
Figure 2. Incoming longwave radiation R L↓ from NARR reanalysis data sets for the 32-km study area in southern Idaho overlaying a Landsat path 39 image on 18 June 2008, overlain by outlines of Idaho counties.The R L↓ legend covers the range 180 W/m 2 -524 W/m 2 for North America on 18 June 2008.

Figure 3 .
Figure 3. Bulk surface resistance in s/m from inverting the full PM equation (r s_PM ) for 18 June 2008.

Figure 4 .
Figure 4. Scatter plots of bulk surface resistance (r s ) from inverting PM equation with full use of surface temperature (T s ) (r s_PM ) versus from inverting the aerodynamic method (r s_aero ) for about 27,000 pixels in agricultural and desert areas of interest (AOI) for 18 June 2008 and 17 May 2008, respectively.(a) Desert AOI on 18 June 2008; (b) agricultural AOI on 18 June 2008; (c) desert AOI on 17 May 2008; (d) agricultural AOI on 17 May 2008.

Figure 6 .
Figure 6.Landsat-based T s (K), simulated T s (K) and differences between Landsat-based T s and T a (K), respectively, on 18 June 2008.

Figure 7 .
Figure 7. Scatter plots of latent heat flux (W/m 2 ) for about 27,000 locations in agricultural and desert areas of interest at Landsat overpass time on 18 June 2008, for different combinations of the relaxation of parameters in the PM method using Landsat (thermal)-based T s and simulated T s .The black symbols represent a higher density of points and the gray symbols a lower density.
June 2008.Estimates of λE increased up to 21 percent for the agricultural AOI on 17 May 2008.The λE PM_∆ra relaxation both increased and decreased estimated λE for different land use classes for different meteorological conditions, but without any distinct pattern.The λE in the desert AOI was underestimated with the relaxation on 18 June 2008, while on 17 May and 18 June 2008, the estimated λE simultaneously increased and decreased for agricultural AOI.

Figure 8
Figure 8 Scatter plots of latent heat flux (W/m 2 ) for about 27,000 locations in agricultural and desert areas of interest at Landsat overpass time on 17 May 2008, for different combinations of the relaxation of parameters in the PM method using Landsat (thermal)-based T s and simulated T s .The black symbols represent a higher density of points and the gray symbols a lower density.

Table 1 .
Parameterization of λE calculations using the Penman-Monteith (PM) method with various relaxations of important parameters and the aerodynamic method.

Table 2 .
Data from NARR reanalysis for the satellite overpass date on 17 May 2008 and 18 June 2008.

Table 3 .
Simulation results for surface energy balance flux parameters on 18 June 2008, using simulated surface temperature for all parameters in the PM.

Table 4 .
Simulation results for surface energy balance flux parameters on 17 May 2008, using simulated surface temperature for all parameters in the PM.

Table 5 .
Summary of simulation results for estimated latent heat flux (λE) and percentage change relative to the full aerodynamic and PM applications (λE aero and λE PM ), based on simulated surface temperature on the 18 June 2008 and 17 May 2008, Landsat overpass date.Impact of Using ∆ approx , but with Stability Correction (ψ m , ψ h ) (λE PM_∆ )

Table 6 .
Statistics for estimated latent heat fluxes (λE) for three combinations of relaxations on Penman-Monteith parameters.