Improved Albedo Estimates Implemented in the METRIC Model for Modeling Energy Balance Fluxes and Evapotranspiration over Agricultural and Natural Areas in the Brazilian Cerrado

In this study we assessed METRIC (Mapping Evapotranspiration at high Resolution with Internalized Calibration) model performance to estimate energy balance fluxes and evapotranspiration (ET) in two heterogeneous landscapes in the Brazilian Cerrado, including fluxes and ET in both agricultural and natural vegetation. The estimates were evaluated by comparing them to flux tower data collected over sugarcane (USR site), woody savanna (PDG site) and stricto-sensu savanna (RECOR site) areas. The selection of the study years (2005–2007 for USR/PDG sites and 2011–2015 for RECOR site) was based on the availability of meteorological data (to be used as inputs in METRIC) and of flux tower data for energy balance fluxes and ET comparisons. The broadband albedo submodel was adjusted in order to improve Net Radiation estimates. For this adjustment, we applied at-surface solar radiation simulations obtained from the SMARTS2 model under different conditions of land elevation, precipitable water content and solar angles. We also tested the equivalence between the measured crop coefficient (Kc_ec) and the reference evapotranspiration fraction (ETrF or F), seeking to extrapolate from instantaneous to daily values of actual evapotranspiration (ETa). Surface albedo was underestimated by 10% at the USR site (showing a better performance for full crop coverage), by 15% at the PDG site (following the woody savanna dynamics pattern through dry and wet seasons) and was overestimated by 21% at the RECOR site. METRIC was effective in simulating the spatial and temporal variability of energy balance fluxes and ET over agricultural and natural vegetation in the Brazilian Cerrado, with errors within those reported in the literature. Net radiation (Rn) presented consistent results (coefficient of determination (R2) > 0.94) but it was overestimated by 8% and 9% in sugarcane and woody savanna, respectively. METRIC-derived ET estimates showed an agreement with ground data at USR and PDG sites (R2 > 0.88, root mean square error (RMSE) up to 0.87 mm day−1), but at the RECOR site, ET was overestimated by 14% (R2 = 0.96, mean absolute error (MAE) = 0.62 mm.day−1 and RMSE = 0.75 mm day−1). Surface energy balance fluxes and ET were marked by seasonality, with direct dependence on available energy, rainfall distribution, soil moisture and other parameters like albedo and NDVI.


