Methodology of spatialization of energy balance components using MODIS and ERA Interim data1

This work aims at proposing a methodology to spatialize the Energy Balance (EB) components from the combined use of MODIS sensor products and ERA-Interim reanalysis data, and to evaluate the performance of the OSEB and METRIC models for estimating EB components in the predominantly humid subtropical climate of Rio Grande do Sul, Brazil. The research was conducted for 3 consecutive years from 2009 to 2011. The temporal variation patterns of the EB components were analyzed by comparing the values estimated by the models against the reference measurements of a micrometeorological tower (Sulflux network). For the spatial variability analysis, reference evapotranspiration (ETo) was calculated at a local scale (data from Weather stations of the National Institute of Meteorology) and a regional scale (ERA-Interim reanalysis data) while the consistency between these estimates was also verified. Afterward, using only the pixels in which the sensible heat flux represented more than 50% of the radiation balance, correlation maps between ETo ERA-Interim and latent heat fluxes were constructed. The estimated BE components are consistent regarding both, the order of magnitude and variation throughout the year, when compared to the reference measurements of the micrometeorological tower. The spatial variability analyses also reveal satisfactory results for both models, however, the METRIC model results showed greater coherence since it provided a greater number of dates with consistent results, thus being indicated as the most suitable model for building the time series.


