Assessing uncertainties in estimating surface energy fluxes from remote sensing over natural grasslands in Brazil

Evapotranspiration (ET) is one of the main fluxes in the global water cycle. As the Brazilian Pampa biome carries a rich biodiversity, accurate information on the ET dynamics is essential to support its proper monitoring and establish conservation strategies. In this context, we assessed an operational methodology based on the Simplified Surface Energy Balance Index (S-SEBI) model to estimate energy fluxes over the natural grasslands of the Pampa between 2014 and 2019. The S-SEBI is an ET model that requires a minimum of meteorological inputs and has demonstrated reasonable accuracy worldwide. Therefore, we investigated the model performance considering radiation data from both ERA5 reanalysis and Eddy Covariance measurements from a flux tower. Furthermore, comparisons from satellite-based estimates with in situ measurements were performed with and without energy balance closure (EBC). Results indicated that the meteorological inputs have low sensitivity on daily ET estimates from the S-SEBI model. In contrast, the instantaneous energy balance components are more affected. The strong seasonality impacts the evaporative fraction, which is more evident in late summer and autumn and may compromise the performance of the model in the biome. The effects in the daily ET are lower when in situ data without EBC are considered as ground truth. However, they are less correlated with the remote sensing-based estimates. These insights are useful to monitor water and energy fluxes from local to regional scale and provide the opportunity to capture ET trends over the natural grasslands of the Pampa.


Introduction
Evapotranspiration (ET) is widely used to measure the amounts of total water loss through several key processes between land and the atmosphere (Wang and Dickinson 2012;Dou and Yang 2018). Its prediction plays an important role in drought analysis, climate change studies, water level balance, agricultural and forest meteorology, long-term decision-making in food and water security policies, and optimum allocation of water resources (Valipour et al. 2019).
ET can be directly measured by using either lysimeter or water balance approach (Kumar et al., 2011). Nonetheless, these methodologies do not allow estimating the land surface fluxes when dealing with large spatial scales (Chen et al. 2005;Courault et al. 2005).
Satellite observations have been used for monitoring surface conditions over the last few decades because they provide the potential to bridge the gap between point measurements and larger scale surface processes (Ma et al. 2015), thus supplying the need for accessible data for the regions lacking land measurements (Zhang et al. 2018). In addition, accurate land use maps, and abundant meteorological reference data have considerably improved the ability to obtain remotely sensed ET and its flux components (Webster et al. 2016).
Several models have been developed using information from different sensors and often in conjunction with ancillary surface and atmospheric data for ET estimation (Bastiaanssen et al. 1998;Roerink et al. 2000;Senay et al. 2007;Allen et al. 2007). They usually vary from purely empirical to more physically based techniques derived from the surface energy balance (SEB) equation and the information obtained from relationships between satellite-derived vegetation index and land surface temperature (LST) (Liou and Kar 2014). Some of these methods are dependent on many variables measured in surface weather stations as input data, which might prevent its application in large areas with insufficient weather data measuring stations. In contrast, there are some simplified approaches that adopt assumptions of meteorological parameters (Mutti et al. 2019). These simplifications and multiple representation forms may lead to uncertainties in the estimates generated (Chen et al. 2016). Consequently, an evaluation of the benefits and inadequacies of these approaches is a fundamental step.
The Simplified Surface Energy Balance Index (S-SEBI) model (Roerink et al. 2000) estimates ET as a fraction of available energy (named evaporative fraction), in which maximum and minimum temperatures are determined over dry and wet surface lines drawn in an LST-albedo relationship. The model has been applied by several researchers that reported reasonable accuracy (Sobrino et al. , 2007Fan et al. 2007;Olivera-Guerra et al. 2014;Bhattarai et al. 2016). One of the main advantages of using the LST-albedo relationship is that the albedo is sensitive to the total vegetation cover including green and senescent vegetation (Merlin et al. 2014). Furthermore, the model requires a minimum of meteorological data, being entirely suitable for zones with few or none in situ measurements (Olivera-Guerra et al. 2014).
Quantifying ET in a variety of terrestrial ecosystems is valuable for a better understanding of different surface types, assessment of local to global water balances, and water management of agricultural lands. Moreover, the specific features of each biome are essential to formulate proper conservation strategies, appropriate infrastructure, and to understand the impacts of land use changes in each environment (Oliveira et al. 2017;Wagle et al. 2017). Natural grasslands are fundamental to preserve water resources and carbon accumulation in the soil (Moreira et al. 2019), and around the world, they cover vast regions including the North American Great Plains, the Eurasian steppes of Russia, China and Mongolia, and the South American Pampas (Chaneton et al. 2012). In Brazilian territory, the Pampa biome represents 2.07% and is recognized as containing a rich biodiversity characterized by a meadow mosaic, with small scrub vegetation areas and forests (Ruviaro et al. 2016). It is known that different human activities may have distinctive impacts on water use due to economic development (Yu et al. 2021). For conservation purposes, accurate information on its dynamics becomes more relevant to support the proper monitoring of the water and energy fluxes.
Over grasslands and pastures in FL, USA, Battarai et al. (2016) suggested that S-SEBI performance for ET retrieving is comparable to the algorithms surface energy balance algorithm for land (SEBAL) and mapping evapotranspiration at high resolution and with internalized calibration (MET-RIC). Similarly, Senkondo et al. (2019) compared daily ET estimates derived from the models operational simplified surface energy balance (SSEBop), SEBAL and S-SEBI in eastern Africa. The authors reported that the S-SEBI model exhibited a statistically similar ET as the ensemble mean of all models tested and highlighted that the model can be applied elsewhere, especially where observed meteorological variables are limited. S-SEBI model has been applied in a large variety of climate conditions. However, a detailed assessment of the S-SEBI model in estimating surface fluxes over natural grasslands has not been well explored, especially considering the uncertainties associated with input data, assumptions made, and initial conditions. Furthermore, there is no consensus on which model performs better under the particular conditions and climate of the Pampa biome. Therefore, the purpose of this paper is to investigate an operation methodology based on the S-SEBI model to accurately estimate the surface energy fluxes over natural grasslands. We established the main advantages and major uncertainty sources of the method and documented how the ET varies seasonally in the humid subtropical climate of the Brazilian Pampa.