Introduction
Brazilian tropical biomes, mainly represented by the Amazon rainforest and the Cerrado (Brazilian savanna), have been the focus of global attention due to their great biological diversity as well as the impact of deforestation on climate, hydrology, and biogeochemical cycles at regional and global scales [1].The Brazilian Cerrado consists of complex vegetation which is ecologically and physiognomically related to other savannas in tropical America, Africa and Australia [2].Its area includes much of central Brazil, part of northeastern Paraguay and eastern Bolivia, making it not only the second largest Brazilian biome (2,036,448 km 2 -approximately 24% of the national territory), but also one of the richest regions in terms of biodiversity on the planet [3].Besides its importance for national agriculture, the Cerrado plays a vital role in the dynamics of water resources in Brazil, including parts of 12 hydrographic regions, in addition to housing some of the largest hydroelectric power plants in the country (accounting for 80% of Brazil's national electrical energy production) [4].
This biome has been showing deforestation rates as high as those for the Amazon region, despite possessing an area that is half of its size.IBAMA's Biome Monitoring and Deforestation Program estimated that in 2010 the deforested area in the Cerrado was 6469 km 2 , while in the Legal Amazon it was 7000 km 2 .Until 2010, the Cerrado's deforested area corresponded to 47% of the biome total area, which demonstrates that in some of its central-south regions, native vegetation has almost disappeared altogether.In São Paulo State, more than 90% of the original Cerrado vegetation has been cleared, as opposed to 16.6% in Piauí State.The states with the largest converted areas are Goiás, Minas Gerais, and Mato Grosso do Sul [3,5,6].For a long time this region was seen as inadequate for agriculture due to the low fertility of its soils, but advances in productive and technological systems, governmental incentives, favorable topographic conditions for mechanization, and low land prices have contributed to turning the Cerrado into a region of intense agricultural activity [7], especially for soybean, corn, sugarcane, beans, cotton, and coffee crops [2,8].Thus, agricultural expansion in the Cerrado is the main reason why Brazilian agriculture occupies an internationally scrutinized position, with a biome that has been considered one of the main global hotspots in the expansion of agricultural crops in the last 50 years [9].
However, land cover and land use changes, triggered by the expansion of agricultural activities, promote the emission of greenhouse gases, causing impacts on water and energy balances [7,8,10], and have the potential to affect ecosystems and diverse economic sectors, such as energy production and water supply.Below average rainfall in some regions of Brazil already requires attention and the careful monitoring of production and food supply throughout the year.Combined with low precipitation, elevated temperatures increase evapotranspiration, which reduces water availability for agricultural crops.Therefore, in the current context of a water crisis, any attempt to improve water use efficiency should be based on reliable ET estimates, including evaporation of soil water or free surfaces and vegetation transpiration [11,12].
Evapotranspiration is the total amount of water transferred from soil, free surfaces, wet vegetation, and plant transpiration to the atmosphere.It is the second largest component of the water cycle (after precipitation), as it returns more than 60% of the precipitable water to the atmosphere, and uses more than half of the solar energy absorbed by the Earth's surface [13][14][15], relating directly to energy, water availability and climate.In addition to its importance to the hydrological cycle, ET has a relevant participation in global climate dynamics and processes of the primary productivity of terrestrial ecosystems [16][17][18].This variable is also related to land cover, as it is controlled by the available solar energy, aerodynamic conditions of the atmosphere, water availability in the soil, and by vegetation characteristics such as phenological stage, canopy geometry and structure/depth of the root system [19].
The spatial and temporal monitoring of ET represents a great challenge to the understanding of the energy and the hydrological partitioning between the surface and the atmosphere in different biomes, especially in tropical areas, where the largest sources of ET are located.In addition, accurate ET estimates are needed to model water resources, assess water balance, assist in irrigation management at the local and regional scale, and monitor water availability and consumption [1,20,21].
In agriculture, the actual crop ET (ETa) is traditionally calculated by multiplying the crop coefficient (Kc) by the reference evapotranspiration (ETr), estimated using field weather-based data and the Penman-Monteith (PM) equation [22,23].ETa consists on the amount of water transferred to the atmosphere under the actual conditions of atmospheric demand, soil moisture, and crop conditions [24].Reference ET is defined as the ET of an hypothetical reference surface, such as alfalfa (ETr) or grass (ETo), with a certain height and surface aerodynamic resistance, under well-watered conditions and representing an extension of at least 100 m with the same or similar vegetation [25].This approach is widely used in agriculture because of its simplicity and low data requirements [26].
Other ways to accurately determine ET in situ are lysimeters and towers equipped with Bowen Ratio or eddy covariance (EC) systems [27].These methods are widely used in many biomes, including deserts, savannas, pastures, forests, and agricultural areas, but they present a high cost of installation and maintenance.Another limitation of these methods is the impossibility of extending the analysis to large areas due to the heterogeneity of their surfaces and the dynamic nature of heat transfer processes [11,28,29].In this sense, remote sensing has the advantage of allowing reliable estimations of spatialized ET for various types of vegetation cover, with low cost and adequate repeatability for monitoring agricultural crops [29].Moreover, the continuous and spatialized quantification of surface physical parameters allows the improvement of agrometeorological models as well as weather and climate predictions [30].
Among the various models developed to estimate ET with remote sensing data, remote sensing energy balance (RSEB) models are currently considered to be the best method for monitoring water use in agricultural crops at local and regional scale [11,17,23,31], requiring surface reflectance and thermal images, as well as complementary surface data.METRIC (Mapping EvapoTranspiration with high Resolution and Internalized Calibration) [32] is a one-source model derived from the Surface Energy Balance (SEBAL) model [33] for applications in irrigated agricultural areas.METRIC utilizes ground-based ETr data acquired at a region of interest to inversely and internally calibrate surface energy balance, which facilitates determining the extreme conditions of heat exchange over agricultural areas (the anchor pixels), aiming to calibrate the sensible heat flux at the pixel-level [23].
Even though latent heat flux is calculated as the residual term of the energy balance, its accuracy depends on the quality of the modeling for other terms and biophysical parameters such as albedo, leaf area index, and momentum roughness length [23].Thus, in this study we evaluated the potential of the METRIC model to estimate the energy balance fluxes and ET over heterogeneous Cerrado regions, and to represent the spatial and temporal pattern of ET in agriculture and natural land covers.To feed METRIC, we applied MODIS/Terra satellite products and surface meteorological data as input.To validate our estimates, measurements from three flux towers located at (i) a sugarcane field, (ii) a woody savanna and (iii) in a stricto-sensu area were employed.As an improvement in this study, we proposed a modification to the Net Radiation estimation by adjusting the broadband surface albedo submodel, using the SMARTS2 model [34] and the method applied by [35,36], which integrate band reflectances according to the intensities of at-surface incoming solar radiation within the domain of the band.

Study Area
This study was conducted for two heterogeneous areas of interest (AOIs) in the Brazilian Cerrado (Figure 1), with simultaneous presence of agriculture and natural vegetation.Each AOI evaluated consists of a 2 • × 2 • plot centralized in the position of the flux towers (described in Section 2.2).For the Cerrado biome, the climate presents well defined dry and humid seasons, which imposes some water stress on vegetation during part of the year.The dry season can last from 3 to 5 consecutive months (between May and September), in which most of the Cerrado vegetation species flower, renew their leaves, and spread their seeds before they germinate at the beginning of the wet season (which lasts from October to April) [8,37,38].In general, soils are deep and have good drainage, but low fertility.ET also has a strong seasonality, varying between 1 mm day −1 in the dry season to 6 mm day −1 in the wet season [39].In agricultural areas, ET varies according to the type of crop and its cultivation cycle, with significant variations due to vegetation structure that may range from bare soil to a fully developed canopy cover [1].season to 6 mm day −1 in the wet season [39].In agricultural areas, ET varies according to the type of crop and its cultivation cycle, with significant variations due to vegetation structure that may range from bare soil to a fully developed canopy cover [1].The AOI in São Paulo State (Figure 1a) is located in the Rio Grande basin, characterized by remnant Cerrado vegetation and intense agricultural production, with large areas of sugarcane, citrus, and eucalyptus [1].This region was constituted by a mosaic of Cerrado and Atlantic Forest areas until the end of the XIX century.Atlantic Forest areas were converted mainly into agriculture in the early twentieth century, while Cerrado areas were used for timber and livestock extraction [8].With the increasing demand for biofuels since the 1970's, Cerrado areas were converted into agricultural crops such as sugarcane.Since 2003, with the increase in demand for biofuels due to the insertion of flex-fuel cars in Brazil, sugarcane has begun to expand intensely in this region, as well as in other regions of São Paulo State.The notable increase in its planted area, especially after 2007, was mainly due to the conversion of agricultural areas and pastures [3,40,41].The experimental sites are located at the center of this AOI and mounted on a natural Cerrado reserve (PDG site) and on a sugarcane plantation (USR site), described below.The climate is temperate with dry winters and warm summers (Cwb) [42].Around the experimental sites, the average annual rainfall (from 1971 to 2007) and the standard deviation correspond to 1517 ± 274 mm, with a maximum in December (274 ± 97 mm) and a minimum in July and August (27 ± 34 mm).The annual mean temperature is 22 °C, ranging from 25 °C (January) to 19 °C (July) [43].
The second AOI (Figure 1b) is located in the region defined as Central Plateau-Central Brazil, characterized by surface of old plains, covering the Federal District (DF) territory and the eastern part of Goiás State.This AOI includes areas of preserved Cerrado (Brasília National Park (PNB), Águas Emendadas Ecological Station (ESECAE), and Brazilian Institute of Geography and Statistics (IBGE) Ecological Reserve) as well as areas of intense agricultural production, which connect two of the largest hydrographic basins of Brazil: Tocantins-Araguaia and the Prata River Hydrographic Basins [44].The AOI climate is marked by a distinct dry season (May to September), and is classified as Tropical Savanna climate (Aw) [42].According to Reference [44], the Federal District region presents The AOI in São Paulo State (Figure 1a) is located in the Rio Grande basin, characterized by remnant Cerrado vegetation and intense agricultural production, with large areas of sugarcane, citrus, and eucalyptus [1].This region was constituted by a mosaic of Cerrado and Atlantic Forest areas until the end of the XIX century.Atlantic Forest areas were converted mainly into agriculture in the early twentieth century, while Cerrado areas were used for timber and livestock extraction [8].With the increasing demand for biofuels since the 1970's, Cerrado areas were converted into agricultural crops such as sugarcane.Since 2003, with the increase in demand for biofuels due to the insertion of flex-fuel cars in Brazil, sugarcane has begun to expand intensely in this region, as well as in other regions of São Paulo State.The notable increase in its planted area, especially after 2007, was mainly due to the conversion of agricultural areas and pastures [3,40,41].The experimental sites are located at the center of this AOI and mounted on a natural Cerrado reserve (PDG site) and on a sugarcane plantation (USR site), described below.The climate is temperate with dry winters and warm summers (Cwb) [42].Around the experimental sites, the average annual rainfall (from 1971 to 2007) and the standard deviation correspond to 1517 ± 274 mm, with a maximum in December (274 ± 97 mm) and a minimum in July and August (27 ± 34 mm).The annual mean temperature is 22 • C, ranging from 25 • C (January) to 19 • C (July) [43].
The second AOI (Figure 1b) is located in the region defined as Central Plateau-Central Brazil, characterized by surface of old plains, covering the Federal District (DF) territory and the eastern part of Goiás State.This AOI includes areas of preserved Cerrado (Brasília National Park (PNB), Águas Emendadas Ecological Station (ESECAE), and Brazilian Institute of Geography and Statistics (IBGE) Ecological Reserve) as well as areas of intense agricultural production, which connect two of the largest hydrographic basins of Brazil: Tocantins-Araguaia and the Prata River Hydrographic Basins [44].The AOI climate is marked by a distinct dry season (May to September), and is classified as Tropical Savanna climate (Aw) [42].According to Reference [44], the Federal District region presents periodic and regular water stress conditions, with approximately 47% of the precipitation falling between December and March.The mean temperature varies between 18 • C (July) and 23 • C (October), and the average annual rainfall (from 1981 to 2010) in this region is 1417 mm year −1 , with a monthly average below 10 mm month −1 in July and a maximum of 243 mm month −1 in December [45].
The flux tower is located in the Brazilian Institute of Geography and Statistics (IBGE) Ecological Reserve, an undisturbed Cerrado area within the Gama Cabeça-de-Veado Environmental Protection Area (APAGCV) [38], known before 1978 as Roncador Ecological Reserve (RECOR-as this experimental tower will be called henceforth).

Flux Tower Sites for Validation
For quality control of the validation data in the three experimental sites, the energy balance closure was considered using the ratio of turbulent fluxes (H + LE) to the available energy (Rn − G) at 30 min intervals during the daytime (8:00-18:00 h) [23].When the 30 min ratios were outside the range between 0.8 and 1.2, the flux data was discarded and replaced by a simple linear interpolation.If the database missed more than 5 of the 30 min intervals during the daytime measurements, the entire day was excluded from the analysis [23,46,47].LE and H values from the eddy covariance (EC) systems were recalculated using the Bowen ratio method (H/LE), as described by Reference [48].After the correction for H and LE, reference ET fraction (F _ec ) and crop coefficient (Kc _ec ) were calculated using the corrected values.Height and sensor characteristics for each experimental site are indicated in Table 1.USR site: The USR flux tower was installed on a sugarcane experimental field situated in Luiz Antonio municipality-São Paulo State (Figure 1a), which belongs to the Companhia Energética Santa Rita, at coordinates 21 • 38 13 S, 47 • 47 25 W, and at an altitude of 552 m.The area was larger than 400 hectares (ha), with a 1.4 m distance between planting rows and a slope of 1.87 ± 0.36 degrees (within the 1 km MODIS pixel centered in the flux tower) (Table 2), surrounded by pasture, citrus, and native savanna forest [43].The soil is characterized as Typic Haplustox, and its texture fractions are 22% clay, 74% sand, and 3% silt.Mean dry bulk density down to 2.6 m depth is 1500 kg m −3 .A denser layer was detected between 10 and 25 cm (d b = 1636 kg m −3 ), due to previous mechanical harvesting.The total ET in the June 2005 cycle was 829 mm (representing 69% of the precipitation), while in the July 2006 cycle it was 690 mm (51% of the precipitation) [43].Sugarcane was planted in 2003, and in 2004 and 2005 there were two harvests with straw burning.The data obtained from this site correspond to the first cycle of regrowth, beginning on 14 April 2005, extending until harvest on 11 May 2006 (totaling 393 days), and the second cycle, with harvest on 20 May 2007 (totaling 375 days).The fetch around the tower consisted of homogeneous sugarcane plantation in all directions within the 1 km MODIS pixel, with a mean canopy height ranging from plant emergence to ~3.5 m (harvest season) [49].Mean NDVI values (calculated from cloud free TM and ETM+ Landsat data for the study period) varied from 0.18 to 0.76, presenting a standard deviation below 0.08 (except for the sugarcane harvest season in 2006-see Figure S4 in the Supplementary Material).The instrumentation and variables measured by the EC system are described by [43,49].PDG site: This micrometeorological tower was installed at the reserve Pé-de-Gigante (PDG), an area of 1213 ha of undisturbed woodland, with an altitude varying between 590 and 740 m), in the Vassununga State Park (Instituto Florestal), in Santa Rita do Passa-Quatro municipality, São Paulo State [50].The tower coordinates are 21 • 37 9.26 S, 47 • 37 56.38 W, at an average elevation of 710 m.The mean slope within the 1 km MODIS pixel centered in the flux tower is 4.84 ± 2.62 degrees (Table 2), with a maximum value in the southeast direction of the flux tower, however, the wind frequency from this sector accounts for less than 15% of the daylight hours during the study period.The arboreal vegetation shows heights from 5 to 8 m, covering 50-70% of the park [49,51].The soil is a Quartzarenic Neosol (by the Brazilian Soil Classification System), with 87% sand, 3% silt and 10% clay [52].The micrometeorological tower provided measurements of Rn, soil heat flux (G), sensible (H) and latent (LE) heat fluxes, as well as auxiliary meteorological data measured at a height of 21 m and recorded every 30 min.The fetch of Cerrado vegetation around the flux tower was greater than 500 m in all directions [52], with mean NDVI (calculated from cloud free TM and ETM+ Landsat data for the study period) of 0.71 ± 0.03 during the study period, and a mean canopy height of 10 m [49].A more detailed description of the site, instrumentation, and data processing is provided by Reference [52].
RECOR site: This flux tower was located at IBGE Ecological Reserve, 26 km away from Brasília, at the coordinates 15 • 55 51.2 S, 47 • 52 21.3 W. and at an average elevation of 1140 m.The area is characterized mainly by cerrado grassland, shrub-land and woodland physiognomies, including patches of cerrado denso and cerradão [38].The tower is mounted over a sensu-stricto savanna vegetation, formed by typical species with an average height between 4-5 m [53].Within the 1 km MODIS pixel where the tower is located, the average slope was 2.69 ± 0.78 (Table 2), and the mean NDVI (calculated from Landsat data from 2011 to 2015) was 0.60 ± 0.04, showing that the Cerrado vegetation around the tower is homogeneous.This configuration represents 65% of the reserve area.Soils are predominantly Red-Yellow Latosols, Red-Dark latosols and Petroplinthic Soils.The average rainfall is 1453 mm yr −1 , with more than 75% of the volume concentrated from November to March.The average annual temperature is close to 22 • C, with monthly maximum and minimum values of 27 • C and 15.4 • C, respectively [54].The EC system was mounted between the height of 7 to 10 m (except the soil sensors) and provided energy balance fluxes and meteorological data every ten minutes for the period between 2011 and 2015.Further information on the experimental site and equipment is presented in the Large-Scale Biosphere-Atmosphere Experiment in Amazonia (LBA) project website [53].

