Surface Albedo and Temperature Models for Surface Energy Balance Fluxes and Evapotranspiration Using SEBAL and Landsat 8 over Cerrado-Pantanal, Brazil

The determination of the surface energy balance fluxes (SEBFs) and evapotranspiration (ET) is fundamental in environmental studies involving the effects of land use change on the water requirement of crops. SEBFs and ET have been estimated by remote sensing techniques, but with the operation of new sensors, some variables need to be parameterized to improve their accuracy. Thus, the objective of this study is to evaluate the performance of algorithms used to calculate surface albedo and surface temperature on the estimation of SEBFs and ET in the Cerrado-Pantanal transition region of Mato Grosso, Brazil. Surface reflectance images of the Operational Land Imager (OLI) and brightness temperature (Tb) of the Thermal Infrared Sensor (TIRS) of the Landsat 8, and surface reflectance images of the MODIS MOD09A1 product from 2013 to 2016 were combined to estimate SEBF and ET by the surface energy balance algorithm for land (SEBAL), which were validated with measurements from two flux towers. The surface temperature (Ts) was recovered by different models from the Tb and by parameters calculated in the atmospheric correction parameter calculator (ATMCORR). A model of surface albedo (asup) with surface reflectance OLI Landsat 8 developed in this study performed better than the conventional model (acon) SEBFs and ET in the Cerrado-Pantanal transition region estimated with asup combined with Ts and Tb performed better than estimates with acon. Among all the evaluated combinations, SEBAL performed better when combining asup with the model developed in this study and the surface temperature recovered by the Barsi model (Tsbarsi). This demonstrates the importance of an asup model based on surface reflectance and atmospheric surface temperature correction in estimating SEBFs and ET by SEBAL.


Introduction
Surface energy balance fluxes (SEBFs) are one of the most important biophysical processes in environmental and hydrological studies [1][2][3]. SEBFs represent the processes of partitioning of available energy on the surface, measured by the net radiation (Rn), to evapotranspiration (ET) and soil and air heating, represented by soil heat flux (G) and parameters for T s recovery models was developed by NASA, using the atmospheric correction parameter calculator (ATMCORR). ATMCORR stands out for its simple operation, considering that its platform is online and requires only some meteorological data and for its robustness, and it can be applied for TM, ETM sensors, and TIRS in different latitudes over long periods [29,36]. This online platform integrates radioactive transfer codes (MODTRAN v4.0) with data from the National Centers for Environmental Prediction (NCEP) [29,36].
Given the importance of estimating SEBFs and ET from the a sup , which is in turn estimated by the surface reflectance and the T s without atmosphere and the emissivity corrections, the objective of this study is to evaluate the performance of the a sup and T s recovery models for the estimation of SEBFs and ET by SEBAL in the Cerrado-Pantanal transition region of the state of Mato Grosso, Brazil. This transition zone consists of upland cerrado vegetation that grades into an extensive wetland complex, with natural woodlands, forests, grasslands, and human land covers, such as agriculture, pasture, and urban areas, which affect Ts, albedo, and other local climatic variables that are important for SEBFs and ET [1].

Study Site
The study area is in the transition region Cerrado-Pantanal, covering path 226 and row 71 of the satellite Landsat 8 in southern Mato Grosso, Brazil ( Figure 1). Data from two flux towers were used, one in the Cerrado and the other in the Pantanal. The Cerrado tower is located in Fazenda Miranda (FMI) (15 • 43 55 S; 56 • 4 19 W), approximately 15 km south of the city of Cuiabá. The vegetation at the FMI is dominated by native and exotic grasses and by the semi-deciduous trees Curatella americana L. and Diospyros hispida A.DC [37], and the soil is classified as Plinthosols [38]. The Pantanal flux tower is in the Baia das Pedras (BPE) of the Estância Ecológica SESC-Pantanal (16 • 29 52 S; 56 • 24 44 W), in the municipality of Poconé, approximately 160 km from Cuiabá. The predominant vegetation in BPE is composed of the tree Combretum lanceolatum Pohl [39], and the soil is classified as Gleysols [38]. The BPE topography is flat, with flooding occurring from January to June. The Köppen-Geiger climate classification of the entire study region is Aw [40]. Annual rainfall is 1372 mm, with a dry season from May to September and a wet season from October to April, and average annual temperature of 26.9 • C [41].
Four types of land uses (agriculture, urban areas, forest, and water bodies) were sampled in the study area to develop a surface albedo model using surface reflectance from the OLI Landsat 8 ( Figure 1). The types of coverage were strategically selected because they represent an area of 9 pure pixels (3 × 3 pixel matrix) to minimize the influence of neighboring types of coverage. The agricultural areas are located northeast of the study area (yellow circles) and comprise a plateau area, with a predominance of soybean and corn planting. The urban areas are inserted in the urban perimeters of the municipalities of Cuiabá and Várzea Grande in densely urbanized regions (red circles). Forests comprise large forest fragments and permanent preservation areas close to rivers (green circles). The areas of water bodies are inserted in the extensive system of Chacororé and Sinhá Mariana bays (blue circles), with areas of up to 64.92 km 2 and 11.25 km 2 , respectively.