INTRODUCTION
Using orbital images in energy balance (EB) studies allows representing almost exclusively its components on a regional scale, continuously and sequentially, over time. This characteristic is especially useful in areas with intense agricultural activity, such as the state of Rio Grande do Sul (RS), which produces approximately 14% of the Brazilian grain production, especially soybeans, rice, and corn in the summer and wheat in the winter (JUNGES et al., 2016;KLERING et al., 2016;SANTOS et al., 2014), having water deficit as its main limitation to obtaining high yields (SENTELHAS et al., 2015).
Currently, several existing methods that present an adequate basis for estimating EB using images, are being applied in different parts of the world (ALLEN; TASUMI; TREZZA, 2007;BASTIAANSSEN, 2000;GAN et al., 2019;KUSTAS et al., 2016;TIMMERMANS;KUSTAS;ANDREU, 2015). Among them, the OSEB ('One-Source Energy Balance') and METRIC ('Mapping Evapotranspiration with Internalized Calibration') are highlighted. However, regional peculiarities require determining the accuracy of the models to define the one to be adopted. Validation is usually done using data measured at experimental sites equipped with 'Eddy Covariance' or 'Bowen Ratio' systems (DICKEN;COHEN;TANNY, 2013;KUSTAS et al., 2016). In RS, the SulFlux network provides these measurements for different locations in the state (http://www.sulflux.ufsm.br/sulflux/?page_id=26) (MOREIRA et al., 2015), but these are scarce and sparse for ensuring an adequate validation of estimates on a regional scale.
Another important aspect to be considered is that almost all models need data measured on the surface, generally in Weather stations. However, due to the limitations on spatial detailing, the ERA-Interim reanalysis data appears as an alternative. This database has been used in conjunction with remote sensing data in several applications, such as soil moisture, atmospheric correction, and in studies related to monitoring of crops and pastures (QIU et al., 2016;RÖTZER;MONTZKA;VEREECKEN, 2015;SÜTTERLIN et al., 2015). The quality of this database has been certified against the measured data (DEE et al., 2011), including in RS (MOREIRA et al., 2017, but its use for estimating the EB components has not been tested so far. In this context, the objective of this work is to propose a methodology for the spatialization of the Energy Balance components from the combined use of MODIS products ('Moderate Resolution Imaging Spectroradiometer') and ERA-Interim reanalysis database while also evaluating the performance of the estimates generated by the OSEB and METRIC models in the humid subtropical climate of Brazil.

MATERIAL AND METHODS
The study area comprises more than 280 thousand km 2 covering the entire state of Rio Grande do Sul, southern Brazil (Figure 1). The region is an important grain producer, especially soybeans, rice, and corn in the spring-summer and wheat in autumn-winter. This work analyzed 138 dates in three consecutive years between 2009 and 2011.
The used data originated from four different sources: a) MODIS Products: the Earth Surface Temperature -MOD11A2, Vegetation Index -MOD13A2, Leaf Area Index MOD15A2, and Albedo -MCD43B3 products were used. The spatial resolution of 1,000 m and temporal resolutions consisting of 16 day-compositions for MOD13A2, and 8 days for MCD4343B3, MOD15A2, and MOD11A2. To cover the study area, data obtained from the LP DAAC ('Land Processes -Distributed Active Archive Center'), https://lpdaac.usgs.gov/ website were used to mosaic the h13v11 and h13v12 tiles. b) ERA-Interim reanalysis data: global solar radiation, air temperature, dew point temperature, and wind speed were used from this database. Records corresponding to synoptic times were acquired with 0.25 x 0.25 spatial resolution and available in a global scope matrix. The data were converted from the original format (NetCDFeach file contains the monthly records of a variable) into matrices for use in an image processing environment and, afterward, the matrix clipping was done to cover the Rio Grande do Sul area. For each variable, the data recorded at the moment of the satellite passage were extracted, the daily average value and the 8-day average, equivalent to the period of the MODIS images, were calculated and the interpolation was performed for a spatial resolution compatible with the MODIS matrix. In the interpolation process, due to the large difference of resolution, the original matrix of the reanalysis data was converted into a point cloud considering the central coordinate of each pixel, and the Kriging interpolation model was applied to this point cloud (ISAAKS; SRIVASTAVA, 1989). c) Reference measurements of the micrometeorological tower in Cruz Alta: the used reference data were measured in an experimental plot cultivated with soy in the summer and wheat in the winter, in the municipality of Cruz Alta, RS, at 28.600859° latitude and 53.672153°l ongitude, and 432 m altitude ( Figure 1). This site is part of the SulFlux Network (http://www.sulflux.ufsm.br/ Methodology of spatialization of Energy Balance components using MODIS and ERA Interim data Performed at a height of 3 m from the soil, the radiation balance (Rn) was measured using the Kipp & Zonen -NR LI TE sensor, the soil heat flux (G) with the Hukseflux -HFP01SC-L sensor, whereas the turbulent sensible (H) and latent (LE) heat fluxes were estimated using the 'Eddy Covariance' method. For this purpose, high-frequency measurements (10 Hz) of a 3D sonic anemometer, CSAT3 (Campbell Scientific Inc.), and an infrared gas analyzer (LI-7500, LI-COR, Inc.) were used. The EB components were computed in 30 min averages following Moreira et al. (2015), and the LE and H data were corrected using the forced closure technique of Bowen ratio (DICKEN; COHEN; TANNY, 2013). d) Meteorological data measured on the surface: the data used were obtained from INMET Weather stations located in different sites in the state ( Figure 1) and referred to air temperature, dew point temperature, global solar radiation, and wind speed.
Using the MODIS and ERA-Interim data, the EB components were estimated by the OSEB and METRIC models, and in both models, the radiation balance is obtained from the Eq. 1.
In Eq. 1, the α, Ts and ɛ s parameters were obtained from the MODIS products and the others, from the ERA-Interim reanalysis data.
The OSEB and METRIC models differ mainly regarding how the sensible heat flux is obtained. In the OSEB model, H is obtained directly from the surface (MOD11A2) and air (ERA-Interim) temperature gradient. In these models the aerodynamic resistance considers the reference surface characteristics as proposed by Allen et al. (2006) and estimated from wind speed data (ERA-Interim).
The METRIC model, on the other hand, uses hot and cold anchor pixels for determining H, based on a self-calibration iterative process of the aerodynamic resistance and vertical temperature differential in the first meters of the atmosphere. In the iterative process, the MOD11A2 surface temperature data and leaf area index (LAI) data from the MOD15A2 product were used to estimate the roughness parameter, besides ERA-Interim air temperature, and wind speed data. The anchor pixels consider the image water limits. The cold pixel is defined as LE ≈ 1.05 * ETo (reference evapotranspiration) (ALLEN et al., 2006), and H = Rn -G -LE. The hot pixel considers the occurrence of residual evaporation from the uncovered soil, LE = E (evaporation in the soil obtained from the soil water balance proposed by Allen et al. (2006)), and H = Rn -G -LE.
Another relevant aspect is that H estimation in the METRIC model requires selecting the extreme temperature pixels in each image (LONG; SINGH, 2013). In the present work, the hot and cold anchor pixels were selected objectively following three criteria. First, the reference values were based on grouping pixels, never on a single pixel. Second, the selection was based on statistical criterio, considering cold pixels those pixels that have temperature valuescorresponding to the accumulated frequency of 2% while the hot pixels correspond to the accumulated frequency above 98%. Third, pixels close to those with the presence of clouds were eliminated from the selection process (a 4-pixel buffer was applied to the cloud mask of the Ts product).
Considering the ecoclimatic conditions of the region and the resolution of the MODIS images, whose 1 km pixel hampers the occurrence for pixels completely without vegetation; a 50% vegetation cover was assumed for determining the residual LE of the hot pixel from the water balance for the METRIC model.
The instantaneous values of the energy components, obtained from the image processing, were converted into daily values following the methodology proposed by Rivas and Carmona (2013), which consists of obtaining linear regression between radiation balance measurements, carried out on the ground, at a time of satellite overpass and the daily values, for later application of the equation in the conversion process of the image data.
After applying the models for estimating the EB components and converting them to daily values, the pixels with inconsistent results were eliminated, remaining only those with Rn and LE values greater than 0 and less than 400 W m -2 , as well as H greater than -100 W m -2 and less than 250 W m -2 . These thresholds were defined from the analysis of variability and occurrence of outliers in the data recorded in the micrometeorological tower of Cruz Alta.
To analyze the results regarding the temporal distribution patterns, graphs of the partitioning LE and H components were constructed throughout the year and the RMSE and MBE errors were calculated between the values estimated by the OSEB and METRIC models and the reference measures (micrometeorological tower).
For evaluating the spatial adequacy of the LE estimates, ETo reference evapotranspiration maps were generated (ALLEN et al., 2006), obtained from the ERA-Interim (ETo ERA) data for the entire state of Rio Grande do Sul. Subsequently, the ETo ERA and LE data matrices obtained from the OSEB and METRIC models were compared. To minimize the differences inherent to these two distinct physical variables, the comparison of the results followed a criterion that consisted of using only pixels with LE values equal to or greater than 50% of Rn. The pixel by pixel correlation coefficient was calculated, and a correlation map between the LE and ETo ERA variables was obtained. In this process, a map of the number of days that met the LE condition greater than 50% Rn was also generated. The LE and ETo ERA data extracted from the matrices (3 x 3 window centered on the stations' coordinates) were used for plotting the dispersion graphs while the correlation coefficients and RMSE and MBE errors were also calculated.

RESULTS AND DISCUSSION
Over the three year-period analyzed, the partition of the EB components determined using the reference measurements (micrometeorological tower) and the local estimates obtained by the OSEB and METRIC models were consistent. The mean values and standard deviations obtained for Rn and G, in both models, were similar (Table  1), which was already expected, since these components are, generally, easier to obtain (TIMMERMANS et al., 2007). On the other hand, the LE and H estimation is considered the most complex (BOEGH; SOEGAARD; THOMSEN, 2002;CAMMALLERI et al., 2012). The LE was obtained as the residual term of the EB equation and had an error similar and less than 50 W m -2 in both models. Furthermore, the greatest discrepancy between the results of the two models was found for the H estimation, with higher accuracy (lower RMSE) observed for the OSEB model.

Radiation balance (Rn); latent evaporation (LE) heat flux; sensible heat flux to the air (H) and ground (G); Standard Deviation (SD); Root mean square error (RMSE); Mean deviation error (MBE)
Some of the uncertainties observed in the models may result from the fact that these models were developed and adjusted for dry or semi-arid climate regions (KUSTAS et al., 2016;TIMMERMANS et al., 2007). It is also noteworthy that, although the tower measurements are being used as the reference data for evaluating the models (measurement location of all EB components), differences are expected. Both model estimates are obtained from data representing the average condition of the surface in a 1 km 2 area based on the best atmospheric conditions over 8 consecutive days. On the other hand, the data from the tower is a point measurement representing the average condition of these 8 days. Despite this, Rn temporal variability ( Figure 2) and order of magnitude (Table 1) calculated for the micrometeorological tower measurements and the estimates of both models were remarkably close. The variability observed over the year is typical of regions with a subtropical climate, where Rn is mainly governed by the incident solar radiation, which causes great variability in the energy supply between the summer and winter periods (MATZENAUER; RADIN; ALMEIDA, 2011).
Additionally, similar temporal patterns were observed for the LE and H fluxes determined by measurements recorded in the micrometeorological tower The greatest discrepancy between data was observed during the winter growing period (mainly at the beginning of the cycle) and with special emphasis on the results of the METRIC model. However, it should be noted that the availability of data for use in the models reduces significantly in the winter. Of the total 136 dates over the 3 years, considering the periodicity of 8 days, OSEB and METRIC presented consistent results on 78 and 107 dates, respectively. This data loss is mainly due to image data inconsistencies, such as cloud cover in the micrometeorological tower area. For the OSEB model, on the other hand, the losses also result from the high sensitivity to air and surface temperature data, which provide inconsistent H values (CAMMALLERI et al., 2012). This information is important and should be considered when choosing a method for estimating EB and building a time series as well as studies on variability and trends.
In the next step, LE was compared to ETo to assess the consistency of the estimates under varying weather or land cover conditions in RS. LE (amount of energy consumed in the evapotranspiration process) and ETo (evaporative demand from the atmosphere) are distinct variables and only present approximate magnitudes under full soil water supply conditions and, therefore, the dataset used in this evaluation includes only the days when the models had LE/Rn greater than or equal to 50%, which tends to minimize, The dispersions (Figure 3) were similar for both models, with correlation coefficients varying between 0.69 and 0.77 (OSEB) and between 0.65 and 0.82 (METRIC). But, in general, a higher correlation coefficient, and lower RMSE and MBE, were obtained for the METRIC model compared to OSEB, except for the Bom Jesus weather station. Also, most of the dispersion points are observed below the 1: 1 line (Figure 3), which is expected, since the evaporative demand of the atmosphere (ETo) tends to be greater than the energy effectively spent in LE.
The maps of the correlation coefficients between the ETo ERA and LE estimated by the OSEB and METRIC models are shown in Figure 4, which, calculated pixel by pixel, demonstrated that the variables for both models are highly associated in practically the entire state. Additionally, in the maps, both models show similar patterns regarding the number of days in which LE / Rn was greater than 50%. Higher values were observed in areas with forest cover and more rugged relief (from 70 to 90 days), while intermediate values (from 40 to 60 days) occurred in Campos de Cima da Serra, Central Depression, and the southern part of the ecoclimatic region of Campanha (border with Uruguay). The greatest differences occurred especially in the Planalto Médio region and those bordering Argentina, where LE/Rn greater than 50% was recorded in 40 and 60 days by the OSEB model, as well as 20 days and 30 days by the METRIC model, respectively, for these two regions (Figure 4c and 4d).
The combined analysis of the maps shows coherence between the correlation coefficients and the number of days for both models. The greatest number of days occurred in the areas with predominant forest cover, where the remnants of the Atlantic forest allow the system to control better the temperature due to the ease of converting energy into LE and consequent reduction in H (DANELICHEN et al., 2014). On the other hand, the LE variability in this type of vegetation is lower over time (compared to other vegetation types), which, combined with the greater availability of days, provides lower correlation coefficients between ETo and LE. In the regions bordering Argentina and the Planalto Médio, the METRIC model provided more consistent results. Close to the border, the predominant countryside vegetation and climate, with higher temperatures and lower relative humidity (MATZENAUER; RADIN; ALMEIDA, 2011), increases the atmospheric evaporative demand which, combined with the shallower soils, increases the water deficit, thus reducing the number of days when LE/Rn is greater than 50%.
In the Planalto Médio region, on the other hand, the more intense agricultural activity allows for greater temporal variability of LE associated with the cycle of annual crops, which combined with the lower availability of data, yields a higher correlation coefficient ( Figure  4b).
In general, both models presented satisfactory results for estimating the EB components. The statistical analysis of the results (Table 01 and Figure  02) indicates that the OSEB and METRIC models alternate between themselves, yielding both lower RMS errors, a higher correlation between LE data obtained from Images, and ETo obtained from data from stations at times. However, generally speaking, the OSEB model had smaller errors, although not significantly different between the two models. On the other hand, the spatial distribution of the METRIC model was more consistent with the state ecoclimatic conditions. Therefore, choosing the most appropriate method for representing EB components in Rio Grande do Sul requires additional information other than the statistical results.
The METRIC model is based on semi-empirical concepts to obtain characteristics of the first meters of the atmosphere, determining H in an interactive process using the anchor pixels that represent the water limits of the image (ALLEN; TASUMI; TREZZA, 2007;BASTIAANSSEN, 2000). Despite this, the definition of residual LE in the pixel that represents the image dry limit introduces greater rigor and coherence to the model for application in humid subtropical climate conditions. Another point to note is that the METRIC model is less dependent on the quality of variables such as air temperature and/or surface temperature than the OSEB model (CAMMALLERI et al., 2012).
Additionally, another necessary consideration is the great diversity of databases and images needed to estimate the EB components in both methods tested, but especially in the METRIC model. This volume of data is inserted in the current scenario of geoprocessing, the "BigData" concept, driven by the increasing processing capacity and programming languages applied to data processing and analysis. In this scenario, and considering that the objective of the proposed methodology is to build a representative time series, even though the METRIC model requires a larger amount of input data, it can be appointed as the most suitable for building the EB time series. Figure 3 -Scatter plots between the LE latent heat flux obtained from the OSEB and METRIC models and the ETo reference evapotranspiration. LE data comprise the estimates in which LE/Rn greater than or equal to 50% Methodology of spatialization of Energy Balance components using MODIS and ERA Interim data Figure 4 -Spatial distribution and histograms of the correlation coefficient between ETo ERA and LE data (a and b for the OSEB; c and d for the METRIC) and the number of days that met the condition, LE/Rn greater than 50% (e and f for the OSEB; g and h for the METRIC). Histograms show the accumulated frequencies of 2% and 98%, and average values in the vertical lines, in blue, red, and black, respectively CONCLUSIONS 1. The OSEB and METRIC models present consistent results for the energy balance components regarding their order of magnitude, as well as their temporal and spatial distribution patterns; 2. The OSEB model provides estimates with lower deviations and errors, but the differences between the errors of both models were not significant; 3. The methodology for estimating the energy balance components based on the combined use of the MODIS and ERA-Interim products was consistent while allowing to conclude that the METRIC model is the most suitable for building a time series due to greater availability of days with valid data and more coherent spatial distribution. Using a larger number of input variables reduces the dependence of the model accuracy on certain variables, as observed in the OSEB model, whose result is highly correlated with the difference between the surface and air temperatures, becoming difficult to accurately determine on a regional scale.