Study site description
The Pampa biome is composed of old grasslands and covers an area of 178,243 km 2 , which includes the whole Uruguay territory, a part of Argentina, and about two-thirds of Rio Grande do Sul State in southern Brazil (Overbeck et al. 2007). The vegetation consists of a mosaic dominated by a grassland matrix with occasional forested areas, and is home to one of the greatest diversities of grassland plants in the world (Boldrini 2009). The experimental area of this study (29°43′27.502″ S; 53°45′36.097″ W; 88 m elevation) is part of the international long term ecological research (ILTER) network and is located in the Federal University of Santa Maria, covering approximately 24 ha of natural vegetation characteristic of the Pampa biome (Rubert et al. 2018) suitable for field validation of different moderate-resolution satellite products (Fig. 1).
In the Pampa biome, cattle production is directly associated with the conservation of its natural grasslands, because it has been carried out for over 300 years and is the central maintainer of Pampa features. These activities are characterized by low environmental impact, with little or no contribution of external input (Viglizzo et al. 2001). Nevertheless, the expansion of the agricultural border together with overgrazing are the most frequent phenomena threatening the Pampa (Oliveira et al. 2017), leading to the loss of soil physical quality and consequent decrease plant diversity (Vargas et al. 2015).
The Pampa biome climate is subtropical humid (Cfa) with hot summers according to the Köppen climate classification system (Maragno et al. 2013). It differs from tropical regions of Brazil in terms of seasonality, which is mainly dominated by solar radiation with regular precipitation rates (Rubert et al. 2018). These characteristics ensure the particular environmental conditions of the Pampa ecosystem. In the summer season, the maximum temperature can be higher than 30 °C, whereas the winter is marked by low temperatures and frosts. The annual average temperature varies between 16 and 18 °C, and the precipitation between 1500 and 1600 mm (well distributed throughout the year). Moreover, this region is constantly subject to sudden changes in weather caused by the passage of the polar front.

Material and methods
According to Chen et al. (2016), there are three major sources of uncertainties related to surface fluxes: inadequate model structure, model input errors, and poorly defined parameters. We focused on studying these aspects for S-SEBI in order to find its fragilities and establish an operational method for retrieving ET over the Brazilian Pampa grasslands. The summary steps of the study are exhibited in the flowchart (Fig. 2).