Micrometeorological Data
The flux towers continuously collected data of incident (Rgi) and reflected (Rgr) solar radiation, net radiation (Rn), soil heat flux (G), air temperature (Ta), relative humidity (RH), and wind speed (u) from 2013 to 2016. The sensors and their installed heights and the used data acquisition system in the towers are shown in Table 1.

Micrometeorological Data
The flux towers continuously collected data of incident (Rgi) and reflected (Rgr) solar radiation, net radiation (Rn), soil heat flux (G), air temperature (Ta), relative humidity (RH), and wind speed (u) from 2013 to 2016. The sensors and their installed heights and the used data acquisition system in the towers are shown in Table 1. Table 1. Description of the equipment used to measure incident solar radiation (Rgi), reflected solar radiation (Rgi), net radiation (Rn), soil heat flux (G), air temperature (Ta), relative humidity (RH), wind speed (u), datalogger, and their respective heights in the Fazenda Miranda (FMI) and Baía das Pedras (BPE) flux towers.

Variable
Equipment Description The SEBFs and ET at the two flux towers were calculated using the Bowen ratio energy balance (BREB) method using the sensor listed in Table 1. This method has been widely applied and has the advantage of requiring few micrometeorological parameters  The SEBFs and ET at the two flux towers were calculated using the Bowen ratio energy balance (BREB) method using the sensor listed in Table 1. This method has been widely applied and has the advantage of requiring few micrometeorological parameters while having a firm physical basis [1,39]. In addition, comparisons between estimates obtained by the BREB and the more direct eddy covariance method provide similar data, which makes the MRB an excellent method for environmental studies in remote and logistically challenging areas, such as the Cerrado-Pantanal ecotone [1,39]. The calculation of the SEBFs and ET is described in detail in [1].

Remote Sensing Data
The study was carried out with data and images obtained between 2013 and 2016 using 27 images of surface reflectance and brightness temperature from the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS) sensors, respectively, from Landsat 8 in path 226 and row 71, and 10 images of surface reflectance of the MOD09A1 product from the MODIS sensor on the TERRA satellite were downloaded from the EROS Science Processing Architecture (ESPA) [espa.cr.usgs.gov accessed on 25 April 2020] of the US Geological Survey (USGS).
The OLI sensor images are composed of 9 bands, with spatial resolutions of 30 m for bands 1-7 and 9, and 15 m for band 8 (panchromatic). The images from the TIRS sensor are composed of bands 10 and 11, with spatial resolution of 90 m. The temporal resolution of the Landsat 8 satellite is 16 days and the radiometric resolution is 16 bits [42]. The images of the surface reflectance without the effect of the atmosphere were processed by the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) hosted on the ESPA platform. LEDAPS is a complex algorithm that integrates internal sensor data (metadata) with external data (NCEP, NOAA, and NASA) to (i) transform the digital number to top of atmosphere (TOA) reflectance; (ii) detect pixels with clouds from TOA reflectance and; and (iii) calculate the corrected surface reflectance from the TOA reflectance [43]. The atmospheric correction of the surface reflectance by the LEDAPS was performed with the radioactive transfer code 6 s (Second Simulation of a Satellite Signal in the Solar Spectrum) [44], integrating (i) meteorological data from the NCEP; (ii) digital elevation models of the GCM (Global Climate Model); (iii) internal aerosol optical thickness (AOT); and (iv) ozone data collected by NASA [42,43,45]. LEDAPS also uses the digital elevation model to correct the parallax error due to the local topographic relief, as well as systematic geometric and precision corrections using surface control chips [42,43,45].
The MOD09A1 surface reflectance product of the MODIS sensor is composed of 7 bands of surface reflectance images with spatial resolution of 500 m, temporal resolution of 8 days, and radiometric resolution of 16 bits. The composition of the images allows the observation of the earth's surface every 8 days due to high spatial coverage, low view angle, the absence or shadow of cloud, and the presence of aerosols [46]. The MOD09A1 product is equivalent to measurements at ground level with no scattering or atmospheric absorption. The product algorithm MOD09A1 corrects the effects of dispersion and absorption of gases and aerosols (atmospheric correction), as well as the adjacency effects caused by the variation of land cover, bidirectional reflectance distribution function (BRDF), and the effects of atmosphere coupling and cloud contamination. The atmospheric correction of this product was also performed by the 6 s algorithm, in which data of ozone concentration, water vapor, and aerosols were obtained from other MODIS products and auxiliary products were obtained from NASA's Data Assimilation Office [46]. The reflectance images of the MOD09A1 product surface used in this study were obtained on the same days, or at the most ±2 days than those obtained by Landsat 8, provided there was no precipitation. A surface albedo (α sup ) model for the OLI Landsat 8 was developed in this study using a multiple linear regression of surface reflectance bands ( Figure 2). The α sup model was based on combining MOD09A1 surface albedo (α MODIS ) with OLI Landsat 8 surface reflectance over different land surface cover types. The α MODIS was used as the dependent variable and surface reflectance data from the OLI Landsat 8 were used as independent variables in the multiple linear regression equation.