Broadband Surface Albedo Submodel Adjustment
As pointed out by Reference [23], variations on surface albedo, as well as on other biophysical parameters, can indirectly affect the accuracy of ET estimations.Therefore, we adjusted the Liang equation [55] for the calculation of broadband surface albedo in the Brazilian Cerrado, a method used by several Land Surface Energy Balance (LSEB) models, in which albedo is calculated by integrating band reflectances across the shortwave spectrum [35,36], as shown in Equation ( 1) where α is the broadband surface albedo, ρ is surface reflectance for a band b, w b is the weighting coefficient representing the fraction of at-surface solar radiation occurring within the spectral for band b [35], calculated as: where R sλ is at-surface spectral hemisphere solar radiation for wavelength λ (µm), UP b and LO b are upper and lower-wavelength bounds assigned to MODIS band b (Table 3), which, in this application, includes all regions between the specific satellite bands (see Figure S3 in the Supplementary Material).According to Reference [36], the inclusion of wavelength regions between sensor bands in the determination of wb provides a more theoretically correct estimate of total surface reflectance.The integration limits were chosen from 300 to 4000 nm, which covers 98% of total at-surface solar radiation [36].For the adjustment of the coefficients wb, solar irradiance spectra were simulated for five combinations of land elevation, precipitable water content, and solar zenith angles (Table 4), following the methodology applied by [35,36], using the Simple Model of the Atmospheric Radiative Transfer of Sunshine (SMARTS2) [34].As observed by [35,36], regardless of the wide range of solar radiation settings, wb did not have significant changes.So, the adjusted model for estimating broadband surface albedo (for MODIS data) applied in this is study is presented in Equation (3). (3)

Satellite Data
To run METRIC, we used two MODIS/Terra products (collection 6, tiles h13v10 and h13v11): MOD09GA-500 m daily surface reflectance [56] and MOD11A1-1 km daily land surface temperature [57].These products are distributed by Land Processes Distributed Active Archive Center (LPDAAC) in HDF format (Hierarchy Data Format), with sinusoidal projection, and converted to GEOTIFF format through MRT (MODIS Reprojection Tool), developed by Earth Resources Observation System Data Center (EDC) and LPDAAC.The output for all products was defined for Geographic Coordinates (Lat-Long) and for Datum WGS84.
A digital elevation model (DEM) was obtained from SRTM (Shuttle Radar Topography Mission) with spatial resolution of 90 m (resampled to MODIS spatial resolution using bilinear interpolation) and utilized to account for changes in surface temperature caused by elevation differences [58,59].

METRIC Model
The METRIC model was implemented in Python language (Arcpy library from ArcGIS 10.3), using MODIS products and surface meteorological data as inputs.The processing dates (Table S1) were selected according to the availability of cloud-free images, low sensor view angle (preferably below 20 • ) in the AOI and soil moisture conditions.Moreover, the energy balance fluxes measured by the EC system were evaluated avoiding days without data.For São Paulo AOI (USR and PDG sites), we selected dates covering the sugarcane phenological stages, for cycles June 2005 and July 2006; and for the Federal district AOI we selected images throughout each year to represent the ET of natural vegetation and agriculture during dry and wet seasons.We highlight that for 2007 and 2014 the smaller number of simulations with METRIC is due to the lack of EC data for validation in USR/PDG and RECOR sites, respectively.
Energy balance is controlled by Rn, which is partitioned in three components (Equation ( 4)): G, H and LE.
where Rn, G, H and LE are all in W m −2 .Rn is the difference between incoming and outgoing radiation fluxes (shortwave and longwave) on the Earth's surface, [32,60] and is calculated by Equation ( 5) where K ↓ is downward shortwave radiation (W m −2 ), α and εs are, respectively, surface albedo (Equation ( 3)) and emissivity (obtained from MOD11 product), L w↓ is the thermal infrared radiation emitted by the atmosphere towards the surface (W m −2 ), and L w↑ is the thermal infrared radiation emitted by the surface (W m −2 ).K ↓ is calculated according to Equation ( 6) where G sc is the solar constant (1367 W m −2 ), θ z corresponds to the solar zenith angle (rad), d 2 is the square of the relative Earth-Sun distance (unitless), and τ is the atmospheric transmittance at the time of satellite overpass.Solar zenith angle is calculated from the solar angle β, obtained from MODIS image metadata.
Relative Earth-Sun distance is calculated as function of the day of the year (DOY), according to Reference [32]: Atmospheric transmissivity is computed for each pixel as a function of altitude (z) [32], applying the SRTM digital elevation model: The incoming longwave radiation emitted by the atmosphere is calculated using the Stefan-Boltzmann law: where ε a is the atmospheric emissivity, σ is the Stefan-Boltzmann constant (567.10 −8 W m −2 K −4 ) and T a is the near surface air temperature (K) [60].εa is estimated as described in [32]: The outgoing longwave radiation emitted from the surface was estimated using the Stefan-Boltzmann law (Equation ( 12)), applying the surface emissivity (εs) and the adjusted land surface temperature (T s_dem ) obtained from MOD11.
Q1 To estimate G with remote sensing data, an empirical approximation is conducted using Rn, surface albedo, Normalized Difference Vegetation Index (NDVI) [61], and T s_dem , as described in Equation ( 13): Sensible Heat Flux (H) corresponds to the heat exchanged between surface and air by conduction (Equation ( 14)) and is determined by the surface properties and the state of the atmospheric boundary layer.
where ρ air is air density (kg m −3 ), C p is air specific heat at constant pressure (J kg −1 K −1 ), dT represents the near-surface temperature difference between two near surface heights z 1 and z 2 and r ah is the aerodynamic resistance (s m −1 ), computed for each pixel as a function of aerodynamic roughness [32].
In the METRIC procedure, r ah and dT values are internally calibrated using a CIMEC process (Calibration using Inverse Modeling at Extreme Conditions), that identifies extreme or near extreme conditions in the image to use as endpoints (anchor points or cold/hot pixels), where ET can be independently estimated and assigned.The CIMEC process is required to produce an image-specific calibration of the energy balance processes to overcome biases associated with the radiometric accuracy of the satellite data, uncertainties in aerodynamic features of the landscape and associated with wind speed, sun angle, background thermal conditions and assumptions made in the model [62].Agricultural areas are usually selected to represent the hydrological extremes of the image, due to the advantage of presenting a better spatial uniformity of the vegetation, making it easy to identify them in orbital images.In addition, the algorithms for r ah , G, and emissivity usually present better results for this type of land cover [62].The intern calibration establishes a linear relationship between Ts and dT, in which ETr is applied to represent the ET of a crop at full vegetative development and with full water availability, where the surface tends to be cool due to evaporative cooling (cold pixel).An adjustment factor of 1.05 is typically used to reflect the actual evaporation amount at this anchor point.For a hot pixel, a non-vegetated agricultural area is selected, and a bare soil evaporation model, such as that presented by Reference [22], can be used to establish the adjustment factor which accounts for the residual soil evaporation from antecedent rainfall [59].
LE is obtained as the residual term of Equation ( 4).The daily ET a is calculated for each pixel from Equation (15), considering that the reference ET fraction (F) is constant throughout daytime (tested in Section 3.4), as described by Reference [46].
where ET a is the daily ET for each pixel (mm day −1 ), F i_M is the reference evapotranspiration fraction at the satellite overpass time (Equation ( 16)) and ET r is the daily accumulated reference evapotranspiration (mm day −1 ), calculated by the Penman-Monteith method [25].
where ET i_M is the instantaneous ET calculated by METRIC for each pixel (Equation ( 17)) and ETr h_i is the hourly reference ET at the satellite overpass time [32,46].
where ρ is water density, λ is the latent heat of vaporization, and 3600 converts seconds to hours.