Theoretical basis and S-SEBI formulation
Among the various energy fluxes within the atmosphere and Earth's surface, the knowledge of both sensible and latent heat fluxes is of fundamental importance for numerical modeling of atmospheric and hydrological processes (Gomez et al. 2005;Liou and Kar 2014). The most applied approaches to estimate ET from remote sensing observations are based on the simplified form of the SEB equation, which is given as: where Rn is the surface net radiation (W m −2 ), G is the soil heat flux (W m −2 ), and H is the sensible heat flux (W m −2 ), and LE is the latent heat flux of evaporation due to ET (W m −2 ). ET in volume units can be computed from LE by the amount of energy needed to evaporate water at a specific temperature and pressure (Glenn et al. 2010). LE is estimated as a residual according to: The Rn can be obtained from remote sensing data according to the following equation (Hurtado and Sobrino 2001): (1) where R s is the shortwave downward radiation (W m −2 ), R ld the longwave downward radiation (W m −2 ), α is the albedo estimated according to Liang (2001), ε is the land surface emissivity (LSE), T s is the LST in Kelvin (K), and σ = 5.67 × 10 −8 W m −2 K −4 is the Stefan-Boltzmann constant. G represents the energy used to heat up the soil or energy released due to cooling of soil and is frequently estimated using an empirical relationship between albedo, LST and normalized difference vegetation index (NDVI) (Bastiaanssen 2000): The LST is a key variable to be retrieved from the TIR data because it reflects the amount of radiation emitted from the surface and subsurface of the earth, and the exchange of energy between the earth surface and atmosphere (Weng et al. 2019). However, atmospheric, angular, and emissivity effects have to be compensated in order to acquire reliable estimates ). The split  Sobrino et al. (1996). Thus, the LST can be computed as where Tj sen and Tj sen are the at-sensor brightness temperatures at the bands i and j (10 and 11 for Landsat 8) in K, ε is the mean LSE, ε = 0.5( εi + εj), Δε is the LSE difference, Δ = ( i − j) and w is the total atmospheric water vapor content (in g·cm −2 ) estimated according to Buck (1981). As a prior knowledge to obtain LST, LSE was calculated from NDVI values according to Sobrino et al. (2008) using the NDVI threshold method. The S-SEBI model does not require calculation of H. On the other hand, the contrast between an albedo-dependent maximum and minimum LST for dry and wet conditions, respectively, is a main base of the model to partition available energy into H and LE fluxes and is used for computing the evaporative fraction (Λ). No additional meteorological data is needed if the surface extremes are available on the scene studied (Liou and Kar 2014). In the case of a wet pixel, H is assumed as zero, so the maximum LE (LE max ) can be estimated by subtracting G from Rn (LE max = R n − G). In contrast, at a dry pixel H is the highest (H max ), which can be estimated by subtracting the G from Rn (i.e., H max = Rn − G). Therefore, Λ is expressed as: where T H and T LE are the temperatures of the dry and wet conditions in K. Afterwards, the daily ET (ET day ) can be computed assuming that Λ is constant for the day under clear-sky conditions (Brutsaert and Sugita 1992), whereas the instantaneous Rn can be extrapolated to a daily scale using the C di ratio. This ratio consists of a relationship between the daily net radiation flux and instantaneous radiation flux (R ni ) and depends on the day of year (DOY) (Sobrino et al. 2007;Olivera-Guerra et al. 2014): where λ is the latent heat of water vaporization (2450 J g −1 ).