2.4.1.
Using Landsat 8 (OLI) A surface albedo ( ) model for the OLI Landsat 8 was developed in this study using a multiple linear regression of surface reflectance bands ( Figure 2). The model was based on combining MOD09A1 surface albedo ( ) with OLI Landsat 8 surface reflectance over different land surface cover types. The was used as the dependent variable and surface reflectance data from the OLI Landsat 8 were used as independent variables in the multiple linear regression equation.
where to are the MOD09A1 surface reflectance for bands 1 to 7, respectively. The surface reflectance images from the OLI Landsat 8 were resampled from 30 to 500 m, to have images with a spatial resolution that is consistent with those of . The

A Conventional Model
Surface albedo ( ) was also estimated using a conventional model (Equation (2)) that was used in a number of studies (e.g., [47,48]). This model has been widely applied in environmental studies and in the estimation of SEBFs and ET algorithms, such as SE-BAL [14]. It consists of a simplified radiative transfer equation that has not been evaluated in complex transition regions, such as the Cerrado-Pantanal ecotone. The surface albedo based on this model can be estimated as: where is the planetary albedo; is the albedo of the atmosphere, which is generally assumed to be about [48]; and is the atmospheric transmittance to global solar radiation, calculated by Equation The α MODIS in this study was estimated following the approach of Liang et al. [17], as explained in Equation (1): where ρ 1 to ρ 7 are the MOD09A1 surface reflectance for bands 1 to 7, respectively.
The surface reflectance images from the OLI Landsat 8 were resampled from 30 to 500 m, to have images with a spatial resolution that is consistent with those of α MODIS . The

A Conventional α con Model
Surface albedo (α con ) was also estimated using a conventional model (Equation (2)) that was used in a number of studies (e.g., [47,48]). This model has been widely applied in environmental studies and in the estimation of SEBFs and ET algorithms, such as SE-BAL [14]. It consists of a simplified radiative transfer equation that has not been evaluated in complex transition regions, such as the Cerrado-Pantanal ecotone. The surface albedo a con based on this model can be estimated as: where α toa is the planetary albedo; α atm is the albedo of the atmosphere, which is generally assumed to be about [48]; and τ oc is the atmospheric transmittance to global solar radiation, calculated by Equation (3) [15]: where P o is the local atmospheric pressure (kPa); K t is the atmospheric turbidity coefficient (K t = 1 if clear sky and K t = 0.5 if cloudy sky; we used K t =1); and W is the precipitable water (in mm; see Equation (4)), obtained by the vapor pressure of water (e a ; in kPa): The albedo of the atmosphere, α toa , was calculated following Equation (5) as a linear combination of the top of atmosphere (TOA) reflectance of the OLI Landsat 8 [48], as: where ρ 2 to ρ 7 are top of atmosphere reflectance of bands 2 to 7 of the OLI Landsat 8.

Surface Temperature (T s ) Correction Models
The surface temperature (T s ) was estimated using four currently available models that include: (i) the atmospheric correction parameter calculator (ATMCORR); (ii) the single-channel (SC); (iii) the radioactive transfer equation (RTE); and (iv) the multichannel split-window (SW). These models aim to recover the radiance attenuated by atmospheric constituents in the spectral window between 10 and 13 µm.

T s Barsi Correction Based on ATMCORR
The ATMCORR (atmcorr.gsfc.nasa.gov, accessed on 10 August 2021) is an initiative by NASA to provide a comprehensive atmospheric correction tool for surface temperature [29,36]. ATMCORR integrates data from the National Center for Environmental Prediction (NCEP) that models the global atmospheric profile for certain dates using the well-known MOD-TRAN v4 code in a set of integration algorithms [29]. The atmospheric profiles generated by the NCEP integrate data from satellites and surface data to model the global atmosphere at 28 altitudes in a spatial grid of 1 • × 1 • . The profile data is generated every six hours with the possibility of resampling the grids. The interpolated data from the NCEP is inserted in MODTRAN v4 and the atmospheric parameters are extracted from the MODTRAN output files, adjusting the data for the moment of the satellite's passage. Due to the robust integration of ATMCORR, this model has been widely applied in studies that demand corrected temperature [49,50]. Thus, the surface temperature obtained using the ATMCORR model as described in Barsi et al. [29], referred to in the present study as T s Barsi , was used as a reference to evaluate the T s as obtained by the other three temperature correction models. The T s Barsi (K) can be calculated using Equation (6) as: where K 1 = 607.76 W m −2 sr −1 µm −1 and K 2 = 1260.76 W m −2 sr −1 µm −1 are calibration constants of the thermal band provided by the TIRS Landsat 8 sensor; and L c is the radiance of a blackbody target of kinetic temperature (W m −2 sr −1 µm −1 ; see Equation (7)): where L TOA is the space-reaching or TOA radiance measured by the TIRS (W m −2 sr −1 µm −1 ); ε is the surface emissivity over TIRS band calculated by Equation (8) [51] and the parameters obtained by ATMCORR; L u is the upwelling or atmospheric path radiance (W m −2 sr −1 µm −1 ); L d is the downwelling or sky radiance (W m −2 sr −1 µm −1 ); and τ is the thermal atmospheric transmission.
where ε s and ε v denote bare soil and vegetation emissivity, respectively, over the TIR band; and FVC is the fraction of vegetation cover (Equation (9)): where NDV I is the normalized difference vegetation index; and NDV I min and NDV I max are the minimum and maximum NDV I, respectively, extracted from the NDVI histogram.

