Correction of Eddy Covariance Based Crop ET Considering the Heat Flux Source Area

: Eddy covariance (EC) systems are being used to measure sensible heat (H) and latent heat (LE) ﬂuxes in order to determine crop water use or evapotranspiration (ET). The reliability of EC measurements depends on meeting certain meteorological assumptions; the most important of such are horizontal homogeneity, stationarity, and non-advective conditions. Over heterogeneous surfaces, the spatial context of the measurement must be known in order to properly interpret the magnitude of the heat ﬂux measurement results. Over the past decades, there has been a proliferation of ‘heat ﬂux source area’ (i.e., footprint) modeling studies, but only a few have explored the accuracy of the models over heterogeneous agricultural land. A composite ET estimate was created by using the estimated footprint weights for an EC system in the upwind corner of four ﬁelds and separate ET estimates from each of these ﬁelds. Three analytical footprint models were evaluated by comparing the composite ET to the measured ET. All three models performed consistently well, with an average mean bias error (MBE) of about − 0.03 mm h − 1 ( − 4.4%) and root mean square error (RMSE) of 0.09 mm h − 1 (10.9%). The same three footprint models were then used to adjust the EC-measured ET to account for the fraction of the footprint that extended beyond the ﬁeld of interest. The effectiveness of the footprint adjustment was determined by comparing the adjusted ET estimates with the lysimetric ET measurements from within the same ﬁeld. This correction decreased the absolute hourly ET MBE by 8%, and the RMSE by 1%.


Introduction
Measuring and modeling ET is difficult due to the nature of water vapor transport into the atmosphere. Allen et al. [1] discussed the factors governing ET measurement accuracy for the following methods: 1. Soil water balance, 2. Lysimetry, 3. Energy Balance Bowen Ratio (EBBR), 4. Eddy Covariance (EC), 5. Scintillometry, 6. Sap flow, and 7. Remote sensing and satellite-based modeling. Soil water balance is an affordable and relatively accurate method that can be used for irrigation scheduling. However, for a novice or a person working outside their specialty, the errors in measurement can be 20-70% [1]. Methods 2-7 require expensive instrumentation and expert operators. Remote sensing and satellite-based modeling holds great potential for practical application, since the process covers large areas that can be instantaneously available via the internet. Sap flow methods directly measure the transpiration through a plant, but large errors are introduced when attempting to scale up measurements to a field or regional scale. Methods 2-5 are primarily used in research, and are being used to calibrate and validate remote sensing models. Properly managed lysimeters have the potential of measuring ET with high accuracy, according to Howell et al. [2]. However, these instruments are large, expensive, and only provide the user with a measurement that represents a rather small area. Methods 3 and 5 are micrometeorological approaches that are based on the conservation of energy. Since EC takes measurements above the transpiring canopy, the temperature and scalar concentrations sampled are actually a mixture of upwind point sources. The EC system has

Flux Footprint Modeling
The 'heat flux source area' (i.e., footprint) is the portion of the upwind surface of the EC instrumentation site containing the effective heat flux sources/sinks contributing to a flux or concentration observation at a certain measurement height [3,4]. If the EC system is located within the constant flux layer over an extended homogeneous surface, then the position of the sensor is not an issue. The extent of the footprint is significantly affected by the stability conditions of the atmosphere, and thus there is better energy balance closure under neutral to unstable conditions [5]. However, over a heterogeneous surface (e.g., a patchwork of agricultural land) the location and size of the flux footprint is needed to interpret the measured flux (i.e., to understand the source area of the heat fluxes). As a result, there has been a proliferation of footprint modeling research since 1990. There are four theoretical approaches: analytical models, Lagrangian-stochastic particle dispersion models, large-eddy simulations, and ensemble-averaged closure models [6,7]. The latter three approaches are mathematically complicated and thus resource intensive. Analytical models can be misleading if they are used in conditions that violate the underlying assumptions of each particular model [6,7]. Nevertheless, analytical models can easily be applied to processing code, and can thus be used to filter and correct flux measurements quickly in post-processing, or possibly in real time if they are programmed in the datalogger. The analytical footprint models of Hsieh et al. [8] and Kormann and Meixner [9] have been used as quality control filters for EC data collected over various land covers [10][11][12][13], and have been used to scale and validate remotely-sensed energy balance models [14][15][16]. The Kormann and Meixner model [9] is also employed by the frequently used open source EC data processing software EdiRe [17] and TK3 [18]. The intercomparison of footprint models was the primary method of validating footprint models in the past [8,19,20]. Foken and Leclerc [21] proposed three methods for the 'in situ' validation of models: (1) the use of artificial tracers, (2) the use of natural tracers, and (3) the effect of isolated heterogeneities. If a suitable artificial tracer is chosen, there is the advantage of there not being any other source or sink for the tracer. The disadvantages of artificial tracers are the difficulty of approximating natural field conditions and the expense of setting up these experiments. Natural tracer validation methods are both inexpensive and easier to use. Cooper [22] used water vapor as a natural tracer, and found good agreement between the point flux measurements of ET, lidar-derived water vapor fluxes, and a footprint model. By comparing the measurements from two adjoining surfaces and those of the two surfaces combined, Beyrich et al. [23] found that Lagrangian stochastic simulation [24] better represented the flux contribution from different fields than Schmid's analytical model [25]. Marcolla and Cescatti [26] compared three analytical models by comparing measurements from an alpine meadow before and after cutting to those during an intermittent time when only a portion of the meadow was cut. As a result, they found that the Schuepp et al. [4] model generally overestimates the footprint, and that none of the footprint models considered perform well in very unstable atmospheric conditions.