In situ flux data
Eddy covariance (EC) method has been used as the typical reference data for validating several ET estimates at the site and pixel level scales (Numata et al. 2017;Li et al. 2018). The data required to estimate ET with EC over the study site were obtained from a flux tower installed at 29.725°S; 53.760°W ( Fig. 1) with the follow sensor: a 3D sonic anemometer (Wind Master Pro; Gill Instruments, Hampshire, UK), to measure wind components and air temperature, and a gas analyzer (LI7500, LI-COR Inc., Lincoln, USA), to measure H 2 O/CO 2 concentrations. Both sensors were installed at 3 m height of the surface and sampled at a 10 Hz frequency, until June 15, 2016. After this period, the gas analyzer and the anemometer were replaced by the sensor integrated CO 2 and H 2 O open-path gas analyzer and a 3D sonic anemometer (IRGASON, Campbell Scientific Inc., Logan, USA), measuring in the same height and frequency. The EddyPro Advanced software (version 5.1, LI-COR) was used to process the LE and H data over half-hour time scales using the EC method (for detailed configurations of EC processing see Rubert et al. 2018). The EC data footprint analyses indicate that about 90% of the flux originated within a circle with a radius of 115 m centered in the flux tower.
Along with the flux tower, the shortwave downward radiation, R s , longwave downward radiation, R ld , and net Previous studies have reported energy balance closure (EBC) problems in EC data due to the systematic biases in instruments, energy sources not considered and losses of turbulent fluxes at high and low frequencies (Twine et al. 2000;Anderson et al. 2008;Tang et al. 2013). A common approach to correct the error is to maintain a constant Bowen-ratio (β) when closing the energy balance (Twine et al. 2000;Wilson et al., 2002;Costa et al. 2010). Thus, the EBC can be forced by using: where H and LE are the sensible and the latent heat flux, respectively (W m −2 ), and Rn is the surface net radiation (W m −2 ). Because there is no consensus about reconciling the surface energy imbalance measured by the EC system , the daily comparisons from satellite data with in situ EC measurements were performed both with and without EBC. To obtain the in situ daily ET (ET in situ with EBC), the corrected LE (LE EBC ) in units of W m −2 was converted to mm/day using the factor 0.0353. We assessed flux data between 2014 and 2019, which included only dates with no clouds over the study site mask. Furthermore, only days with less than 2 h missing LE and H data in situ were considered when Landsat-8 OLI/TIRS products were available. These gaps were filled using a simple interpolation method.

Reanalysis data
Uncertainties of coarse reanalysis data may impact the quality of the derived ET (Vinukollu et al. 2011;Badgley et al. 2015). To evaluate the meteorological inputs for the model, we were acquired from the ECMWF ERA5 reanalysis dataset (Hersbach 2016). The product represents the fifth ECMWF reanalysis generation for the global climate and weather and was obtained at a spatial resolution of 0.25° × 0.25° and hourly temporal resolution. The inputs employed included hourly the average of R s and R ld between 13:00 and 14:00 a.m. GMT time in order to match it with the satellite overpass. To obtain daily values, the hourly ERA5 data were posteriorly converted through the method previously mentioned (see 2.3.1).

Landsat data
The Landsat program has the longest record of Earth observation from space (Pahlevan et al. 2014), providing a 40-year mission of acquiring global, moderate resolution images of the Earth's surface every 16 days (Barsi et al. 2014;Ke et al. 2015). One of the most relevant products derived from the analysis of a combination of the Landsat visible/near-infrared (NIR) and thermal infrared (TIR) bands is the measure of the ET (Reuter et al. 2015). Landsat 8 satellite is the latest member of Landsat family and carries two sensors, the operational land imager (OLI) and the thermal infrared sensor (TIRS). OLI collects data at a 30-m spatial resolution with eight bands located in the visible, NIR, and in the shortwave infrared regions of the electromagnetic spectrum, whereas TIRS measures the thermal radiance at 100 m spatial resolution using two bands located in the atmospheric window between 10 and 12 μm. For the previous satellites (Landsat 4, Landsat 5, and Landsat 7), there was only one TIR channel. Landsat 8 launch in 2013 was a valuable innovation once its two thermal channels and very high sensitivity allow the application of the widely known split-window (SW) algorithm to retrieve LST (Jiménez-Muñoz et al. 2014) thus providing an opportunity to assess the climate-terrestrial interactions through combinations of both thermal and vegetation remote sensing (Webster et al. 2016). We acquired twenty-eighth Landsat 8 scenes (between 2014 and 2019) with clear-sky conditions over the study area from the US Geological Survey website in level 1 (L1) product. Landsat L1 data are radiometric, geometric, and terrain-corrected.
To obtain the normalized difference vegetation index (NDVI) and albedo, Landsat 8 OLI surface reflectance product was also downloaded from the Landsat Data collection. These products are generated at the earth resources observation and science (EROS). The EROS science processing architecture (ESPA) on-demand interface corrects satellite images for atmospheric effects to create level 2 data products. The data are generated from the land surface reflectance code (LaSRC) that uses a unique radiative transfer model (Vermote et al. 2016).