T s SC Correction Based on the Single-Channel (SC) Model
The single-channel (SC) model consists of the correction of surface temperature (T s SC ; see Equation (10)) based on correction functions γ, ψ, and δ that can be estimated by the parameters L u , L d , and τ [34]. The SC model can be applied to any of the bands in the Landsat 8 TIRS. This study used band 10 to correct T s [referred to as T s SC ]: where ψ 1 , ψ 2 , and ψ 3 are atmospheric correction functions, calculated by Equations (11)-(13) from parameters obtained from ATMCORR; and γ and δ are functions of L toa , brightness temperature (T b ; K), and b γ , which is equal to K 1 [51]:

T s RTE Correction Based on the RTE Model
The corrected T s using the radiative transfer equation is referred to in this article as T s RTE (K), and was calculated following Equation (16) based on the L toa and the parameters obtained by ATMCORR [51]: where C 1 = 1.19104 × 10 8 W µm 4 m −2 sr −1 and C 2 = 14387.7 µm K are constant; and λ is the effective wavelength of the band.

T s SW Correction Based on the Split-Window (SW) Model
The split-window surface temperature correction model is one of the simplest techniques, in which the radiation attenuation by atmospheric absorption is proportional to the difference in radiance measured simultaneously by the two thermal bands [28,34]. The surface temperature (T s SW ; K) based on the SW model can be calculated as: where T b10 and T b11 are the brightness temperature of bands 10 and 11 (K) of TIRS; c x is constant with the following values c 0 = −0.268, c 1 = 1.378, c 2 = 0.183, c 3 = 54.30, c 4 = −2.238, c 5 = −129.20, and c 6 = 16.40 [34]; ∆ε is the difference in emissivity of the thermal bands 10 and 11 of TIRS; and w is the water vapor concentration (g cm −2 ) calculated by Equation (18) [52].

Estimation of SEBFs and ET Using SEBAL
The SEBAL algorithm was processed according to the flow chart shown in Figure 3. It was proposed to estimate the daily evapotranspiration (ET) from the instantaneous latent heat flux (LE; W m −2 ) obtained as a residue of the energy balance equation (Equation (18)): where Rn is the net radiation (W m −2 ); G is the soil heat flux (W m −2 ); and H is the sensible heat flux (W m −2 ).

Estimation of SEBFs and ET Using SEBAL
The SEBAL algorithm was processed according to the flow chart shown in Figure 3. It was proposed to estimate the daily evapotranspiration (ET) from the instantaneous latent heat flux ( ; W m −2 ) obtained as a residue of the energy balance equation (Equation (18)): where is the net radiation (W m −2 ); is the soil heat flux (W m −2 ); and is the sensible heat flux (W m −2 ).