Validation
The estimated energy balance fluxes and ET were compared with EC measurements in each experimental site, evaluating the model performance by regression analysis.The estimated values were extracted from the pixel (1 km × 1 km) centered on each flux tower (Figure 1), as performed by [60].For energy balance estimates, the comparison was made at the time of satellite overpass.For ET a, the estimates were compared to daily accumulated values.The source area (footprint) for turbulent fluxes was calculated using the 1D analytical model proposed by [63].Considering all METRIC simulation dates, the cumulative normalized flux (CNF) measurements of 90% were 780, 1120 and 723 m for USR, PDG and RECOR sites, respectively (Figure S6).The average values of the height-to-fetch ratio under near-stable conditions were 1:120, 1:59 and 1:90 for USR, PDG and RECOR sites, respectively.These values are close to those indicated by Reference [1] for the sugarcane experimental site and by Reference [52] for PDG site.
To assess the model's accuracy, we used the ratio of estimated and observed values (b), coefficient of determination (R 2 ), mean absolute error (MAE) (Equation ( 18)), and root mean square error (RMSE) (Equation ( 19)).RMSE is the overall error in the predictions relative to the actual measured value, while MAE measures the average magnitude of the errors in a set of predictions (without considering their direction), and R 2 determines the strength of the linear relationship between estimates and measurements [60,64].Furthermore, we used the Student's t-test (with a 95% confidence level) to evaluate the significance of the linear relationship between observed and estimated values [1,65].
where X obs is the flux tower observation, X est is the value estimated by METRIC, and n is the number of samples.

Energy Balance Closure for Flux Measurements
According to References [48,66,67], measures of H and LE tend to be underestimated by EC systems due to sources of errors such as non-homogeneous surface coverage, soil characteristics, instrumental errors, topography, divergence or flux dispersion, among others.In this sense, for METRIC processing dates (Table S1) the accuracy of EC data was checked using the energy balance closure for USR, PDG, and RECOR experimental sites for daytime flux measurements (8:00-18:00), as shown in Figure 2.  The hourly ratio of (H + LE) to (Rn − G) was unbalanced by 7%, 11%, and 4% for USR, PDG and RECOR sites, respectively, and the t-test indicated that the slope of the regression line through the origin was statistically different from unity, with a 95% confidence level.According to References [48,65] these are small uncertainties and can be corrected using the Bowen-ratio approach.In this manner, we applied this method to correct H and LE observed values, which were later used to validate METRIC estimations.

Broadband Surface Albedo
At the USR site, Equation (3) underestimated surface albedo by 10% (R 2 = 0.94, MAE = 0.02 and RMSE = 0.03), with a 95% confidence level (Table 5), presented in Section 3.3.The temporal profile presented in Figure 3a shows that the observed values vary during the sugarcane cycle, with higher values at the initial stage, when there is a greater contribution of soil in surface reflectance [49].During the crop development (indicated by an increase in NDVI-Figure 3a  The hourly ratio of (H + LE) to (Rn − G) was unbalanced by 7%, 11%, and 4% for USR, PDG and RECOR sites, respectively, and the t-test indicated that the slope of the regression line through the origin was statistically different from unity, with a 95% confidence level.According to References [48,65] these are small uncertainties and can be corrected using the Bowen-ratio approach.In this manner, we applied this method to correct H and LE observed values, which were later used to validate METRIC estimations.

Broadband Surface Albedo
At the USR site, Equation (3) underestimated surface albedo by 10% (R 2 = 0.94, MAE = 0.02 and RMSE = 0.03), with a 95% confidence level (Table 5), presented in Section 3.3.The temporal profile presented in Figure 3a shows that the observed values vary during the sugarcane cycle, with higher values at the initial stage, when there is a greater contribution of soil in surface reflectance [49].During the crop development (indicated by an increase in NDVI-Figure 3a), observed surface albedo decreases from 0.18-0.20 (at initial stage) to values around 0.15, due to the thickening of sugarcane canopy, which favors the extinction/absorption of solar radiation.At the maturation stage (after the wet to dry season transition) albedo increases again (reaching values of 0.20) due to a reduction in leaf area index and to the presence of senescent leaves, which alters the crop reflectance patterns [68].The reduction in observed and estimated albedo on 5 May 2005 is an effect of soil moisture reducing reflectance, caused by isolated precipitation in the experimental area in days prior to simulation [49].For June 2005 and July 2006 cycles the estimated surface albedo for sugarcane ranged from 0.13 to 0.18, similar to that described by Reference [69] for sugarcane in Red-Yellow Latosol and mechanical harvesting with straw burning.The adjusted equation showed a better performance for vegetative growth stage (full crop coverage), but at initial and maturation stages we observed a tendency to underestimate surface albedo, possibly due to a greater contribution of soil in the reflectance.
In addition, at the beginning of the two evaluated cycles, burning of harvest straw, soil plowing, and fertilization were carried out (procedures that decrease the reflectance observed by satellite).In the maturation stage, underestimates also occur, possibly due to residual soil moisture which decreases reflectance.
At the PDG site, surface albedo was underestimated by 15% (R 2 = 0.94, MAE = 0.02, and RMSE = 0.02), with a 95% confidence level.The observed albedo presents a seasonal variation, with an average value (α = 0.13) smaller than that observed at the USR site (α = 0.17).The lowest values occur during the dry season (α ∼ = 0.12), increasing gradually up to 0.16 in the dry to wet season transition.In the dry season, Ki penetrates more easily into the canopy, and thus the soil covered by dark plant litter absorbs more radiation, reducing the albedo [49].In the rainy season, on the other hand, the albedo is controlled by tree canopy reflectance, which creates shade on the ground and dominates the absorption of photosynthetically active and near-infrared radiation.During the year the adjusted albedo equation tends to underestimate this variable, especially during the dry season, which could be an effect of the interaction of Ki with tree canopy and substrate, which presents a dynamic of alternation in albedo control between wet and dry seasons.A better performance is observed for the wet season, when there is a greater vegetative vigor (indicated by higher NDVI values) and the woody vegetation canopies dominate the reflectance.be an effect of the interaction of Ki with tree canopy and substrate, which presents a dynamic of alternation in albedo control between wet and dry seasons.A better performance is observed for the wet season, when there is a greater vegetative vigor (indicated by higher NDVI values) and the woody vegetation canopies dominate the reflectance.S1).
For the sensu-strictu savanna at RECOR site, the adjusted model overestimated the surface albedo by 21%, with a 95% confidence level (R 2 = 0.96, MAE = 0.02 and RMSE = 0.03).The observed albedo shows lower values (0.09 to 0.12) than those obtained for PDG woody savanna (0.12 to 0.16) and does not show the same seasonal variability.The estimated albedo values, on the other hand, follow the NDVI series for the study period, varying between 0.09 and 0.16, with an increase in the wet season caused by a higher reflectance due to a greater amount of green vegetation, both in woody and grassy areas.S1).
For the sensu-strictu savanna at RECOR site, the adjusted model overestimated the surface albedo by 21%, with a 95% confidence level (R 2 = 0.96, MAE = 0.02 and RMSE = 0.03).The observed albedo shows lower values (0.09 to 0.12) than those obtained for PDG woody savanna (0.12 to 0.16) and does not show the same seasonal variability.The estimated albedo values, on the other hand, follow the NDVI series for the study period, varying between 0.09 and 0.16, with an increase in the wet season caused by a higher reflectance due to a greater amount of green vegetation, both in woody and grassy areas.