Performance metrics
The image processing was automated through the development of the algorithms in interactive data language (IDL). Ground validation was carried out for the pixel where the tower flux is located (Fig. 2). The criteria employed to assess the performances of the flux estimates were the coefficient of correlation (R), root mean square error (RMSE), and mean absolute error (MAE). The RMSE is defined by the square root of the sum of the variances and describes the accuracy of estimations, which given as: where Pi is a predicted/simulated value, Oi is an observed value, and n is the number of observations. The MAE is the average of the absolute difference between the predictions and the observations, defined as:

Uncertainties of ERA5 reanalysis data
Measurements of energy radiative fluxes are essential in assessing theoretical treatments of radiative transfer in the atmosphere. ET models usually demonstrate different behavior depending on the meteorological input data used (Badgley et al. 2015). In this section, we evaluated the agreement of ERA5 meteorological reanalysis with the in situ R s and R ld measured at the flux tower at the Landsat 8 overpass time (Fig. 3).
R s determines the surface radiative energy balance during the daytime, and significant uncertainties have been reported in its global reanalysis products, which are usually related to the modeling schemes used in the reanalysis systems (Yi et al. 2011). R s followed an apparent seasonal pattern, with higher values in warmer seasons and lower in the colder ones, which is directly related to the seasonal difference of solar azimuth angle (Moriwaki and Kanda 2004).
In general, the errors have a similar range during all seasons in which it was found that ERA5 product underestimated the R s , with RMSE of 69.80 W m −2 and MAE of 62.37 W m −2 (Fig. 3a). Opposite results were found by Tall et al. (2019) that reported that R s from ERA5 tended to be overestimated particularly during monsoon time over Western Africa. . Even though, the authors highlighted that ERA5 product has extensive changes compared to the previous product ERA-Interim (such as higher spatial and temporal resolutions and a generally improved representation), because the bias was slightly reduced in that case. Urraca et al. (2018) also compared estimates from ERA5 and ERAinterim products and observed better R s results with ERA5.
R dl from ERA5 relative to in situ data was also underestimated, but exhibited a less characteristic seasonal pattern compared to R s . In general, the underestimations are similarly distributed through the year, with RMSE and MAE of (11) |Pi − Oi| 29.11 W m −2 and 28.46 W m −2 , respectively (Fig. 3b). It is already documented that the measurement of R dl can be dubious (Sugita and Brutsaert 1993;Duarte et al. 2006), mainly because it is greatly influenced by the cloud coverage and water vapor content. In other words, a low ratio of R s at the surface to clear-sky solar radiation indicates the occurrence of clouds, which also means greater R dl (Yao et al. 2014). Figure 4 illustrates the comparisons of instantaneous G, Rn and LE estimated from the S-SEBI model (S-SEBI ERA5 and S-SEBI in situ ) against ground measurements obtained from EC data. The tower footprint was utilized to evaluate the performance of the model in comparison to the in situ data at the Landsat 8 overpass.

Comparison of instantaneous surface energy components
G component relative to in situ data was overestimated by both satellite-based measurements (S-SEBI-ERA5 and S-SEBI in situ ), and generated very similar correlations. Nonetheless, RMSE and MAE were found to be higher when G was obtained using in situ parameters (S-SEBI in situ ) and produced 54.73 and 52.65, respectively. For S-SEBI ERA5 , RMSE, and MAE had values of 47.40 and 45.11, respectively. In addition, it should be noted that G in situ from EC measurements has very low values, reaching negative values in colder periods. Limitations of the empirical formulation used to calculate G from satellite observations were already mentioned in the literature (Allies et al. 2020;Käfer et al. 2020). The G component is commonly sensitive to changes in vegetation cover because as vegetation grows and shades the surface, less incoming solar energy can reach the ground, attenuating the temperature transfer to the soil (Yang et al. 1999;Zimmer et al. 2020). Especially in the Pampa biome,  cattle production is associated with the conservation of the grasslands (Oliveira et al. 2017) as part of its natural ecosystem. Therefore, besides the climatic seasonality, the Pampa biomass is directly influenced by the different forms of grazing management to which it is subjected.
The biomass is significantly low in winter when there is less solar radiation available, which clearly affects the components used to calculate G and it results in low values. In addition to LST and albedo, Rn is also required in G calculation from remote sensing and this component is dependent on R s and R dl . As the sum of the incoming and outcoming shortwave and longwave radiations, Rn typically achieves its maximum values in summer and minimum in winter (Kofroňová et al. 2019). Rn had a strong correlation for both S-SEBI ERA5 and S-SEBI in situ with values of 0.929 and 0.944, respectively. RMSE and MAE were considerably lower when S-SEBI was running with in situ values (S-SEBI in situ ), producing 42.48 and 22.65, respectively. On the other hand, for S-SEBI ERA5 both metrics are higher, resulting in 96.31 and 89.01, respectively. These results demonstrate that even considering a more accurate Rn, a worse G is produced by the model. In other words, G formulation (Eq. 4) is not really sensitive to the Rn component, and consequently, R s and R dl (Table 1).
LE had correlation coefficients of 0.629 and 0.601 for S-SEBI ERA5 and S-SEBI in situ , respectively. RMSE and MAE were also higher for S-SEBI carried out with in situ data (S-SEBI in situ ), exhibiting values of 93.92 and 77.41 for S-SEBI ERA5 and 118.47 and 101.55 for S-SEBI in situ. These results are in accordance with other studies performed in the Brazilian Pampa (Schirmbeck et al. 2018;) and suggest that the instantaneous LE from S-SEBI model can be more affected by other factors (such as the Λ, obtained from remote sensing through the Eq. 6) than G or Rn, which will be discussed in the next section.