The
(Equation (19)) represents the balance of short-wave and long-wave radiation on the surface: where ↓ is the measured incident solar radiation (W m −2 ); is the surface albedo; ↓ is the long-wave radiation emitted by the atmosphere in the direction of the surface (W m −2 ); ↑ is the long-wave radiation emitted by the surface to the atmosphere (W m −2 ); and ε is the surface emissivity. The ↑ and ↓ were calculated by Equations (20) and (21): where and are the surface and atmosphere emissivity; are the Stefan-Boltzmann constant ( = 5.67.10 −8 W m −2 K −4 ); is the surface temperature (K); and is the air temperature (K). The ↑ was calculated using the surface temperature calculated by the models described in item 2.6.
The was calculated by Equation (22) [12]: The Rn (Equation (19)) represents the balance of short-wave and long-wave radiation on the surface: where R s↓ is the measured incident solar radiation (W m −2 ); α is the surface albedo; R L↓ is the long-wave radiation emitted by the atmosphere in the direction of the surface (W m −2 ); R L↑ is the long-wave radiation emitted by the surface to the atmosphere (W m −2 ); and ε is the surface emissivity. The R L↑ and R L↓ were calculated by Equations (20) and (21): where ε sup and ε atm are the surface and atmosphere emissivity; σ are the Stefan-Boltzmann constant (σ = 5.67.10 −8 W m −2 K −4 ); T s is the surface temperature (K); and T a is the air temperature (K). The R L↑ was calculated using the surface temperature calculated by the models described in item 2.6. The G was calculated by Equation (22) [12]: where T s is the surface temperature (K) calculated by the different models described in Section 2.6; α sup is surface albedo calculated by the models described in Sections 2.4 and 2.5; NDV I is the normalized difference vegetation index; and Rn is the net radiation calculated by the different T s models described in Section 2.6 and α sup described in Sections 2.4 and 2.5.
H is the central variable in the SEBAL algorithm and estimated by the classic aerodynamic (Equation (23)) [8]: where ρ is the specific air mass (kg m −3 ); c p is the specific heat of air at a constant pressure (1004 J kg −1 K −1 ); dT is the temperature difference near the surface (K); and r ah is the aerodynamic resistance to the transport of sensible heat flux (s m −1 ) between two heights (z 1 = 0.1 m and z 2 = 2.0 m). The r ah is obtained as a function of the friction speed after an iterative correction process based on atmospheric stability functions [8].
The dT was calculated from a linear relationship with the T s (Equation (24)), and the values of the coefficients "a" and "b" were obtained using information from two "anchor" pixels [8]: In SEBAL, the "anchor" pixels represent conditions of hydrological extremes, in which "cold" represents surfaces with H close to zero and "hot" surfaces with LE close to zero. In general, the cold pixel can be represented by a body of water or a well-irrigated crop, and the hot pixel can be represented by a severe surface water restriction, such as exposed soils [8].
In non-agricultural environments, as those of concern in this study, the conditions for choosing the cold pixel may not be properly satisfied, restricting the choice of the cold pixel in areas of native forest. In this study, an approach similar to that used in METRIC was used, using the values of Rn and G of the cold pixel of a known surface and the actual evapotranspiration (ETr) from an estimate reference evapotranspiration (ETo), with local weather station data and the cultivation coefficient (Kc) of the cold pixel surface [15]. Then, the ETr was converted to LE to obtain the H of cold pixel. Thus, it was possible to find the coefficients of Equation (24) and solve the dT by the system formed by Equations (23) and (24) in an iterative process.
After obtaining the LE of each pixel by Equation (18), the daily evapotranspiration (ET; mm d −1 ) of each pixel was calculated by Equation (25), from the instantaneous evaporative fraction (FE i ; see Equation (26)) and daily Rn (Rn 24h ; W m −2 ) of each pixel and the latent heat of vaporization of water (λ; kg m −3 ) [12]:

Evaluation Approach and Performance Indicators
This study followed four steps to evaluate the effects of surface albedo and temperature models on SEBFs and ET that include:

1.
Developing a surface albedo model by combining MODIS and Landsat 8 dataset. A subset of the data was used for model development and the remaining was used to evaluate the model performance over different land cover types. In this analysis, the MODIS surface albedo by Liang et al. [17] was assumed to be as a reference against which to compare the developed and existing models.

2.
Comparing the performance of the of the developed surface albedo model with the currently used conventional model.

3.
Retrieving and evaluating land surface temperature based on four different methods.
In this analysis, the model by Barsi, et al. [29] was assumed to be the reference against which to compare other retrieval methods. The comparison between the different retrieval methods was conducted over the sample sites.

4.
Evaluating the combined effects of the surface albedo models and the brightness temperature and temperature retrieval methods on SEBFs and ET. Since both variables (i.e., α and T s ) are used in SEBAL model to estimate SEBFs and ET, a set of combinations of the two variables were developed as shown in Table 2 to identify these effects. The averages of all variables were calculated with a confidence interval (CI) of ±95% using bootstrapping of 1000 iterations of random resamples with substitution [54]. The accuracy of surface albedo models analyzed in this study as well as the estimated SEBFs and ET were assessed using the Willmott coefficient (d; see Equation (27)), the root mean square error (RMSE; see Equation (28)), the mean absolute error (MAE; see Equation (29)), the mean absolute percentage error (MAPE; see Equation (30)), and the Pearson's correlation coefficient (r):
A comparison of the surface albedo between a MODIS and a sup as well as between a MODIS and a con indicated that a sup performed better than a con , as shown in Table 3. The summary of the comparison shown in Table 2 was based on surface albedo values from all selected sites. The average of a sup was not significantly different from that of a MODIS , while the average of a con was 49% higher than the that of a sup ( Table 3). The RMSE of a sup was 5.6-fold lower and the Willmott and correlation coefficients were approximately 2-fold higher for α sup than a con . Table 3. Average (±95% confidence interval) of the surface albedo estimated by MODIS (a MODIS ) used as reference values, and the average (±95% confidence interval), mean absolute error (MAE), mean absolute percent error (MAPE, %), root mean square error (RMSE), Willmott coefficient (d), and Pearson correlation coefficient (r) of the surface albedo estimated by the model developed in this study (a sup ) and the surface albedo estimated by the conventional model (a con ). Values with (***) indicate p-value < 0.001. All units are dimensionless. Regarding the performance of a sup over the different land use types, it appears that a sup had better performance than a con over the different sampled land uses. The averages a sup and a MODIS were similar in pasture and urban areas, and they were close in the forest and water bodies, while the means of a con were from 36% to 64% higher than a MODIS (Table 4). Table 4. Average (±95% confidence interval) of the surface albedo estimated by MODIS (a MODIS ), used as reference values, surface albedo estimated by the model developed in this study (a sup ) and surface albedo estimated by the conventional model (a con ) in agriculture, urban area, forest, and water bodies on the study area. All units are dimensionless.