Validation of Energy Balance Fluxes
Table 5 presents a comparison between observed and estimated values of α, Rn, G, H, and LE, at the time of satellite overpass.For USR site (Figure 4a), the t-test indicated that Rn was overestimated by 9% (R 2 = 0.93, MAE = 59.84W m −2 , RMSE = 75.36W m −2 ), which is mainly due to underestimation of surface albedo at the crop initial stage.In vegetative growth and maturation stages the estimated Rn values were closer to 1:1 line, due to a better adjustment of the albedo submodel for conditions of full vegetation cover.It was found that the land surface temperature of MOD11 was underestimated by up to 10% at the dates when Rn was overestimated, which reduces the long wave radiation flux emitted by the surface and consequently increases Rn.Considering the study period, Rn ranged from 408 W m −2 to 728 W m −2 , close to those presented by Reference [70] for the same crop.As observed by Reference [49], Rn presents a seasonal pattern, mainly controlled by Ki, surface albedo, and land surface temperature.Regarding the longwave component, L w↓ ranged from 310 to 365 W m −2 , following the seasonal variations of air temperature.L w↑ , in a similar way, varied according to the surface temperature, ranging from 403 to 489 W m −2 (see Figure S7 on the Supplementary Material).G estimates for USR site did not show good agreement with the observed values, being underestimated by 31% (R 2 = 0.76), demonstrating that this submodel does not satisfactorily represent the reality for this landscape context.Sensible heat flux was overestimated by 26% (R 2 = 0.86), showing low correlation with the observed values.According to Reference [70], values of H for sugarcane in the wet season are controlled by Rn, while in dry season the maturation or harvest of the crop favors the H contribution in the energy partitioning, due to a lower leaf area index.During the wet season, the relative errors for H varied between −1% and 18%, while Reference [70] obtained relative errors up to 70% using SEBAL/METRIC.For LE, the estimated/observed values present a balanced distribution around the 1:1 line (Figure 4a), with R 2 = 0.94, MAE = 16.93W m −2 and RMSE = 21.43W m −2 (Table 5).The t-test indicated that the ratio of the estimated / observed values was significantly different from unity, suggesting that METRIC underestimated LE by 5%.This result shows that the errors in previous estimates were caused by the internal calibration in H estimates [46,71].Instantaneous LE estimates ranged from 197 to 569 W m −2 during the two sugarcane growth cycles, near the values obtained by Reference [70] for the same sugarcane area, that ranged from 146 W m −2 to 529 W m −2 .
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 26 At the PDG site (Figure 4b) Rn was overestimated by 9% (R 2 = 0.84, MAE = 58.83W m −2 , RMSE = 67.95W m −2 ), due to albedo underestimates and deviations in Ts, as verified for the USR site.METRIC estimates for this variable ranged from 454 to 772 W m −2 , with values generally higher than those obtained for the USR site at each date.This is a result of the albedo lower values in the woody savanna (Figure 3), as well as lower values of Ts (which reduce the longwave radiation emitted by the surface).Lw↓ and Lw↑ ranged from 290 to 357 W m −2 and from 394 and 477 W m −2 , respectively (Figure S7b), following the seasonal variation of air and surface temperatures.According to the t-test, G was overestimated by 39% (R 2 = 0.54) and H was underestimated by 42% (R 2 = 0.75).However, for LE the t-test indicated that b was significantly equal to unity, indicating that LE estimated by METRIC and in-situ measurements were similar (R 2 = 0.92 MAE = 48 W m −2 and RMSE = 56 W m −2 ).METRIC estimated LE values from 107 to 520 W m −2 at this site, following the seasonal variation of available energy and NDVI (Figure 3b).
The ratio between estimated and observed values of Rn at RECOR site (Figure 4c) was significantly equal to unity, with a 95% confidence level (R 2 = 0.96, MAE = 37 W m −2 , RMSE = 47 W m −2 ), showing that METRIC is effective in estimating this variable with MODIS data.Instantaneous At the PDG site (Figure 4b) Rn was overestimated by 9% (R 2 = 0.84, MAE = 58.83W m −2 , RMSE = 67.95W m −2 ), due to albedo underestimates and deviations in Ts, as verified for the USR site.METRIC estimates for this variable ranged from 454 to 772 W m −2 , with values generally higher than those obtained for the USR site at each date.This is a result of the albedo lower values in the woody savanna (Figure 3), as well as lower values of Ts (which reduce the longwave radiation emitted by the surface).L w↓ and L w↑ ranged from 290 to 357 W m −2 and from 394 and 477 W m −2 , respectively (Figure S7b), following the seasonal variation of air and surface temperatures.According to the t-test, G was overestimated by 39% (R 2 = 0.54) and H was underestimated by 42% (R 2 = 0.75).However, for LE the t-test indicated that b was significantly equal to unity, indicating that LE estimated by METRIC and in-situ measurements were similar (R 2 = 0.92 MAE = 48 W m −2 and RMSE = 56 W m −2 ).METRIC estimated LE values from 107 to 520 W m −2 at this site, following the seasonal variation of available energy and NDVI (Figure 3b).
The ratio between estimated and observed values of Rn at RECOR site (Figure 4c) was significantly equal to unity, with a 95% confidence level (R 2 = 0.96, MAE = 37 W m −2 , RMSE = 47 W m −2 ), showing that METRIC is effective in estimating this variable with MODIS data.Instantaneous Rn estimates ranged from 504 to 779 W m −2 , with an average 7% higher than that obtained at the PDG site.L w↓ varied between 308 and 370 W m −2 , while L w↑ ranged from 387 to 465 W m −2 (Figure S7c), values close to those obtained for woody savanna at PDG site.For G estimations, the model showed poor performance by overestimating this flux by 210% (R 2 = 0.67).Sensible heat flux was underestimated by 30% with R 2 = 0.92) W m −2 , absorbing previous deviations in its internal calibration.This led to accurate estimates of LE, whose values were considered equal to in-situ measurements by the t-test (b = 1), with a 95% confidence level.

Reference Evapotranspiration Fraction (F) and Actual Evapotranspiration (ETa)
The METRIC model assumes that the reference evapotranspiration fraction (F) is constant throughout the day for well-watered agricultural crops [65], and can be considered equal to its crop coefficient (Kc).Thus, we tested this assumption using the statistical analysis presented in Table 6.The test consists of (i) comparing the instantaneous (at satellite overpass time) F obtained from the eddy covariance system (F i_ec ) with the mean values of daytime hourly F (F mean ) and; (ii) comparing F i_ec with the Kc obtained from EC system (Kc _ec ), as performed by Reference [46].Table 6 shows that the assumption was considered valid for the three experimental sites, which indicates that F was approximately constant during the daytime under clear sky conditions (with a confidence level of 95%), for different phenological stages of both sugarcane and natural vegetation.Thus, we applied Equation (15) to extrapolate our ETa estimates from instantaneous to daily values.Table 6.Statistical analysis results for reference evapotranspiration fraction (F) and daily evapotranspiration (ETa).F i_ec is the instantaneous observed F, F mean is the diurnal (8:00-18:00) mean value of F, Kc_ ec is the Kc observed for each experimental site, Kc_ metric is the Kc estimated by METRIC, ETa obs is the observed daily evapotranspiration, and ETa metric is the daily evapotranspiration simulated by METRIC.Table 6 indicates that the crop coefficients estimated by METRIC for sugarcane (USR site) were significantly equivalent to those measured by the EC system (Kc_ metric = Kc _ec ), according to the t-test, with a 95% of confidence level.The temporal profile of sugarcane Kc in the June 2005 and July 2006 cycles is shown in Figure 5, and grouped according to the crop phenological phases, as described in References [22,72].The deviations obtained in Kc _metric at some dates may have been caused by spectral mixing of different land cover types within the same MODIS pixel (1 km for Kc), with different reflectance values and surface temperature [73].The mean estimated Kc values per phenological stage were 0.16 for germination, 0.57 for tillering, 0.57 for vegetative growth and 0.50 for maturation stages, while Reference [22] indicates values of 0.40, 1.25, and 0.75 for the initial, development and late stages, respectively.Reference [74] estimates the sugarcane Kc at the initial stage (up to 50 days after planting) between 0.35 and 0.40, which was obtained in this study for the June 2005 cycle.During this period, the rainfall was above average, which increased the water availability allowing a more accentuated growth than that for the July 2006 cycle.In the woody savanna of PDG site, the t-test also indicated that the estimated F(Kc_metric) as well as the ETa estimates were equivalent to the measured values (b = 1), with a 95% confidence level.Figure 7b shows that the estimated ETa varied between 2.48 and 6.56 mm day −1 with the same seasonal pattern verified on the USR site, despite having higher values.For this area, Reference [52] observed ETa values varying from 1.0 mm day −1 (dry season of 2010) to 7.1 mm day −1 (summer season of 2012), and Reference [4] obtained an average value of 2.81 mm day −1 for the wet season and 1.70 mm day −1   ETa estimates for the USR site (Figure 6a) are close to the 1:1 line, with some estimate deviations during the rainy season.It is important to point out that for these dates it was difficult to select suitable hot pixels for the internal calibration, because the agricultural areas were almost entirely covered with the sugarcane crop.But in general METRIC performed well in estimating ETa (R 2 = 0.94, MAE = 0.21 mm day −1 and RMSE = 0.35 mm day −1 ), as indicated by the t-test (b = 1, with a 95% confidence level).The temporal profile of ETa for sugarcane (Figure 7a) shows that during the two cycles ETa varied according to crop development (indicated by Kc in Figure 5 and NDVI in Figure 3a), the rainfall distribution, and energy availability (data not shown).Estimated ETa ranged from 0.29 mm day −1 to 5.21 mm day −1 , with lower values occurring after harvest (when ETa is dominated by soil evaporation), gradually increasing to the highest values at peak development (higher water availability and leaf area index), followed by a decrease in the maturation stage (senescence of leaves reduces transpiration and increases H) [75].