Comparison of daily estimates
The S-SEBI model assumes that Λ is stable during the daytime, which is used as a basis to convert the surface energy fluxes from instantaneous to daily (Brutsaert and Sugita, 1992). We compared the Λ estimated by S-SEBI with Λ obtained in situ considering both scenarios with and without EBC (Fig. 5). Λ in situ was determined from EC measurements by using its general formulation: According to Zhang and Lemeur (1995), the concept of Λ cannot be valid under non-clear sky conditions because the diurnal constancy of Λ may not be fulfilled under cloudy circumstances. Despite the two correlations are weak, regression analysis showed that Λ is much more correlated with ET in situ when EBC is performed. It indicates that even  if the EBC is not properly achieved, using EC data with EBC as ground truth produces better estimates. Zahira et al. (2009) pointed out that errors in the determination of T H and T LE lines, and consequently in Λ (Eq. 6), impact the estimation of LE. Although the use of the relative diagram LST-albedo to estimate the evaporative regimes can be considered a potential strength of the S-SEBI , it requires enough wet and dry pixels in the scene for hydrological contrast. Hence, in areas where the relationship is not well correlated, a more sophisticated algorithm may be necessary. The LST is very important in diagnosing many of the major surface energy balance components, including sensible heat, net radiation, and soil heat flux (Anderson et al., 2012). LST-dependent models assume that ET can cool land surfaces under the condition of homogeneous atmospheric forcing (Sun et al. 2016). Overall, retrieving LST from satellite data with high accuracy has been a research problem at least in the last three decades, especially because of the effects introduced by the atmosphere in the thermal infrared region that must be corrected (Qin et al. 2001;Sobrino et al. 2005;Tardy et al. 2016); otherwise, they may affect the accuracy of both H and LE. Uncertainties in LST retrieval algorithms are generally related to the high w content in the atmosphere, and most of them tend to perform badly under conditions of w > 3 g.cm −2 (Prata, 1994). However, the SW algorithm is able to perform well even in conditions of high-water vapor, as reported by several researchers (Jiménez-Muñoz et al. 2014;Yu et al. 2014), which ensures the reliability of the methodology employed. Little has been discussed about the effects of LST accuracy in the ET estimations. Particularly over the Pampa, Rocha et al. (2020) studied the LST influence in the S-SEBI model between 2015 and 2019 and found that an error of up to 2 K in the LST produces an average variation of 0.18 mm/day. Consequently, the authors concluded that LST itself does not actually influence the average of the retrieved ET (see  for more extensive discussion).
Besides being the dominant controlling factor of climate and hydrology, ET is one of the main fluxes in the global water cycle. Figure 6 depicts the comparison between daily ET between satellite-based measurements derived from the S-SEBI model (S-SEBI ERA5 and S-SEBI in situ ) against field measurements with and without the EBC. Coefficient of correlation was high for both S-SEBI ERA5 and S-SEBI in situ when compared to in situ data considering EBC (Fig. 6a). On the other hand, when contrasted with in situ data without EBC, the correlation did not show significant variations between S-SEBI ERA5 and S-SEBI in situ.
RMSE and MAE had superior agreement with in situ data with EBC than without it, producing 1.38 and 1.12 for S-SEBI ERA5 and 1.06 and 0.79 for S-SEBI in situ , respectively. In contrast, when considering EC data without EBC as ground truth, the values found were 1.71 and 1.30 for S-SEBI ERA5 and 1.44 and 1.06 for S-SEBI in situ , respectively. Furthermore, these performance metrics are notably more sensitive to the components than the coefficient of correlation which did not show much differences. The results found are in agreement with other validation exercises reported in the literature for the S-SEBI model (Roerink et al. 2000;Gomez et al. 2005;Verstraeten et al. 2005;Sobrino et al. 2007;Galleguillos et al. 2011;Käfer et al. 2020). In this Fig. 6 Comparative analysis of the daily evapotranspiration (ET) between satellite measurements (S-SEBI ERA5 and S-SEBI in situ ) and field measurements a with and b without energy balance closure (EBC) context, the Pampa grasslands do not behave so differently from other vegetation types in terms of accuracy of the final daily ET.
Assuming EC data with EBC as ground truth, the RMSE of the daily ET is impacted in 0.32 mm/day if ERA5 radiation data are used instead of in situ ones. It indicates that the radiation components (R s and R dl ) may play a smaller role in the S-SEBI model. These variations are even less pronounced when EC data without EBC are considered as ground truth, which are affected in 0.27 mm/day. In general, the use of in situ or ERA5 data as meteorological inputs have low sensitivity on S-SEBI model accuracy to estimate the daily ET over the Brazilian Pampa. Moreover, the results may suggest that when the satellite observations are compared against EC data without EBC, the fact of using reanalysis or in situ radiation data tends to have even lower influence in the daily ET. These findings may be helpful in situations where some of the in situ energy balance components are missing or are not trustable enough to perform a proper EBC (Table 2).
Transpiration, soil evaporation, and interception evaporation components are frequently embedded within remote sensing-based ET models. However, comparisons have demonstrated that although the total evaporative flux from different models agrees relatively well, the different components diverge substantially (Talsma et al. 2018). Even with an instantaneous G and consequently LE less accurate for S-SEBI in situ in comparison to S-SEBI ERA5 , the use of Λ together with the method used for extrapolating from instantaneous to daily (Sobrino et al. 2007), clearly contributes to produce more reliable daily ET estimates. Upscaling methods of instantaneous ET to daily ET usually have their own bias (Singh and Senay 2015) and variation depending upon seasonality and cloud conditions. The temporal scaling is usually a weakness of remotely sensed data (Zahira et al. 2009); nevertheless, using radiation data (either from reanalysis or in situ) for this transformation increases the accuracy of the estimates, especially because the global radiation behavior throughout the day is considered. Figure 7 shows the seasonal variation of the daily ET between both satellite-based measurements (S-SEBI ERA5 and S-SEBI in situ ) and field EC data with and without EBC. The daily ET followed a seasonal behavior with the highest values in the summer, period of fast growing for vegetation, and lowest in winter due to the cold weather. Over the Pampa, the ET is marked by strong seasonality, because it is highly dominated by precipitation rates, which directly affect the partitioning of the turbulent flux (Rubert et al. 2018). Overall, estimates from S-SEBI in situ fitted better in comparison to field data with and without EBC than the ones from S-SEBI ERA5 . ET from S-SEBI ERA5 was underestimated in most cases, but generally in warmer periods, it tended to have better accordance, especially when in situ data with EBC are considered as ground truth. Commonly in the Pampa biome, during cold periods, the small growth of the natural grasslands is very difficult to detect by satellite imagery because all of the biomass production is being consumed by cattle in a situation named overgrazing, which occurs when the consumption of biomass is larger than the aboveground net primary production of the vegetation utilized to feed the animals (Scottá and Fonseca 2015).