T s Retreival Models
Based on a comparison with T s barsi , the results indicated that T s SC and T s RTE had much lower discrepancies based on the obtained MAE, MAPE, and RMSE, and higher agreement based on the Willmott coefficient (d) and Pearson correlation (r), compared to T s SW and T b (Figure 4 and Table 5). The averages of T s barsi , T s SC , T s RTE , and T s SW were not significantly different; however, T b was lower than T s barsi by about 2%. The largest correction error was observed when comparing T b with T s barsi , while T s RTE had the least errors compared to T s barsi . The surface temperatures (T s ) corrected by the different models had MAE and RMSE up to 86% lower than the T b . much lower discrepancies based on the obtained MAE, MAPE, and RMSE, and higher agreement based on the Willmott coefficient ( ) and Pearson correlation ( ), compared to and (Figure 4 and Table 5). The averages of , , , and were not significantly different; however, was lower than by about 2%. The largest correction error was observed when comparing with , while had the least errors compared to . The surface temperatures ( ) corrected by the different models had and up to 86% lower than the .

SEBFs and ET Estimates Based on α and T s Combinations
A summary of the comparison between estimated and measured Rn based on all model combinations (Table 2 over both flux towers, i.e., FMI and PBE) is shown in Table 6.  (Table 2) of a sup with all T s as well as the combination of a con with T b did not show a significant difference from the values measured at the flux towers, but the average of estimated Rn based on the combination of a con and all T s were 15% lower than the measured Rn ( Table 6). The estimated Rn with the combination of a sup and T b had the lowest errors and the highest Willmott's d and r, while the highest errors and lowest coefficient d and r were observed with the combination of a con and T s SW (Table 6). Table 6. Average (±95% confidence interval) of the measured net radiation (Rn; W m −2 ) in the flux towers, and the average (±95% confidence interval), mean absolute error (MAE), mean absolute percent error (MAPE), root mean square error (RMSE), Willmott coefficient (d), and Pearson correlation coefficient (r) of the estimated net radiation using the conventional (a con ) and parameterized (a sup ) surface albedo models combined with brightness temperature (T b ) and the surface temperature corrected by the Barsi model (T s barsi ), the single-channel model (T s SC ), the radiative transfer equation model (T s RTE ), and the split-window model (T s SW ). Values with (***) indicate p-value < 0.001. Unlike Rn, the averages of estimated G based on a sup with all T s retrieval methods, including T b , did not differ between each other, but were between 35% and 54% higher than the measured G ( Table 7). The values of d and r changed significantly, but the errors in estimated G with T b were 18% less than with T s . The average of estimated H based on all combinations of a sup and T s as well as the combination of a con with T b did not show a significant difference from those of the measured values, while the averages of estimated H based on a con and the different T s were between 26-35% lower than those of the measured values ( Table 8). The MAE and RMSE in estimating H with a sup were between 7-47% less than those based on a con . The estimated H with the combination of a sup and T s barsi had the smallest MAE, MAPE, and RMSE, while the largest MAE, MAPE, and RMSE, and smallest d and r were obtained with the combination of a con and T s SW .
Opposed to what was observed with the Rn, G and H, there was no difference in the averages of estimated LE and ET based on the different combinations of a sup and T s (Tables 9 and 10). It should be noted that the MAE, MAPE, and RMSE in LE and ET estimates with a sup were on average 28% and 20% lower, respectively, and the coefficients were slightly higher than those estimated with a con , with emphasis on the combination of a sup and T s barsi . Table 7. Average (±95% confidence interval) of the measured soil heat flux (G; W m −2 ) in the flux towers, and the average (±95% confidence interval), mean absolute error (MAE), mean absolute percent error (MAPE), root mean square error (RMSE), Willmott coefficient (d), and Pearson correlation coefficient (r) of the estimated soil heat flux using the conventional (a con ) and parameterized (a sup ) surface albedo models combined with brightness temperature (T b ) and surface temperature corrected by the Barsi model (T s barsi ), the single-channel model (T s SC ), the radiative transfer equation model (T s RTE ), and the split-window model (T s SW ). Values with (**) indicate p-value < 0.01.     Table 9. Average (±95% confidence interval) of the measured latent heat flux (LE; W m −2 ) in the flux towers, and the average (±95% confidence interval), mean absolute error (MAE), mean absolute percent error (MAPE), root mean square error (RMSE), Willmott coefficient (d), and Pearson correlation coefficient (r) of the estimated latent heat flux using the conventional (a con ) and parameterized (a sup ) surface albedo models combined with brightness temperature (T b ) and surface temperature corrected by the Barsi model (T s barsi ), the single-channel model (T s SC ), the radiative transfer equation model (T s RTE ), and the split-window model (T s SW ). Values with (***) indicate p-value < 0.001.