Site
In the woody savanna of PDG site, the t-test also indicated that the estimated F(Kc_ metric ) as well as the ETa estimates were equivalent to the measured values (b = 1), with a 95% confidence level.Figure 7b shows that the estimated ETa varied between 2.48 and 6.56 mm day −1 with the same seasonal pattern verified on the USR site, despite having higher values.For this area, Reference [52] observed ETa values varying from 1.0 mm day −1 (dry season of 2010) to 7.1 mm day −1 (summer season of 2012), and Reference [4] obtained an average value of 2.81 mm day −1 for the wet season and 1.70 mm day −1 for the dry season between 2001 and 2003.In the woody savanna of PDG site, the t-test also indicated that the estimated F(Kc_metric) as well as the ETa estimates were equivalent to the measured values (b = 1), with a 95% confidence level.Figure 7b shows that the estimated ETa varied between 2.48 and 6.56 mm day −1 with the same seasonal pattern verified on the USR site, despite having higher values.For this area, Reference [52] observed ETa values varying from 1.0 mm day −1 (dry season of 2010) to 7.1 mm day −1 (summer season of 2012), and Reference [4] obtained an average value of 2.81 mm day −1 for the wet season and 1.70 mm day −1    For the sensu-stricto savanna of the RECOR site, unlike the other experimental sites, the t-test indicated that b values for Kc and ETa validation were significantly different from unity, then METRIC overestimated Kc by 11% (R 2 = 0.96, MAE = 0.07 and RMSE = 0.08) and ETa by 14% (R 2 = 0.96, MAE = 0.62 mm day −1 and RMSE = 0.75 mm day −1 ), with a 95% confidence level.Figure 6c shows that the ETa estimates were always above the 1:1 ratio line, which indicates that METRIC tended to overestimate daily ET during both wet and dry seasons over the sensu-stricto savanna.The ETa profile presented in Figure 7c shows that the estimated ETa varied between 2.71 and 6.91 mm day −1 , following the same seasonal pattern obtained in PDG site: higher values during the wet season and lower values during the dry season.For the sensu-stricto savanna of the RECOR site, unlike the other experimental sites, the t-test indicated that b values for Kc and ETa validation were significantly different from unity, then METRIC overestimated Kc by 11% (R 2 = 0.96, MAE = 0.07 and RMSE = 0.08) and ETa by 14% (R 2 = 0.96, MAE = 0.62 mm day −1 and RMSE = 0.75 mm day −1 ), with a 95% confidence level.Figure 6c shows that the ETa estimates were always above the 1:1 ratio line, which indicates that METRIC tended to overestimate daily ET during both wet and dry seasons over the sensu-stricto savanna.The ETa profile presented in Figure 7c shows that the estimated ETa varied between 2.71 and 6.91 mm day −1 , following the same seasonal pattern obtained in PDG site: higher values during the wet season and lower values during the dry season.vegetative vigor (bright red in Figure 8d), presenting higher NDVI (mean of 0.78) and maximum values of ETa reaching 4.5 mm day −1 .The area around PDG site shows the highest values of NDVI and ETa both in dry and wet seasons, since it is a region surrounded by reforestation and Cerrado woodlands.

Spatial Distribution of ETa
At the Federal District AOI, the agricultural areas (see Figure S1) on 4 January 2014 (Figure 9a) are at full vegetative vigor (soybean crop in yellow and maize in red, mainly), with NDVI and ETa values above 0.9 and above 6.0 mm day −1 , respectively.In these areas there are rain fed and irrigated fields, causing variability in ETa values estimated by METRIC.In an irrigated field, soil water availability is higher, which associated to a high evaporative demand (characteristic of this region) leads to higher ETa values.On 13 July 2015 (Figure 9b) NDVI values are generally smaller throughout this AOI (even for natural vegetation areas), but high values of NDVI (0.9) and ETa (above 4 mm day −1 ) occur in areas of irrigated agriculture (predominantly maize crop).However, for a large part of the AOI, low ETa values occur (<1.0 mm day −1 ), evidence of a reduction in soil moisture during the dry season.

Discussion
The albedo submodel was adapted for the study areas and applied with MOD09 data as input, but the t-test indicated, with a 95% confidence level, that albedo was underestimated at the USR and PDG sites by 10% and 15%, respectively, and overestimated by 21% at RECOR site (Table 5).For the USR site, Equation (3) presented a better result when compared to those obtained by Reference [76] for this experimental site with the equations proposed by Reference [36,55], which underestimated the surface albedo by 19% and 24%, respectively.We must also consider that during the harvest season there may be spectral mixing of harvested areas and sugarcane crop within the MODIS pixel (Figure S5 on the Supplementary Material).Moreover, for taller vegetation (as woody savanna), deviations in estimated albedo may occur in situations where the sun is not at its zenith, which causes the mixing of sunlit and shaded areas in the same pixel, particularly in the winter months, when the solar zenith angle is higher [59].
The method applied to estimate surface albedo does not correct for any bidirectional reflectance distribution function (BRDF) caused by off-nadir sensor view angles, such as those occurring in MODIS images (sensor view angles up to 55 degrees), where backscattering or forward scattering directions causes deviations in the reflected signals to sensor [36].Then, for dates where the sensor view angle was close to 20 degrees off-nadir (maximum sensor view angle recommended by [77]), the accuracy of MOD09 surface reflectance may be affected due to lack of consideration of BRDF, leading to deviations in the estimated surface albedo.In addition, MOD09 atmospheric correction does not take into consideration local atmospheric parameters, which can lead to errors [78].
NDVI values for PDG site (woody savanna) ranged from 0.60 to 0.80, while for RECOR (sensustricto savanna) it varied between 0.50 and 0.75, but both presented the same seasonal pattern, At the Federal District AOI, the agricultural areas (see Figure S1) on 4 January 2014 (Figure 9a) are at full vegetative vigor (soybean crop in yellow and maize in red, mainly), with NDVI and ETa values above 0.9 and above 6.0 mm day −1 , respectively.In these areas there are rain fed and irrigated fields, causing variability in ETa values estimated by METRIC.In an irrigated field, soil water availability is higher, which associated to a high evaporative demand (characteristic of this region) leads to higher ETa values.On 13 July 2015 (Figure 9b) NDVI values are generally smaller throughout this AOI (even for natural vegetation areas), but high values of NDVI (0.9) and ETa (above 4 mm day −1 ) occur in areas of irrigated agriculture (predominantly maize crop).However, for a large part of the AOI, low ETa values occur (<1.0 mm day −1 ), evidence of a reduction in soil moisture during the dry season.