Motivation and Objectives
Analytical footprint models have been extensively studied for different atmospheric conditions, and have been compared to other more sophisticated models (e.g., Lagrangian stochastic). Only a few studies have explored the accuracy of such models over irrigated cropland surrounded by rainfed crops and/or fallow land. Furthermore, there has not been a study that has explored the correction of flux measurements for a footprint that extends beyond the area of interest. Therefore, the objectives of this study are: (a) to compare the performance of the Schuepp et al. [4], Hsieh et al. [8] and Kormann and Meixner [9] analytical footprint models over irrigated cotton adjacent to dryland cotton; and (b) to develop and evaluate a method for improving the LE flux, measured by an EC system, by accounting for the footprint fraction that extends beyond the irrigated field boundaries.

Site Description
For this study, data from the Bushland Evapotranspiration and Agricultural Remote Sensing Experiment 2008 (BEAREX08) was used [27]. The geographic coordinates of the USDA-ARS Conservation and Production Research Laboratory (CPRL), located at Bushland, TX, USA are 35 • 11 N, 102 • 06 W, and its elevation is 1170 m above mean sea level. The soils in and around Bushland are classified as slowly permeable Pullman clay loam. The major crops in the region are corn, sorghum, winter wheat, and cotton. The wind direction is predominantly from the south/southwest direction. The average precipitation that occurs during the cotton growing season (May-October) is 350 mm [28]. About 600 mm of irrigation, precipitation, and soil water are needed to grow cotton [29]; thus, irrigation needs to provide about 250 mm of timely water for a successful cotton harvest. The typical growing season grass reference ET is 6.0-8.2 mm/day −1 [28]. In addition, the long-term annual microclimatological conditions indicate that the study area is subject to very dry air and strong winds. The growing season averages at Bushland for air temperature and horizontal wind speed are 20 • C and 3.9 m s −1 , respectively [28].

Large Monolith Weighing Lysimeters
Four precision weighing lysimeters [30], 3 × 3 × 2.3 m deep, were used to measure the cotton ET. Each lysimeter contained a monolithic Pullman clay loam soil core. The lysimeters were located at the center of four fields (210 m east-west by 225 m north-south), two (east) of which were irrigated by a linear move system, and two (west) of which were not irrigated. The change in the lysimeter mass was measured by a load cell (SM-50, Interface, Scottsdale, AZ, USA) and recorded by a datalogger (CR7-X, Campbell Scientific Inc., Logan, UT, USA) at a 0.17 Hz frequency. The high frequency load cell signal was averaged for 5 min and composited to 15-min means. The lysimeter was calibrated as explained in Howell et al. [2]. A simple soil water balance using the change in water storage from four neutron probes, irrigation, and precipitation data showed that the northeast (NE) lysimeter had a larger ET than the surrounding field [31]. According to that study, the ET measurements from the NE lysimeter were 18% larger than the surrounding field from Day of Year (DOY) 182 to 220, as determined by a soil water balance using the neutron probe soil water content readings, due to the greater leaf area index (LAI) on the lysimeter than the surrounding field. Therefore, the NE lysimeter ET measurements, from that period, were reduced by 18% in order to ensure that they were representative of the entire field, as indicated by Evett et al. [31]. Pictures of a lysimeter box and an EC system are shown in Figure 1.

Eddy Covariance Energy Balance System
Five EC systems (EC1, EC2, EC3, EC5, and EC8) out of the nine from the BEAREX08 experiment were used to monitor the exchange of heat fluxes at different parts of the CPRL site. The instrument positions are shown in Figure 2. The instrumentation details of all of the systems used are given in Table 1. The time series data consisted of the horizontal (u),