Seasonal and spatial pattern of ET
Differences between curves with and without EBC evidenced that measurements with EBC can better capture the variations of the ET seasonal pattern in the natural grasslands of the Pampa biome. Most variations found between ET in situ with and without EBC are seen in late winter and spring seasons, which are noticeable in the DOYs of 202, 266, 278, 282, and 331 (Fig. 7), demonstrating that  the imbalance issues in EC measurements without EBC can become more pronounced in these periods. Rubert et al. (2018) reported that the overestimation of ET from EC data normally is noticeable in the hydrological regime where water availability is the limiting factor. The authors mentioned that although the relationship between the ET and soil water content depends on soil type, vegetation type, and vegetation adaptation to dryness, the role of soil water content near the surface is expressive. According to Fontana et al. (2018), the general climate condition in the Pampa is considered uniform in terms of temporal distribution of rainfall throughout the year. However, events such as La Niña cause high vegetation water stress. Results indicated that the daily ET from remote sensing observations produced better agreement with in situ EC data in late winter and spring seasons. On the other hand, late summer and autumn exhibited worse estimates. Figure 8 investigates the spatial behavior of albedo, LST, Λ, and daily ET for two different DOYs. The DOY 58 represents a day with less accurate daily ET estimates in late summer, whereas DOY 278 has better daily ET results and is a spring day. Although DOY 278 exhibited a very accurate daily ET estimate, the S-SEBI model underestimated the daily ET for both cases (Fig. 8). Differences of 2.26 and 0.2 mm/ day against in situ measurements with EBC, and 3.02 and 0.52 mm/day in comparison to the in situ measurements without EBC, were found in the DOYs 58 and 278, respectively. LST and albedo play key roles in the energy balance. Typically, the lower the albedo is, less solar radiation ends up being reflected back into space, and more is absorbed by the surface, thereby increasing the LST (Findell et al. 2007). Thus, the ET rate increases with the increase of water content in soil due to this lower albedo. It is noticeable that in spring when the albedo starts to exhibit higher values the opposite is observed. Once the surface reflects more energy strongly reducing the LST, the albedo-LST relationship used in the S-SEBI can be compromised, directly affecting T H and T LE definition, and consequently Λ (Eq. 6). Thus, the model better represents the pattern of the relationship in spring where the higher data accuracy is seen. Despite the general behavior and accuracy of the S-SEBI in the Pampa biome is not that different from other grassland types, the climate condition commanded by the strong seasonality might prevent the suitable operation of the model in specific DOYs (such as 58), since it directly influences the relationship that is the basis of the model. As the increasing water scarcity allies to the ever-growing population has created urgency for integrated and sustainable water resources management (Haghverdi et al. 2021), this understanding is useful to establish conservation strategies over the Pampa in addition to monitoring the climate variability and its impact on the water resources availability.