Discussion
The albedo submodel was adapted for the study areas and applied with MOD09 data as input, but the t-test indicated, with a 95% confidence level, that albedo was underestimated at the USR and PDG sites by 10% and 15%, respectively, and overestimated by 21% at RECOR site (Table 5).For the USR site, Equation (3) presented a better result when compared to those obtained by Reference [76] for this experimental site with the equations proposed by Reference [36,55], which underestimated the surface albedo by 19% and 24%, respectively.We must also consider that during the harvest season there may be spectral mixing of harvested areas and sugarcane crop within the MODIS pixel (Figure S5 on the Supplementary Material).Moreover, for taller vegetation (as woody savanna), deviations in estimated albedo may occur in situations where the sun is not at its zenith, which causes the mixing of sunlit and shaded areas in the same pixel, particularly in the winter months, when the solar zenith angle is higher [59].
The method applied to estimate surface albedo does not correct for any bidirectional reflectance distribution function (BRDF) caused by off-nadir sensor view angles, such as those occurring in MODIS images (sensor view angles up to 55 degrees), where backscattering or forward scattering directions causes deviations in the reflected signals to sensor [36].Then, for dates where the sensor view angle was close to 20 degrees off-nadir (maximum sensor view angle recommended by [77]), the accuracy of MOD09 surface reflectance may be affected due to lack of consideration of BRDF, leading to deviations in the estimated surface albedo.In addition, MOD09 atmospheric correction does not take into consideration local atmospheric parameters, which can lead to errors [78].
NDVI values for PDG site (woody savanna) ranged from 0.60 to 0.80, while for RECOR (sensu-stricto savanna) it varied between 0.50 and 0.75, but both presented the same seasonal pattern, following the rainfall distribution (Figure S2) [39].At the time when NDVI reaches its maximum value, the estimated albedo begins to decrease, due to the shade provided by taller tree canopies, which favors the extinction of solar radiation at the surface and reduces the reflected shortwave radiation.NDVI decreases during the dry season can be attributed to leaf abscission that occurs during the wet to dry season transition for many Cerrado tree species [79].The observed values of albedo in the RECOR site do not present the same seasonal variation observed for the woody savanna of the PDG site, possibly due to the spatial distribution of the different vegetal species in stricto-sensu savanna, where the woody and grassy vegetation mixture does not present the same alternation in albedo control observed for woody savanna.
Net radiation was satisfactorily estimated by METRIC at the stricto-sensu savanna but overestimated by 8% and 9% at sugarcane and woody savanna experimental sites, respectively, possibly due to deviations in the estimated surface albedo or to the land surface temperature obtained from MOD11.Our results showed a better fit than those obtained by Reference [1] (R 2 = 0.83 for sugarcane and R 2 = 0.88 for savanna), but with higher RMSE ( [1] obtained RMSE = 37 W m −2 for sugarcane and RMSE 47 W m −2 for savanna).Rn is typically higher during the wet season, which is a typical pattern for the Cerrado biome [1,79], although it is the opposite of that observed for the humid evergreen forest of the Amazon Basin, where prolonged cloud cover tends to reduce wet-season Rn [80].For Cerrado vegetation, the decline in Rn during the dry season is due not only to lower solar radiation, but also in part to an increase in the surface albedo (Figure 3) caused by a reduction in vegetation leaf area and/or greenness during the dry season, when soil water availability is low [79].The longwave components of Rn (L w↓ and L w↑ ) are estimated using the Stefan-Boltzmann law (Equations ( 10) and ( 12)), therefore, they are directly influenced by air and surface temperatures, respectively [60].Thus, for the three experimental sites, L w↓ and L w↑ curves have a slight increase during summer, when both air and surface temperatures are higher.Comparing these two components (Figure S7), L w↑ presented higher values than L w↓ , once surface temperature obtained from MOD11 was higher than air temperature (at satellite overpass time).Hence, the resulting longwave component of Rn is always negative, with a seasonal variation less pronounced than that observed for the shortwave component.Considering the three experimental sites evaluated, Rn ranged from 272.34 W m −2 (sugarcane) to 837.80 W m −2 (sensu-stricto savanna).As these estimates are obtained at satellite overpass time (~10:30), and since measurements by remote sensing are possible for clear sky days, our results for METRIC processing days represent the instances of greater energy availability to the biophysical processes [70].
Regarding G, METRIC estimates showed low agreement with field measurements for all experimental sites.Generally, this flux is the one that presents the biggest modeling errors with satellite data, since it is an empirical approximation specifically calibrated for semi-arid areas and may not accurately represent the local conditions of our study area [1,18,33,70,81], and these errors can be associated with the difference of spatial and temporal scales between satellite data and in situ measurements [23,65].LE estimates were significantly equal to observations for PDG and RECOR sites, but underestimated for the sugarcane experimental site, although only by 5%.The estimated values of LE ranged from 18.8 W m −2 to 376.53 W m −2 in the sugarcane field (USR site), 107 to 520 in PDG site and 197 to 569 W m −2 in RECOR site (sensu-stricto savanna), following the seasonal pattern indicated by Reference [1] for sugarcane and by Reference [79] for a mixed grassland savanna in Mato Grosso, Brazil.The accuracy in LE estimations depends on the anchor pixel selection, which determines the self-calibration of H [65]. Thus, the biases in H are necessary because its intern calibration absorbs errors derived from other steps in previous estimates (such as albedo, Rn, and G) [65,71,82].
Our results corroborate to the assumption that F is constant throughout the day (and is equivalent to Kc) at the three experimental sites, showing that the METRIC method for extrapolating from instantaneous to daily ET is valid for heterogeneous surfaces containing agricultural and natural Cerrado vegetation.According to Reference [23] Kc values for heterogeneous surfaces depend on canopy geometry, plant density, and on adopted agricultural practices, which determine the ET partition between transpiration and evaporation.For sugarcane it was verified that Kc varied in each phenological stage according to increases in leaf area index and water availability.Kc values below those referenced by Reference [22], especially for vegetative growth stage, may be due to the spatial resolution for Kc maps (1 km), which can mix different land cover types with sugarcane in the same pixel.
Overall, METRIC was effective at estimating ETa at USR and PDG sites, with a 95% confidence level.For the RECOR site, ETa was overestimated by 14%, with R 2 and RMSE close to those reported in the literature.Reference [83] reported values of R 2 ranging from 0.45 to 0.95 and RMSE varying between 10 to 30% for ETa estimates in natural ecosystems using remote sensing based models, and stated that the uncertainty associated with remote sensing estimates of ET is constrained by the accuracy of ground measurements, for which for the flux tower data are in the order of 10-30%.For fully covered crops, several studies [11,17,65,84] have indicated that METRIC could present absolute errors between 5-20% in the daily ETa estimate.Reference [59] obtained R 2 = 0.67 and RMSE = 0.81 mm for a Tropical semideciduous forest, in which METRIC tended to overestimate daily ET in comparison to EC flux tower data.They report that the error can be associated with the selection of cold pixels, since in their study area there are no ideal agricultural areas to be selected for cold pixels.
Sugarcane ETa presented a higher variability (ranging from 0.29 to 5.21 mm day −1 ) due to crop variation from bare soil in the beginning of the cycle, when ET is controlled mainly by soil evaporation, to the maximum vegetative vigor, when the NDVI and leaf area index reach their maximum values.However, the woody savanna ETa was always higher, even in the dry season, because of the greater density of the woody component (predominantly deep rooted), that can draw water stored in deeper soil compartments.According to Reference [39] deep soil compartment can contribute to as much as 83% of the total water used in a tree-dominated Cerrado landscape, also providing water for vegetation in other periods of high demand, such as in the transition from dry to rainy season.Reference [52] showed that, at the PDG site between 2009 and 2012, despite a 11% increase in the annual rainfall observed in an analysis year, water use increased by 4%, while a 16% decrease in rainfall (in the subsequent year) did not affect ET.They also observed that in woody savanna the vegetation exhibited a strong surface control in parallel with a decrease in leaf area index, thus reducing the demand for water on a gradually declining water availability, an evidence of its adaptation to this seasonally dry climate.
Seasonality in ETa estimates at the RECOR site was comparable to those obtained for the PDG site, with minimum values occurring at the end of the dry season, indicating that the reduced water availability constrained dry season ETa at both sites.Differences in maximum and minimum values may be a result of the savanna vegetation spatial heterogeneity, enhanced by the distinctive leaf phenology of the species [52,85,86].Comparing two phytophysiognomies (cerrado denso and campo sujo) during the dry and wet seasons, Reference [39] estimated ETa ranging from 1.4 to 5.8 mm day −1 for cerrado denso and from 0.9 to 4.5 mm day −1 for campo sujo.They stated that the differences in the two ecosystems are associated not only with differences in root distribution, but also with differences in tree densities and phenology of full-leaf canopies.In general, METRIC showed a good performance, considering that Reference [17] had suggested that remote sensing-based ETa estimates with similar biases are considered as acceptable.