Eddy Covariance Energy Balance System
Five EC systems (EC1, EC2, EC3, EC5, and EC8) out of the nine from the BEAREX08 experiment were used to monitor the exchange of heat fluxes at different parts of the CPRL site. The instrument positions are shown in Figure 2. The instrumentation details of all of the systems used are given in Table 1. The time series data consisted of the horizontal (u), lateral (v), and vertical (w) wind vectors (m s −1 ); calculated sonic temperature (T s , • C); water vapor concentration (H 2 O, g m −3 ); carbon dioxide concentration (CO 2 , mg m −3 ); and atmospheric pressure (P, kPa), all measured at a frequency of 20 Hz. The Campbell Scientific 3-dimensional (3D) sonic anemometer (CSAT3) sensor was oriented toward the predominant wind direction, with an azimuth angle of 225 • from true north for EC8, and 180 • for systems EC1, EC2, EC3, and EC5. Installed within and between the crop rows, about 4 m east from each EC system, were instruments for measuring soil heat flux, soil temperature, and soil volumetric water content. The temperature and water content data were used to calculate the soil heat storage in the layer between the surface and the depth of the soil heat flux plate installation. The soil heat flux plates (SHFP) were installed at 0.08 m depth within and between the crop rows. The soil thermocouple pairs were installed at 0.02 and 0.07 m depths close to the SHFP locations. The soil water content reflectometers (CS616, Campbell Scientific Inc., Logan, UT, USA) were installed slanted at an approximate angle of 13 degrees across the 0.02-0.1 m depth, in order to measure the volumetric soil water content in this depth zone. The water content reflectometers were field cross-calibrated against the water contents reported by a conventional time domain reflectometry (TDR) (TR-100, Dynamax, Inc., Houston, TX, USA) system that used a soil-specific calibration that minimized soil temperature influences on the TDR water content readings [32]. Soil temperature was sensed and recorded during the cross-calibration, and a soil temperature correction was developed for the CS616 data (Personal communication, [33]).