Conclusions
The Pampa biome carries a rich biodiversity and is considered priority for conservation. In this context, accurate information on the ET dynamics is essential to support its proper monitoring. As different remote sensing-based models for ET estimation have different input data requirements, it is fundamental to understand its strength and limitations in addition to comprehend the uncertainties of the components and their impacts on final ET. Because a detailed evaluation of the S-SEBI in retrieving surface fluxes over the natural grasslands of the Pampa has not been properly addressed, we propose to investigate an operation methodology based on the model to estimate the water and energy fluxes.
The use of the S-SEBI model in conjunction with ERA5 (S-SEBI ERA5 ) radiation data allowed to obtain ET over natural grasslands with RMSE and MAE of 1.38 mm/ day and 1.12 mm/day, respectively. In contrast, by using in situ radiation data as input (S-SEBI in situ ), it produced RMSE and MAE of 1.06 mm/day and 0.79 mm/day, respectively. Thus, the RMSE is impacted in 0.32 mm/ day when reanalysis data are used to replace in situ radiation data, indicating that meteorological inputs have low sensitivity on daily ET estimates by the S-SEBI model. In contrast, the instantaneous components are more affected. The impact in the daily ET is lower when in situ data without EBC are considered as ground truth, despite they are less correlated with the remote sensing-based estimates. It can be helpful particularly when in situ energy balance components are missing and the EBC cannot be properly performed.
Our findings demonstrated that the strong seasonality of the region influences the evaporative fraction, which is the basis of the S-SEBI. This behavior is more evident in the transition between late summer and autumn and may compromise the performance of the model. The methodology proposed can guide the generation of accurate long-term records of ET, since it includes remote sensing products that can be easily obtained for different freely available satellites, such as the Landsat series that provides a continuously acquired collection of data. Moreover, these insights are useful to guide future development and application of remote sensing-based models over grasslands. Therefore, the operational method provides the opportunity to capture ET trends through the time over the natural grasslands of the Brazilian Pampa. Once selecting a model from the wide range available to retrieve ET under a given set of circumstances is challenging for users, in further work, we intend to compare the S-SEBI performance with other SEB models.