Surface Albedo Models Performance
The surface albedo model (a sup ) developed in this study performed well compared to the conventional one (a con ). The RMSE based on a sup was less than 0.03 required by climate forecasting models [18] and within the range of 0.01-0.02 found in previous studies [18,56]. The largest discrepancies shown by a con as indicated in the reported MAE and RMSE can be due a number of factors that include (i) the broad spectrum or broadband transmittance being inadequate for the atmospheric correction of the composition of discrete bands; (ii) not considering the differences in atmospheric transmissivity for each band; and (iii) the non-correspondence of the narrow and wide bands with solar radiation on the surface [14].
The a sup in the four land types was within the range found in other studies. The a sup over agricultural areas ranged from 0.14 to 0.18 [56,57]; from 0.15 to 0.20 over urban areas due to the complexity of mixtures of built-up area and vegetation in backyards and streets [58]; from 0.11 to 0.13 over the Cerrado forest [57,59]; and from 0.05 to 0.07 over water bodies depending on the composition of the water [58,59]. The a con of water bodies was greater than 0.10, which is above the values obtained in the lakes of the region [59].

Evaluation of T s Retrieval Models
In general, the difference between T b and T s varies between 1 and 5 K in the 10-12 µm spectral region, subject to the influence of atmospheric conditions and surface emissivity [60]. For mathematical convenience, the TOA thermal radiation is generally expressed in terms of T b with an emissivity of 1.0 [61]. The TOA radiance is the result of radiation emitted by the Earth's surface, upward radiation emitted by the atmosphere, and downward radiation emitted by the atmosphere [62]. The TOA radiation is mostly attenuated by water vapor and, to a lesser extent, by trace gases and aerosols [63,64].
In this study, T s barsi was used as a reference to validate T b , T s SC , T s RTE , and T s SW , because previous studies indicated T s barsi based on Landsat had relatively low MAE and RMSE values low (ranging between 0.2 and 2.5 K) when compared to the field measurements of T s that were taken over different land surface types represented by the Surface Radiation Budget (SURFRAD) Network (https://gml.noaa.gov/grad/surfrad/ accessed on 10 August 2021) during different atmospheric profiles [65][66][67][68]. Furthermore, the improved performance of T s barsi is associated with the robustness of the MODTRAN algorithm and the integration of atmospheric data from the NCEP to generate the atmospheric correction parameters (L u , L d , and τ) [36]. These parameters estimated by ATMCORR were also used in the other models, which may justify the good relationship between the T s estimated by the T s SC , T s RTE , and T s SW .
The good relationships between T s SC , T s RTE , and T s SW with T s barsi obtained in this study agreed with other validation and simulation studies, which indicated that the MAE and RMSE obtained in this study are within those limits reported in the literature. The typical MAE and RMSE of T s SC and T s RTE vary between 1 and 3 K [31,69], and the T s SW is around 1.5 K [33]. Using low spatial resolution data, T s SC and T s RTE presented MAE and RMSE from 1.6 to 2.4 K [70], and T s SW from 1.5 to 2.9 K [71].
The good agreement of T s RTE with T s barsi maybe due to both models using the radiative transfer equation of Planck's inverse equation [29,30,35,51]. The main difference of T s RTE and T s barsi is on the conversion of thermal radiance into T s , since T s RTE is converted by the inverted Plank equation and T s barsi by a specific Planck curve equation with calibration constants determined for the TIRS Landsat 8 [35,36]. T s RTE has been widely used in studies of water bodies with an accuracy of around 0.2 K and in studies of terrestrial bodies with errors of up to 2 K [35,72].
The RMSE of T s SC around 1.3 K showed its good agreement with T s barsi , at the lower limit of the range from 1.2 to 2 K obtained under different conditions of atmospheric water vapor [30,34]. The biggest errors of T s SW can be attributed to the model being multichannel, which introduces greater noise if using only one thermal channel [28,34,73]. However, T s SW is obtained by combining thermal bands with defined coefficients, considering different emissivity for each band and requiring only knowledge of the atmospheric water vapor [28,34].

The Effects of α and T s Retreival Models on SEBFs and ET
In general, RMSE of Rn is typically found to be between 20 and 80 W m −2 with different orbital sensors (TM Landsat 5, TM+ Landsat 7, and MODIS) [59,[74][75][76][77][78][79][80]. The RMSE obtained in this study were close to those reported by [59] over the Cerrado zone and by [10] on the Cerrado-Pantanal transitional zone in Brazil, which highlight the relatively acceptable accuracy of estimated Rn obtained based on all combinations. The better performance of the Rn estimated with the T b maybe due to the shortwave and longwave radiation balance [10]. The a sup can be overestimated by up to 15%, which leads to an underestimation of Rn [11,81], while T b is generally lower than T s , leading to an underestimation of long-wave radiation emitted by the surface (R L↑ ), which therefore leads to overestimation of Rn. Despite the better performance of Rn with T b , the MAPE of Rn estimated with a sup and all T s were less than 2%, and the RMSE less than 20 W m −2 . In addition, the difference in MAE and RMSE of the estimated Rn with all T s and the same surface albedo model was less than 5 W m −2 and MAPE less than 1%.
The obtained MAE and RMSE values of G were within the range of 15-32 W m −2 , which was similar to those obtained in other studies [82,83]. The low performance of G has been reported in other studies with different land uses [82][83][84]. Probably, the low performance of the G estimate is due to the low sensitivity of the model to the high spatial complexity of the study area. G tends not to have a high impact on the SEB and ET of densely vegetated surface, due to the lesser part of the available energy used to heat the soil, but G tends to impact the SEB and ET of surfaces with low vegetation cover, as the pastures and some natural grasslands in Cerrado and Pantanal [13,82,85].
The MAE and RMSE of H estimated based on all combinations with a sup and the combination of a con with T b were less than the 50 W m −2 that was reported by [6,78,82,86].
Estimates of H with T s barsi , T s SC , and T s RTE were on average 3% lower than that with T b , indicating that the differences between T s and T b do not significantly impact H. This is because the internal calibration process of SEBAL alleviates impacts of low T b values [11]. The estimation of H by SEBAL is a function of the linear relationship "dTs = a + bT s ", using two extreme pixels to calculate the constants "a" and "b" [8,15]. The initial value of these constants is obtained from meteorological information, satellite estimates (Rn − G; SAV I), and the operator's choice (anchor pixels), and these constants are adjusted by iterations [11,15,87]. The estimation of these constants by numerical iterations eliminates the effects of the negative bias of T b and transmits the calibration effect for all other pixels in proportion to the inserted T s . Therefore, the differences between T s and T b tend not to significantly affect SEBF estimates [11,16,87].
The MAE and RMSE of LE estimates were within the range of 30-70 W m −2 found in previous studies based on measurements with flux towers and lysimeters, and the MAPE was less than 20% [78,82,84,86,88]. The MAE and RMSE of the ET estimates were also within the range of 0.3 mm to 0.6 mm day −1 , and the MAPE within the range of 8% to 20% found in other studies [6,78,[89][90][91]. In this study, SEBAL was applied in areas with grasses and shrubland typical of the Cerrado-Pantanal transition region under different natural water conditions and obtained errors between 11% and 12.5%, which represented absolute errors of ET less than 0.35 mm day −1 .
The slight difference between the MAE, MAPE, and RMSE and the correlation and Willmott coefficients of the LE and ET estimated with T b and T s shows that the recovery of T s by the models does not significantly impact the estimation of these parameters. This effect was also observed in studies by [11] and [16]. This reinforces that the internal calibration of SEBAL keeps the "dTs = a + bT s " stable and minimizes the impacts of the insertion of T s in the LE and ET estimates, since the T s of the anchor pixels represent the extreme conditions of water availability, regardless of the removal of the effect of the atmosphere and the emissivity of the surface on the thermal band [11,87]. In contrast, LE and ET performed better with a sup instead of a con . A similar result was also observed in the work of [11], which proposed an internal SEBAL calibration to remove the effect of the atmosphere in each band of the Landsat 5 sensor, whose estimate of ET with a con introduced random errors of ±1 mm day −1 .

Conclusions
In this study, a model of surface albedo (a sup ) with the OLI Landsat 8 surface reflectance was developed. Surface temperature (T s ) was recovered by different models from the brightness temperature (T b ). The performance of surface energy balance fluxes (SEBF) and evapotranspiration (ET) estimates, based on different combinations of surface albedo and temperature models, was evaluated against ground-based observations of SEBF. The a sup model performed better than a conventional surface albedo model (a con ) as it provided lower MAE, MAPE, and RMSE and higher Willmott coefficients (d) and Pearson correlation (r) when compared with surface albedo data based on MODIS (a MODIS ). In addition, average values of a sup were similar to those found by a MODIS , while those of a con were about 36-64% higher than a MODIS . Additionally, a con showed some limitations over water bodies. Minimizing these errors in spatially complex areas, such as the Cerrado-Pantanal transition, is important for accurate estimates of SEBFs and ET.
The retrieval of surface temperature (T s ) by the different models combined with a con significantly influenced estimates of the net radiation (Rn) and the sensible heat flux (H). Estimates of the Rn were on average 15% lower and those of H, which were about 26-35% lower than the measured Rn and H, respectively. However, estimates of Rn and H based on the combination of T s with a sup were not significantly different from those measured. Moreover, the averages of latent heat flux (LE) and evapotranspiration (ET) were also not significantly different from those measured based on all combinations.
The determination of the a sup model, with the OLI Landsat 8 surface reflectance for the studied Cerrado-Pantanal transition region, improved the performance of SEBAL in estimating the Rn, H, LE, and ET, when combined with both T s and T b . SEBFs and ET estimated by SEBAL with a sup had lower errors (i.e., RMSE) and higher agreement and correlation coefficients d and r. It is noteworthy that the SEBFs and ET estimated by the combination a sup and T s barsi presented the best performance. The combination of a con and T s SW worked well to estimate ET over the mixed shrub-grass site of the PBE, while combination of a sup and T b worked well to estimate ET over the grassland site of the FMI. The evaluation conducted in this analysis over the spatially complex gradient of natural ecosystems in southern Brazil provided a robust test of the performance of these surface albedo and temperature algorithms and can help to guide future studies on the use of appropriate models for the estimation of SEBFs and ET over other regions with similar complex environments.