Eddy Covariance Data Processing
The major components of the energy balance at the land surface are net radiation (R n ), soil heat flux (G), sensible heat flux (H), and latent heat flux (LE), all in units of W m −2 , and can be expressed as: where the left side in Equation (1) is defined as the available energy (R n − G), and the right side is defined as the turbulent fluxes (H + LE), with the positive sign convention away from the surface for H and LE and negative for R n and G. The EC method is based on the direct measurement of high-frequency vertical wind (w) and a scalar concentration (c), such as water vapor or air temperature, producing LE through Equation (2) and H through Equation (3), respectively, assuming that the mean vertical velocity is negligible: where q is the air-specific humidity (kg kg −1 ); T is the air temperature ( • C); ρ a is the moist air density (kg m −3 ); C p is the specific heat of dry air at constant pressure (J kg −1 K −1 ); λ is the latent heat of water vaporization (J kg −1 ), which varies T; and the primes (apostrophes) denote the deviation or fluctuation from the mean (e.g., Time series, high-frequency EC data were post-processed with the EdiRe software package [17] following the guidelines described in Lee et al. [34] and Burba and Anderson [35], as summarized in Table 2. Before the covariances were calculated, spikes of six standard deviations from the population mean were removed from the time series. If four or more consecutive data points were detected with values larger than the standard deviation then they were not considered as a spike [36]. The time delay between the CSAT3 and LI-7500 was removed using a cross-correlation analysis [37]. Although the terrain for the site was practically flat, the CSAT3 cannot be perfectly leveled, such that the vertical component (w) is perpendicular to the mean streamline plane. For this reason, the coordinates were rotated using the double rotation method [38]. According to Lee et al. [34], this method is suitable for ideal sites with little slope and fair weather conditions. The effects of density fluctuations induced by heat fluxes on the measurement of the eddy fluxes of water vapor using the LI-7500 were corrected using the procedure outlined by Webb et al. [39]. The spectral loss in the high frequency band due to path-length averaging, Atmosphere 2021, 12, 281 6 of 22 sensor separation, and signal processing was corrected after Moore [40]. The data from the fine wire thermocouple was intermittent due to equipment failure, and thus the sonic temperature (T s ) was used in sensible heat flux calculation. Schotanus et al. [41] recommended correcting T s for crosswind and humidity effects, which are commonly referred to as the heat flux correction. The CSAT3 implements the crosswind correction online, and therefore the heat flux only needed to be corrected for humidity fluctuations. The sonic temperature flux was converted to actual air temperature flux following the method of Schotanus et al. [41]. A dimensionless parameter that characterizes the processes in the surface layer is the atmospheric stability parameter (ζ), Equation (4), which is the ratio of the convective production to the mechanical production of turbulent kinetic energy [42]: where z m is the horizontal wind speed observation height above the zero-plane displacement height (i.e., z m = z − d); g is the gravitational acceleration (9.8 m s −2 ); u* is the friction velocity (m s −1 ); and H, T, ρ a , and C p are as defined above. A positive ζ represents a stable stratification of the atmosphere, a negative ζ represents an unstable stratification, and |ζ | < 0.02 represents a neutral stratification.
The canopy heights (h c , m) and leaf area index (LAI, m 2 m −2 ) for the NE and southeast (SE) fields were measured five times during the study, on the following days of the year (DOY): 171, 182, 200, 210, and 220. Field measurements were taken for the northwest (NW) and southwest (SW) fields three times, on DOYs 200, 210, and 220. The crop emergence for the east and west fields occurred on DOY 150 and 165, respectively. Curves were fitted to the h c vs. DOY and LAI vs. DOY data in order to determine the h c and LAI as functions of the DOY; these are shown in Appendix A. The analytical expression from Raupach [43][44][45] was used to estimate the zero-plane displacement height (d, m): where c d1 is an empirical constant estimated from laboratory and field data to be on the order of 7.5 [44], and Λ is the canopy area index, which for unstressed, growing canopies is equivalent to the LAI. Quality control criteria were set for the wind direction, stationarity, and integral turbulence characteristics. The flow from behind the sensor can be distorted by the instrumentation. Therefore, any data with a wind direction beyond ±90 • of the orientation of the sonic anemometer were excluded from further analysis. The estimates of fluxes via the eddy covariance method are based on simplified forms of the Navier-Stokes equations for certain atmospheric conditions [46]. These conditions are not always met, and therefore must be evaluated. The tests for stationarity and integral turbulence were performed following the methods outlined in Foken and Wichura [47], and Thomas and Foken [48].

Footprint Modeling Methodology
The analytical footprint models proposed by Schuepp et al. [4], Hsieh et al. [8], and Kormann and Meixner [9] were used to determine the footprint weight for each hourly flux per unit source area (i.e., 1 m × 1 m), hereafter referred to as S90, H2000, and KM01, respectively. This was achieved by first rotating into the mean wind direction the respective x and y coordinates of a 450 × 420 grid, with each point representing one square meter (i.e., the size of the entire four-field site) centered at the position of each EC system. The rotated coordinates were then used to calculate the relative footprint contribution of each point on the grid. The footprint contribution of each field was then found by creating a 450 by 420 matrix for each field, which consisted of 1 s and 0 s, with the 1 s representing the area of interest; multiplying each of those respective matrices by the matrix of footprint density weights; and then summing all of the elements of the final matrix. The MATLAB codes, following this procedure for each footprint model, are shown in Appendices B and C.

Schuepp Model
The S90 model views the transpiring vegetative surface as a continuum of upwind line sources, each occupying an infinitesimal strip of width δx. As proposed by Gash [49], the approximate analytical solution by Calder [50] was applied to the basic advection-diffusion equation, assuming neutral stability and uniform wind velocity. After differentiating with respect to z m and then integrating with respect to the upwind distance (x), an equation for the vertical concentration gradient at z m is reached: where k is von Kármán's constant (0.41), Q 0 is the area flux density, and U is the uniform wind velocity (m s −1 ), which is defined as the average wind velocity between the surface and z m (m), assuming a logarithmic profile for u(z): where z 0 is the surface roughness length (m). The relative contribution to the vertical flux at height z m , f(x,z m ), is then obtained by taking the derivative of Equation (6) and multiplying by ku*z m : where φ M is the momentum correction for stability, which is a function of Monin-Obukhov stability parameter ζ [51].

Hsieh Model
The H2000 model is a hybrid approach which fits Calder's analytical solution [50] to the results of the numerical Lagrangian model of Thomson [52]. In the analysis of their results, Hsieh et al. [8] scaled Gash's [49] effective fetch with the Monin-Obukhov stability length (L, m), and accounted for the effects of stability by introducing the similarity parameters D and P, so that: with where D and P are found by the regression of Equation (10) (10) and differentiating with respect to x results in the analytical footprint expression:

Crosswind Function
Both S90 and H2000 are one-dimensional (1D) footprint models (i.e., they are expressed along the mean wind direction (x)). The KM01 model is a two-dimensional (2D) model that expresses the footprint in the mean wind direction and the crosswind direction (y). The S90 and H2000 models were expanded into 2D models, so that the results could be compared to KM01. The diffusion in the lateral direction is commonly assumed to have a Gaussian distribution centered on the mean wind direction: where D y is the crosswind concentration distribution function, and σ y is the standard deviation of the lateral spread of the source area, which can be related to the plume travel time, x/ū p , and the standard deviation of lateral wind fluctuations, σ v as σ y ≈ σ v ·x/ū p [9,16,53,54]. The Gaussian distribution was combined with H2000 and S90 to expand them into 2D footprint models:

Kormann and Meixner Model
The KM01 model uses the solution of the resulting two-dimensional advectiondiffusion equation for the power law profiles of the mean wind velocity and the eddy diffusivity. The mathematical framework is a stationary gradient diffusion formulation with height-independent crosswind dispersion. The profile parameters are determined by equating the power law to Monin-Obukhov similarity profiles. The complete mathematical description is given by Kormann and Meixner [9]. However, a condensed set of formulas similar to that presented by Neftel et al. [55] is presented here. First, the profiles of the horizontal wind speed and the vertical eddy diffusivity are described by power law relationships: The proportionality constants α u and α κ , and the exponents m and n are determined by comparing Equation (15) and (16) to the respective Monin-Obukhov similarity profiles at height z m : where φ H is the dimensionless gradient function of the heat profile defined by Dyer [51] as: Based on these quantities, the shape factor r and the constant µ of Equations (22) and (23) are defined as: Finally, the parameters A-E, shown below, are related to the above quantities: The 3D flux footprint can then be expressed as: The first half of the equation is the Gaussian crosswind distribution, and the second half is the crosswind-integrated longitudinal distribution.

Footprint Validation Procedure
The premise of the validation of footprint models using water vapor as a natural tracer is that if the footprint model correctly estimates the footprint weight for each element of the heat flux source area, then those elemental footprint weights represent a fraction of the EC-derived ET. If, then, the ET for each element of the heat flux source area can be accurately measured by another method, the sum of the products of the footprint weights and their respective EC-based ET rates should equal that of the independently-measured ET value. Following that premise, the cumulative footprint within each of the fields was used to calculate a composite ET (ET composite ) which was then compared to the measured ET value at EC8: where F is the cumulative footprint fraction for each field (the subscripts indicate the field), and ET EC1 , ET EC3 , ET EC5 , and ET EC2 are the EC-based ET (mm h −1 ) from EC1, EC3, EC5, and EC2, respectively. The underlying assumption is that the surface conditions within each field are homogeneous. In order to ensure that the EC-based ET values from each field were representative of the field, an infield cumulative footprint fraction limit was set. The selection of the footprint fraction limit was based on the optimization of the exclusion of the EC data that had significant influence from areas outside of the field, while still retaining enough data for the analysis. A minimum 80% of the footprint for EC1, EC3, EC5, and EC2 must come from the NE, NW, SW, and SE fields, respectively.

ET Correction Using Footprint Fractions
When the footprint for a heat flux measurement extends beyond the field of interest, the flux (measured by the EC system) is then influenced by the surrounding area heat fluxes. During the summer of 2008, the east fields were irrigated, and as a result the ET from those fields was greater than that from the west fields. Therefore, any contribution to the flux from the west fields would be cause for an underestimation of the EC-based ET for the east fields. Using the same footprint limit as in the validation procedure, the ET at EC8 was corrected for the footprint fraction extending beyond the NE field using the model shown below: where ET FC is the corrected ET (at EC8 location), considering footprint weights, and dET i is the difference in ET between the ith EC and EC8 systems: The contribution from adjacent fields to the ET measured by EC8 is removed by the terms dET EC3 , dET EC5 , and dET EC2 in Equation (32), while the latter part of the equation removes any contribution to the ET at EC8 that is not accounted for in the adjacent NW, SW, and SE fields, based on the assumption that there is little to no ET from the area extending beyond the four fields, and thus dET = ET EC8 -0. The effectiveness of the adjustment was determined by comparing ET FC to the ET measurements from the NE lysimeter.

Statistical Analysis
The mean bias error (MBE), the root mean square error (RMSE), and the linear regression analysis based on the least squares method for the comparison of a fitted equation slope and intercept were used for the comparison of ET values [56,57]: where N is the number of pairs compared, P i is the predicted or corrected value, and O i is the observed value.

Surface Roughness
Because the crop height and LAI were derived from five spot measurements throughout the study period, it is necessary to evaluate the spatial homogeneity of each field. An underlying assumption for all three of the analytical footprint models presented here is that the surface is spatially homogenous for the extent (fetch) of the footprint. This condition is rarely met when measuring fluxes over agricultural lands due to the typical patchwork of fields, each with a different crop, surface roughness, and water availability. Although it is the objective of this study to ascertain the effectiveness of analytical footprint models over such terrain, each field within this study needs to be reasonably homogeneous in order to properly estimate the composite ET. There are many factors that contribute to the spatial variability of a field (e.g., soil type, irrigation and/or rainfall spatial uniformity, topography, soil fertility, agro-chemical spatial application uniformity, plant germination, etc.). However, a good indicator of the combined effect of these conditions is the vegetation biomass status. Using remote sensing data, the amount of vegetation cover can be determined via the normalized difference vegetation index (NDVI). The NDVI is calculated using visible (VIS) red (R) and near-infrared (NIR) light reflected by the vegetation, and it varies from −1 to 1 [58]. The pigment in plant leaves, chlorophyll, strongly absorbs visible light (in the bandwidth 0.4 to 0.7 µm of the electro-magnetic spectrum) for use in photosynthesis. The cell structure of the leaves, on the other hand, strongly reflects near-infrared light (from 0.7 to 1.1 µm). The more leaves a plant has, the more these wavelengths of light are affected, respectively. The USA National Aeronautics and Space Administration's (NASA) Landsat 5 thematic mapper satellite operates on a 16-day acquisition schedule, in which an image will be captured of any given surface every 16 days. The satellite reflectance images of the site for the stated study period were obtained using the online Earth Explorer tool [59], and NDVI calculations were performed. Three of the nine images covering the crop growth period could not be used due to excessive cloud cover. The remaining six images are shown in Figure 3 (spatial pixel size of 30 m × 30 m). The vegetative cover was uniform over all of the fields, with a mean NDVI of 0.24 and a standard deviation (SD) of 0-0.01 (0.0-4.7%) during the period between DOYs 155 and 171 ( Figure 4). By DOY 187, the east fields began to show greater vegetation coverage, but the SD of the NDVI within each field remained small, at 0.01 (2.5%) and 0.02 (4.0%) for the west and east fields, respectively. The variability in the NDVI (and therefore in the vegetation cover and surface roughness) within each field increased significantly by DOY 203, with the north fields showing the greatest relative variability, with SD of 0.05 (10.4%) and 0.04 (11.3%) for the NE and NW fields, respectively. The variability in the surface vegetation in the NW field was mostly due to poor seed germination in this field, possibly due to chemical residues from the previous crop.
Atmosphere 2021, 12, x FOR PEER REVIEW 12 of 23 tion is rarely met when measuring fluxes over agricultural lands due to the typical patchwork of fields, each with a different crop, surface roughness, and water availability. Although it is the objective of this study to ascertain the effectiveness of analytical footprint models over such terrain, each field within this study needs to be reasonably homogeneous in order to properly estimate the composite ET. There are many factors that contribute to the spatial variability of a field (e.g., soil type, irrigation and/or rainfall spatial uniformity, topography, soil fertility, agro-chemical spatial application uniformity, plant germination, etc.). However, a good indicator of the combined effect of these conditions is the vegetation biomass status. Using remote sensing data, the amount of vegetation cover can be determined via the normalized difference vegetation index (NDVI). The NDVI is calculated using visible (VIS) red (R) and near-infrared (NIR) light reflected by the vegetation, and it varies from -1 to 1 [58]. The pigment in plant leaves, chlorophyll, strongly absorbs visible light (in the bandwidth 0.4 to 0.7 μm of the electro-magnetic spectrum) for use in photosynthesis. The cell structure of the leaves, on the other hand, strongly reflects near-infrared light (from 0.7 to 1.1 μm). The more leaves a plant has, the more these wavelengths of light are affected, respectively. The USA National Aeronautics and Space Administration's (NASA) Landsat 5 thematic mapper satellite operates on a 16-day acquisition schedule, in which an image will be captured of any given surface every 16 days. The satellite reflectance images of the site for the stated study period were obtained using the online Earth Explorer tool [59], and NDVI calculations were performed. Three of the nine images covering the crop growth period could not be used due to excessive cloud cover. The remaining six images are shown in Figure 3 (spatial pixel size of 30 m × 30 m). The vegetative cover was uniform over all of the fields, with a mean NDVI of 0.24 and a standard deviation (SD) of 0-0.01 (0.0-4.7%) during the period between DOYs 155 and 171 (Figure 4). By DOY 187, the east fields began to show greater vegetation coverage, but the SD of the NDVI within each field remained small, at 0.01 (2.5%) and 0.02 (4.0%) for the west and east fields, respectively. The variability in the NDVI (and therefore in the vegetation cover and surface roughness) within each field increased significantly by DOY 203, with the north fields showing the greatest relative variability, with SD of 0.05 (10.4%) and 0.04 (11.3%) for the NE and NW fields, respectively. The variability in the surface vegetation in the NW field was mostly due to poor seed germination in this field, possibly due to chemical residues from the previous crop.

Footprint Validation
The statistical results of the comparison of ETcomposite to ET from EC8 are shown in Table 3. All three footprint models showed large errors. The cause for such large errors was the inter-instrument variability of 26.5±43.0 W m -2 and -87.6±103.3 W m -2 for the NE and SE EC systems (EC1 and EC2), respectively [60]. An alternative method was explored in which the lysimeter data from the east fields (NE and SE) was used, instead of data from EC1 and EC2, to calculate ETcomposite. The SE lysimeter ET was representative of the field ET for the SE field, as was the NE lysimeter after correction, as shown by Evett et al. [31]. Since the variability between the sensors was found to be significant, the ETcomposite was recalculated using the data from the east lysimeters (NE and SE), EC3, and EC5. Due to the variability in the meteorological and surface conditions with respect to time, the ETcomposite was evaluated for two separate periods of time. The statistical results of the comparison of ETcomposite calculated using a combination of lysimeter and EC data to the ET from EC8 are shown in Table 4. During the initial growth stage, the surface roughness was uniform but small, with an NDVI less than 0.27 and 0.32 for the west and east

Footprint Validation
The statistical results of the comparison of ET composite to ET from EC8 are shown in Table 3. All three footprint models showed large errors. The cause for such large errors was the inter-instrument variability of 26.5±43.0 W m −2 and −87.6±103.3 W m −2 for the NE and SE EC systems (EC1 and EC2), respectively [60]. An alternative method was explored in which the lysimeter data from the east fields (NE and SE) was used, instead of data from EC1 and EC2, to calculate ET composite . The SE lysimeter ET was representative of the field ET for the SE field, as was the NE lysimeter after correction, as shown by Evett et al. [31]. Since the variability between the sensors was found to be significant, the ET composite was recalculated using the data from the east lysimeters (NE and SE), EC3, and EC5. Due to the variability in the meteorological and surface conditions with respect to time, the ET composite was evaluated for two separate periods of time. The statistical results of the comparison of ET composite calculated using a combination of lysimeter and EC data to the ET from EC8 are shown in Table 4. During the initial growth stage, the surface roughness was uniform but small, with an NDVI less than 0.27 and 0.32 for the west and east fields, respectively. The smoother surface and high average wind velocity of 5 m s −1 caused the footprints to extend farther upwind than those for later in the study, when the surface was rougher and the wind velocity calmer, at an average of 3 m s −1 (Table 5). The minimum cumulative footprint limit caused more data to be filtered out during this early period, which was the reason for there being less data for it. Only four data points could be obtained for the S90 model for the early growth stage, which is not enough data to draw a good conclusion on its performance during this period of time. However, both H2000 and KM01 showed larger errors during the initial growth stage than they did later in the growing season. As the surface roughness increased with the growth of the crop, the discrepancy between the three models and between ET composite and ET from EC8 narrowed. During the latter period of the study, the MBE was −0.03 to −0.04 mm h −1 (−3.30% to −4.76%) and the RMSE was 0.10 mm h −1 (10.19% to 11.39%). For the entire study period, all three models performed well, with the S90 and H2000 performing slightly better than the KM01.  EC8  69  11  1  3  85  EC3  4  1  74  5  83  EC5  0  3  0  76  79  DOY 195-228  EC8  80  6  1  1  88  EC3  2  1  84  3  90  EC5  0  1  0  85  87  DOY 158-228  EC8  77  7  1  2  87  EC3  3  1  80  4  87  EC5  0  2  0  82  84 During the early growth stage, the surface roughness is very small. The KM01 uses wind velocity, u, instead of the roughness length, z 0 [55]. Under ideal conditions, the use of u and z 0 are equivalent, as they are related. The advantage of using u is that it is measured in situ, along with the other input parameters, whereas the z 0 is derived using the Monin-Obukhov similarity theory or estimated using empirical relationships with the canopy height. The size of the KM01 footprint is heavily dependent on the ratio of u/u*, which is interpreted as the relative strength of the horizontal advection vs. vertical diffusion [55]. As shown in Figure 5, the u/u* was higher during the early growth stage, decreasing steadily until about DOY 200, where it mostly leveled off. The higher u/u* is why the footprint extended farther upwind during the early growth stage, which is the reason why less of the footprint was within the field of interest, as shown in Table 5.
Monin-Obukhov similarity theory or estimated using empirical relationships with the canopy height. The size of the KM01 footprint is heavily dependent on the ratio of u/u*, which is interpreted as the relative strength of the horizontal advection vs. vertical diffusion [55]. As shown in Figure 5, the u/u* was higher during the early growth stage, decreasing steadily until about DOY 200, where it mostly leveled off. The higher u/u* is why the footprint extended farther upwind during the early growth stage, which is the reason why less of the footprint was within the field of interest, as shown in Table 5. The difference in each footprint model under both stable and unstable atmospheric conditions is shown in Table 6. For the stable condition, the S90 performed better, with an MBE of -0.03 mm h -1 (-3.16%) and an RMSE of 0.12 mm h -1 (12.75%). However, Figures 6  and 7 show that the S90 and KM01 models yielded similar shapes, whereas the H2000 model has a wider and longer overall footprint, with a smaller concentration near the tower. For the unstable condition, all three models performed similarly, and considerably better than the stable condition, with the RMSE ranging from 0.05 to 0.07 mm h -1 (7.20 to 9.13%), which is consistent with the shapes of the footprint illustrated in Figures 8 and 9. However, it should be noted that all of the models did show strong influence close to the tower, especially for the unstable atmospheric condition. Because a more accurate estimate of the ET was obtained using the lysimeters in the east fields compared to the west fields in which the EC systems were used, the ETcomposite for the unstable condition should be more accurate as well.  The difference in each footprint model under both stable and unstable atmospheric conditions is shown in Table 6. For the stable condition, the S90 performed better, with an MBE of −0.03 mm h −1 (−3.16%) and an RMSE of 0.12 mm h −1 (12.75%). However, Figures 6 and 7 show that the S90 and KM01 models yielded similar shapes, whereas the H2000 model has a wider and longer overall footprint, with a smaller concentration near the tower. For the unstable condition, all three models performed similarly, and considerably better than the stable condition, with the RMSE ranging from 0.05 to 0.07 mm h −1 (7.20 to 9.13%), which is consistent with the shapes of the footprint illustrated in Figures 8 and 9. However, it should be noted that all of the models did show strong influence close to the tower, especially for the unstable atmospheric condition. Because a more accurate estimate of the ET was obtained using the lysimeters in the east fields compared to the west fields in which the EC systems were used, the ET composite for the unstable condition should be more accurate as well.

Correction using Footprint Fractions Evaluation
On average, 16-30% of the footprint extended beyond the east fields during the study period. The EC ET footprint correction, following the procedure outlined in Equations (31) and (32), was systematically applied to the fluxes of good data quality and compared to the lysimetric ET rates. As shown in Table 7, the footprint corrected EC ET resulted in an MBE of 0.01 mm h -1 (0.63%), an RMSE of 0.10 mm h -1 (11.63%), and a slope of 0.99 for KM01. The correction tended to overcorrect during the early growth stage. Each model performed well in correcting EC ET, with the S90 and KM01 models performing slightly better than H2000. The results show that the area surrounding the field had a significant influence on the heat fluxes measured at site EC8. The EC ET corrected for the footprint beyond the area of interest results in a reliable estimate of ET.

Correction using Footprint Fractions Evaluation
On average, 16-30% of the footprint extended beyond the east fields during the study period. The EC ET footprint correction, following the procedure outlined in Equations (31) and (32), was systematically applied to the fluxes of good data quality and compared to the lysimetric ET rates. As shown in Table 7, the footprint corrected EC ET resulted in an MBE of 0.01 mm h −1 (0.63%), an RMSE of 0.10 mm h −1 (11.63%), and a slope of 0.99 for KM01. The correction tended to overcorrect during the early growth stage. Each model performed well in correcting EC ET, with the S90 and KM01 models performing slightly better than H2000. The results show that the area surrounding the field had a significant influence on the heat fluxes measured at site EC8. The EC ET corrected for the footprint beyond the area of interest results in a reliable estimate of ET.

Conclusions
The validation of footprint models is only as accurate as the estimate of the natural tracer (e.g., water vapor flux, ET) for each ground cover condition within the footprint. In this study, there was considerable variability in the estimate of ET between the EC systems deployed, which brought the validation procedure into question because the errors in the footprint estimate and instrument variability could not be differentiated for the EC1 and EC2 systems. The ET composite estimated using a combination of the EC and lysimeter-based ET estimates showed evidence that all three tested 'heat flux source area' models seem to accurately estimate the footprint. The footprints extended farther upwind early in the growing season due to the higher relative strength of advection vs. vertical diffusion, which is related to the smoother surface when the crop is small. The somewhat larger errors between the composite ET and measured ET at EC8 during this early growth stage were probably due to the inaccurate assumption that there was no ET contribution from the fallow land to the west of the study site. The correction of EC-based ET considering the footprint fraction that extends beyond the field of interest decreased the absolute MBE by about 8%, and the RMSE by about 1%. Each of the three footprint models yielded a good estimate of the footprint over the highly advective and heterogeneous agricultural land of the Texas Panhandle. The H2000 model gave slightly more consistent results across all of the growth stages and atmospheric stability conditions. Depending on the surface roughness and field dimensions, the footprint could be primarily confined within the boundaries of the field. For a cotton field, once the crop height and NDVI reaches 0.25 m and 0.3, respectively, then a fetch of 350 m is sufficient to confine more than 80% of the footprint within the field for an EC system installed 2.5 m above the ground.
As a result from this study, it is recommended that care be taken when measuring ET during the early growth stage of a crop, when the surface is smoother. An option during this period of time would be to deploy the EC system at a minimum height of two meters, which increases the influence of the surface roughness on the turbulence. The system height could then be raised after the initial growth stage. The S90, H2000, or KM01 footprint models should primarily be used as a tool to interpret the source area contribution to heat fluxes to an EC system, and thus as a tool to verify the validity or representativeness of the data. In addition, the correction of EC ET using the proposed model, when the heat flux source area extends beyond the field of interest boundaries, is an acceptable method to obtain a more accurate estimate of the eddy covariance-derived hourly ET.  Acknowledgments: The authors are grateful for the data provided by Terry A. Howell, William P. Kustas, John H. Prueger, and Steven R. Evett from the USDA ARS. The authors would like to thank the USDA-ARS and CSU CEE as well as CSU -Extension for providing the funding for this research. In addition, we are grateful to Jay Ham (CSU) for his assistance, and to many technicians from the USDA-ARS, Conservation and Production Research Laboratory (CPRL), Bushland, Texas, who helped with the instrumentation and data manipulation, etc.

Conflicts of Interest:
The authors declare no conflict of interest.