Conclusions
We evaluated METRIC model performance to estimate the energy balance and ET in heterogeneous areas in the Brazilian cerrado.Flux tower data from sugarcane, woody savanna, and sensu-stricto savanna were used to validate METRIC estimates.The surface albedo model was adapted applying the SMARTS2 model, leading to deviations between 10-21% and maximum RMSE = 0.03.
Our results suggest that METRIC is effective for simulating the energy balance and ET spatial and temporal variability over agricultural and natural land covers in the Brazilian Cerrado with relative accuracy.Overall, METRIC performed well for Rn (R 2 > 0.94) and LE (R 2 > 0.93), leading to accurate ETa estimates for USR and PDG sites (R 2 = 0.94 and R 2 = 0.88, respectively) but showed an overestimation of 14% for the RECOR site.
Furthermore, the assumption that F is constant throughout the day and can be used to extrapolate instantaneous to daily ET values was tested and validated for our three experimental sites.For sugarcane we obtained Kc curves covering its phenological stages, which can be applied in yield models for this crop in future studies.
Sugarcane energy balance fluxes and ET varied throughout the crop cycles, showing dependence on available energy, rainfall, albedo, and NDVI.For the natural cerrado vegetation, ET was mainly controlled by the seasonal variation of net radiation, rainfall, and surface water availability.Our results reproduced a typical pattern for the Brazilian cerrado, which is marked by the seasonality of leaf area index and green biomass.
Many species of woody and stricto-sensu savannas have deep root systems, which allow transpiration to be maintained even during the dry season.Thus, the conversion of undisturbed cerrado to grass-dominated or agricultural crops (combined with a higher fire frequency) has the potential to change the energy balance by increasing surface albedo and temperature [10], as well as water balance, by increasing runoff and groundwater recharge, decreasing canopy interception and reducing ET due to less water uptake by deep roots (especially during the dry season).Therefore, assuming the local precipitation is reliant on regional rates of ET, the cerrado climate could become drier [4,39,52,87].Lastly, as the Brazilian Cerrado is a large biome with different conditions of land cover, soil types, vegetation, and hydrometeorology, it is necessary to carry out more studies in different regions in order to improve our understanding of ET dynamics.

Figure 1 .
Figure 1.Study Area in the Brazilian Cerrado-(a) São Paulo area of interest and (b) Federal District area of interest-Color composites from MOD09GA (left), OLI/Landsat-8 (5 January 2014 ) for Brasília study area (upper right) and TM/Landsat-5 (14 April 2016) for São Paulo study area (bottom right).The flux tower sites are indicated on the map.

Figure 1 .
Figure 1.Study Area in the Brazilian Cerrado-(a) São Paulo area of interest and (b) Federal District area of interest-Color composites from MOD09GA (left), OLI/Landsat-8 (5 January 2014 ) for Brasília study area (upper right) and TM/Landsat-5 (14 April 2016) for São Paulo study area (bottom right).The flux tower sites are indicated on the map.
), observed surface albedo decreases from 0.18-0.20 (at initial stage) to values around 0.15, due to the thickening of sugarcane canopy, which favors the extinction/absorption of solar radiation.At the maturation stage (after the wet to dry season transition) albedo increases again (reaching values of 0.20) due to a reduction in leaf area index and to the presence of senescent leaves, which alters the crop reflectance patterns [68].The reduction in observed and estimated albedo on 5 May 2005 is an effect of soil moisture reducing reflectance, caused by isolated precipitation in the experimental area in days prior

Figure 4 .
Figure 4. Comparison at the time of satellite overpass between estimated and observed energy balance values of net radiation (Rn), soil heat flux (G), sensible heat flux (H), and latent heat flux (LE) for (a) USR (sugarcane), (b) PDG (woody savanna), and (c) RECOR (sensu-stricto savanna) experimental sites.The dotted line represents the 1:1 relation.

Figure 4 .
Figure 4. Comparison at the time of satellite overpass between estimated and observed energy balance values of net radiation (Rn), soil heat flux (G), sensible heat flux (H), and latent heat flux (LE) for (a) USR (sugarcane), (b) PDG (woody savanna), and (c) RECOR (sensu-stricto savanna) experimental sites.The dotted line represents the 1:1 relation.

26 Figure 5 .
Figure 5. Daily values of crop coefficient estimated by METRIC (Kc_METRIC) and observed by the eddy covariance system (Kc_ec) in the June 2005 and July 2006 sugarcane cycles at USR site.Shaded areas separate the phenological stages of (a) germination, (b) tillering, (c) growth, and (d) maturation.DAP = days after planting.

Figure 5 .
Figure 5. Daily values of crop coefficient estimated by METRIC (Kc _METRIC ) and observed by the eddy covariance system (Kc _ec ) in the June 2005 and July 2006 sugarcane cycles at USR site.Shaded areas separate the phenological stages of (a) germination, (b) tillering, (c) growth, and (d) maturation.DAP = days after planting.

Figure 5 .
Figure 5. Daily values of crop coefficient estimated by METRIC (Kc_METRIC) and observed by the eddy covariance system (Kc_ec) in the June 2005 and July 2006 sugarcane cycles at USR site.Shaded areas separate the phenological stages of (a) germination, (b) tillering, (c) growth, and (d) maturation.DAP = days after planting.

Figure 7 .
Figure 7. Temporal profiles of ETa (mm day −1 ) observed by the eddy covariance system (blue) and estimated by METRIC (orange), during the study period in each experimental site: (a) USRsugarcane, (b) PDG-woody savanna, and (c) RECOR-sensu-stricto savanna.

Figure 7 .
Figure 7. Temporal profiles of ETa (mm day −1 ) observed by the eddy covariance system (blue) and estimated by METRIC (orange), during the study period in each experimental site: (a) USR-sugarcane, (b) PDG-woody savanna, and (c) RECOR-sensu-stricto savanna.

Figures 8
Figures 8 and 9 present the spatial distribution of ETa for São Paulo State and Federal District AOIs, illustrating the difference in the patterns observed during the dry and wet seasons.For the São Paulo AOI ETa varies between the dry and the wet seasons, reflecting the conditions of soil water availability, precipitation distribution (Figure S2a,b), and vegetative vigor (the increase of ETa when the sugarcane crop is at full development during the wet season).The 3 October 2005 image is from the dry to wet season transition, when the sugarcane crop is still in the harvest period, with several areas of bare soil (dark blue) or covered with straw left on the soil after harvesting (cyan), which are the same areas that appear in bright red (fully vegetated) in the 1 April 2006 image (Figure8d).In the dry season, NDVI and ETa values (Figure8b,c, respectively) are generally lower for agricultural areas than those in the wet season (Figure8e,f, respectively).On 1 April 2006, sugarcane is at maximum vegetative vigor (bright red in Figure8d), presenting higher NDVI (mean of 0.78) and maximum values of ETa reaching 4.5 mm day −1 .The area around PDG site shows the highest values of NDVI and ETa both in dry and wet seasons, since it is a region surrounded by reforestation and Cerrado woodlands.

Figure 8 .
Figure 8. RGB Color composites for MODIS bands (RGB261), Normalized Difference Vegetation Index (NDVI), and actual evapotranspiration (ETa) simulated by METRIC on (a-c) 3 October 2005 (transition from dry to wet season) and (d-f) 1 April 2006 (end of wet season) for the study in northern São Paulo State (near USR and PDG sites).Gray pixels are urban areas extracted from TerraClass thematic map.

Figure 8 .
Figure 8. RGB Color composites for MODIS bands (RGB261), Normalized Difference Vegetation Index (NDVI), and actual evapotranspiration (ETa) simulated by METRIC on (a-c) 3 October 2005 (transition from dry to wet season) and (d-f) 1 April 2006 (end of wet season) for the study in northern São Paulo State (near USR and PDG sites).Gray pixels are urban areas extracted from TerraClass thematic map.

Figure 9 .
Figure 9. RGB Color composites for MODIS bands (RGB261), Normalized Difference Vegetation Index (NDVI), and actual evapotranspiration (ETa) simulated by METRIC for dates of (a-c) 4 April 2014 (wet season) and (d-f) 13 July 2015 (dry season) for the study area near Brasília-DF.Gray pixels are urban areas extracted from TerraClass thematic map.

Figure 9 .
Figure 9. RGB Color composites for MODIS bands (RGB261), Normalized Difference Vegetation Index (NDVI), and actual evapotranspiration (ETa) simulated by METRIC for dates of (a-c) 4 April 2014 (wet season) and (d-f) 13 July 2015 (dry season) for the study area near Brasília-DF.Gray pixels are urban areas extracted from TerraClass thematic map.

Supplementary Materials:
The Supplementary Materials are available online at http://www.mdpi.com/2072-4292/10/6/1181/s1.Author Contributions: B.S.O., E.C.M. and M.C.-B.design the research.B.S.O., M.C.-B., G.B. and G.A.V.M. processed the data.B.S.O., E.C.M. and M.C.-B.contributed to the analysis and interpretation of the results.B.S.O. and E.C.M. wrote the paper and all the other members revised it.

Table 1 .
Sensor heights and configuration inside USR, PDG and RECOR experimental sites.

Table 2 .
Zonal statistics for slope (degrees) obtained for the 1 km MODIS pixel centered in each micrometeorological tower (USR, PDG and RECOR).

Table 3 .
MODIS band widths and ranges for UP b and LO b applied in Equation (2).

Table 4 .
Levels of solar Zenith angle, Precipitable water, and Local elevation applied as input in SMARTS2 to simulate at-surface solar radiation spectra.
for the dry season between 2001 and 2003.
for the dry season between 2001 and 2003.