Evaluation of OSEB and SEBAL models for energy balance of a crop area in a humid subtropical climate

Bragantia, Campinas, v. 77, n. 4, p.609-621, 2018 ABSTRACT: This study evaluated the adequacy of the One-Source Energy Balance (OSEB) and Surface Energy Balance Algorithm for Land (SEBAL) models to estimate evapotranspiration in grain growing areas with humid subtropical climate in Rio Grande do Sul, southern Brazil. The dataset was obtained from a micrometeorological station (Eddy Covariance) and MODIS (Moderate Resolution Imaging Spectroradiometer) products during 84 observation days between 2009 and 2011. The OSEB and SEBAL models were used to estimate the partition of net radiation (Rn) into the components latent heat flux (LE), sensible heat flux (H), and ground heat flux (G) estimated from the MODIS images while the experimental data measured in situ AGROMETEOROLOGY Article

ABSTRACT: This study evaluated the adequacy of the One-Source Energy Balance (OSEB) and Surface Energy Balance Algorithm for Land (SEBAL) models to estimate evapotranspiration in grain growing areas with humid subtropical climate in Rio Grande do Sul, southern

INTRODUCTION
The correct understanding of the partitioning of energy and mass fluxes on the Earth surface is extremely important for climate, hydrological and agrometeorological studies.Nevertheless, scalar measurements of energy fluxes are generally restricted to specific points such as micrometeorological sites since measurements are costly and reliant on expensive instruments and specialized staff.On a regional scale, Remote Sensing (RS) is an excellent tool to monitor the spatial distribution and temporal evolution of natural vegetation or agricultural crops.RS images acquired at different wavelengths allow determining the surface physical properties necessary to measure the energy and mass flux between the surface and atmosphere on a regional scale (Boegh et al. 2002;Kustas et al. 2004;Timmermans et al. 2007).
Most studies estimating the energy balance (EB) components are based on one-dimensional flux model describing the exchange mechanisms of radiation and heat fluxes between the surface and atmosphere while observing the energy conservation principle (Brutsaert 1984).The EB can be defined by how the net radiation (Rn) of the surface is divided into latent heat flux (LE), sensible heat flux to the air (H), and ground (G) (Friedl 2002;Timmermans et al. 2007), but most studies neglect the horizontal heat advection and heat storage at the canopy layer.
Most models for estimation of EB components by RS are based on the combined use of meteorological and image data.The Rn and G components are easily estimated from the RS data.Rn can be spatialized from albedo and thermal image while G is estimated as a portion of Rn, proportional to the vegetation index.Other components such as the turbulent fluxes, LE and H, are more complex to estimate.In general, LE is obtained as a residual component from the EB equation whereas H is estimated from the thermal images.The models differ greatly in the methodology to estimate H.
The Surface Energy Balance Algorithm for Land (SEBAL) (Bastiaanssen 2000) model has been applied using images from different sensors and locations across the globe.In Rio Grande do Sul, Brazil, the SEBAL model using high-resolution images has already been applied (Santos et al. 2010;Monteiro et al. 2014).However, because it has been parameterized for semi-arid climate, its suitability for different climatic conditions must still be verified, since the method presupposes the existence of extreme water conditions.The selection of points with extreme moisture conditions, the hot and cold pixels, is a requirement for determining H.However, this water condition requirement may not be fulfilled in many cases, especially in humid climate.Furthermore, defining these extremes is even more difficult when using images of moderate or low spatial resolution.
This study aimed at verifying the adequacy of the One-Source Energy Balance (OSEB) and SEBAL models to estimate the latent heat flux in grain producing areas in the humid subtropical climate of Rio Grande do Sul, using low-resolution MODIS (Moderate Resolution Imaging Spectroradiometer) images.The performance of the mathematical models was evaluated by comparing the results to the experimental data of the energy balance components obtained from the experimental area in Cruz Alta, Rio Grande do Sul.

MATERIAL AND METHODS
The study site is located in northwestern Rio Grande do Sul, Brazil (Fig. 1).This grain producing area is greatly important to the state and country.The climate in the region is humid subtropical, classified as Cfa (according to Köppen) with highly variable rainfall and frequent droughts, which are responsible for the yield differences observed in the state.
The study period lasted three years, from 2009 to 2011.The reference data from the micrometeorological station in Cruz Alta were used to check the performance of the OSEB and SEBAL mathematical models to estimate the EB component estimates.
Over the period, 84 days were selected for analysis.The selection criteria were completely clear sky day over the area throughout the day (verified from the global solar radiation data) and the availability of MODIS products with adequate quality.
The analyzed dates were distributed over the period of three years, which allowed analyzing the models in various weather conditions.Of the analyzed dates/images, 20 corresponded to summer crops, 32 winter crops, and 32 partial vegetation cover.The NDVI profile pattern \and analyzed dates (Fig. 1c) were distributed to represent the variability of results according to time of year/season, crop type, and degree of surface coverage.
Of the EB components, Rn and G are easily estimated from RS images (Kustas et al. 2004;Sánchez et al. 2008).On the other hand, the estimation of the turbulent fluxes, LE and H, is more complex.In most EB models, LE is obtained as a residual term in the EB equation.
The main difference between models to estimate the EB components from remote sensing data lies in the way the sensible heat flux (H) is determined.Among them, three approaches are highlighted.In fact, the OSEB models are based directly on the radiometric temperature difference between the vegetated surface and air temperatures (Boegh et al. 2002;Friedl 2002, Kustas et al. 2004, Wang et al. 2006;Sánchez et al. 2008;Timmermans et al. 2007;Tang et al. 2013).Likewise, the TSEB (Two-Source Energy Balance) models treat differently the heat exchange between the atmosphere and vegetated areas, and between the atmosphere and areas with bare soil based on different equations to determine the components in soil and vegetated areas (Sánchez et al. 2008, Cammalleri et al. 2012, Tang et al. 2013).Finally, the third approach, the SEBAL (Surface Energy Balance Algorithm for Land) (Bastiaanssen 2000), the SSEBI (Simplified Surface Energy Balance Index) (Roerink et al. 2000, Mattar et al. 2014) and METRIC (Mapping Evapotranspiration with Internalized Calibration) (Allen et al. 2007) models could be considered as a subfamily of OSEB models, seeking to represent the spatial variations of the vegetation and atmosphere water availability by mapping the cold and hot image pixels.Thus characterizing, respectively, the maximum latent heat flux (LE = maximum and H = 0) and maximum sensible heat flux (LE = 0 and H = maximum).
The studied (OSEB and SEBAL) models are based on the principle of energy conservation, Brutsaert (1984) (Eq.1): The first term of the EB equation, net radiation (Rn), can be estimated from the satellite images using Eq.2: where: R G↓ is the global solar radiation (MJ·m -2 ·day -1 ); a is the surface albedo (dimensionless); ε s is the surface emissivity (dimensionless); ε a is the atmosphere emissivity (dimensionless); s is the Stefan-Boltzmann constant (4.9 × 10 -9 MJ·m -2 ·K -4 ·day -1 ); T a is the air temperature (K); and T s is the surface temperature (K).
The a, T s , and ε s parameters were determined from remote sensing thermal images while the other components from the experimental data obtained at the micrometeorological station (Boegh et al. 2002;Sobrino et al. 2005;Rivas and Caselles 2004).
For the OSEB and SEBAL models, the latent heat flux, LE, was estimated as the residual term of Eq. 1.
The main difference between models OSEB and SEBAL models is how the third term of the EB equation, the sensible heat flux to the atmosphere H, which represents the portion of absorbed radiation released as heat from the surface to the air, is determined.Both models estimate H from a temperature differential and aerodynamic resistance, according to Eq. 3: where: ρ is the air density (1.15 Kg·m -3 ); cp (1004 J·Kg -1 ·K -1 ) is the specific heat of humid air at constant pressure; dT is the differential temperature between two levels on the surface; and ra is the aerodynamic resistance to heat transport (s·m -1 ).
In the OSEB model, the differential temperature, dT, is the difference between surface (obtained from satellite image) and air temperatures, normally obtained from meteorological stations at 2 m above the surface (Boegh et al. 2002;Kustas et al. 2004;Wang et al. 2006), and ra is obtained from wind speed according to Allen et al. (1998).
For the SEBAL model, the differential temperature and aerodynamic residence are obtained from an iterative calculation procedure, which seeks atmospheric stability condition for each image pixel, based on anchor pixels that represent the extreme soil moisture of the image.The detailed methodology for obtaining these parameters is explained by Bastiaanseen (1995).
In both models, ground heat flux G is estimated as a fraction of Rn, inversely to the vegetation indices (Allen 1998;Boegh et al. 2002;Qi et al. 1994;Moran et al. 1989).In the OSEB model, G is estimated from Eq. 4 (Moran et al. 1989): The SEBAL model estimates G from NDVI, but it adds the albedo (α) and surface temperature (T s ) data using Eq. 5 (Bastiaanssen 2000): Ts/α)(0.0038α + 0.0074α)(1-0.98NDVI 4 )]Rn (5) The EB components were obtained from the MODIS products: Earth surface temperature -MOD11A2; vegetation index -MOD13A2; Albedo -MCD43B3; and leaf area index -MOD15A2.All used products had 1,000 m spatial resolution while temporal resolutions consisted of temporal compositions of 16 days for MOD13, 8 days for MCD43B3 and MOD15A2, and daily for MOD11A2.
These products are available as daily scenes covering a specific area of the globe.A mosaic of the H13V11 and H13V12 sinusoidal tiles, obtained from the Land Processes Distributed Active Archive Center website (https://lpdaac.usgs.gov),was prepared to cover the study site.
The EB components were obtained experimentally at 3 m high, using the data collected from the microteorological station, equipped with net radiation sensor (Rn), ground heat flux sensor (G), measuring in low frequency 1 Hz, 3D sonic anemometer and infrared gas analyzer measuring in high frequency (10 Hz).The turbulent fluxes of sensible (H) and latent heat (LE) were estimated by the Eddy Covariance method, obtaining 30 min averages.The EB sensor has been described in Table 1.The station have other sensors to measure atmospheric pressure, rainfall, shortwave incident radiation, photosynthetically active radiation, and soil temperature.
The Eddy Covariance system estimates all EB components independently, generally the (LE + H) is different from the available energy (Rn -G) and should be adjusted to close the EB.The method applied to close the EB is based on maintaining the Bowen Ratio while the available surface energy is re-partitioned (Tang et al. 2013).
The micrometeorological measurements were performed in a 40 × 60 m experimental plot in Cruz Alta, Rio Grande do Sul, located at -28.6036 latitude, -53.6736 longitude, and 432 m altitude.The plot was planted with soybean (summer crop) and wheat (winter crop), the main cultures of the area, representing the agricultural lands surrounding the station.
The analyses were separated into three periods: the summer and winter crops, and partial vegetation cover.The results are presented in sets with temporal NDVI vegetation index profiles to understand the dynamics of the crops and the partitioning of the energy components.

Energy balance in humid subtropical climate
The model results (average of 3 × 3 pixel window centered at the micrometeorological station) were compared with the reference measurements performed at the micrometeorological station.Those analyses were performed using dispersions for both component and model.The quality of the estimates was checked by the Root Mean Square Error (RMSE), Mean Bias Error (MBE, computed as observed -estimated), and the analysis of the EB component ratios in relation to Rn.

Partition of the EB components -experimental data
The pattern of EB components measured in the micrometeorological station is shown in Fig. 2.These plots show the mean data recorded every half hour during 20 random days in the summer crops (Fig. 2a), 32 in the winter crops (Fig. 2b), and 32 days over partial vegetation cover (Fig. 2c).The EB closure has been applied to this data.
The results for the energy balance components showed higher Rn values for the summer crops (Fig. 2a), lower and similar for the winter crops (Fig. 2b) and partial vegetation cover (Fig. 2c).The results for the energy balance components showed that Rn was higher for the summer crops (Fig. 2a), and lower and similar for the winter crops (Fig. 2b) and partial vegetation cover (Fig. 2c).The highest energy input varied between 796 W•m -2 and 580 W•m -2 , with mean of 698 W•m -2 in the summer, and between 744 W•m -2 and 302 W•m -2 , with mean of 474 W•m -2 in the winter; these measurements were recorded at 12:30 PM (local time) when available radiation is maximized.The observed values for Rn and LE, H and G fluxes are associated with the subtropical climate of the study area and varying coverage over time.Figure 1c shows that the summer crop cycles were shorter and had higher NDVI, which resulted from the higher energy availability, rapid vegetation development, and higher biomass accumulation.The opposite was observed in the winter crops.In the partial vegetation coverage period, the increasing albedo resulted from more exposed soil, senescent vegetation, which increased short-wave radiation loss and provided Rn value similar to winter crops (Fig. 2), with lower energy input.
At the time of the satellite overpass (Table 2), the Rn partition was similar to that observed throughout the day.This effect is called evaporative fraction self-preservation, providing, under some circumstances, the tool to extrapolate from instantaneous satellite retrievals to daily estimates of evapotranspiration (Cragoa and Brutsaert 1996;Cammalleri et al. 2014).The largest partition of energy is consumed by the latent heat flux, followed by sensible heat flux in the air and soil (Fig. 2).
As already mentioned, the evapotranspiration process consumes the most Rn due to the humid subtropical climate of the region.However, despite the damp climate, soil water availability is often insufficient to meet the evaporative demand of summer crops, and drought becomes responsible for yield variability hence the need to quantify it.For these same crops, the sensible heat flux, H, was the second component with the highest portion of available energy.
The highest H values were observed during the partial vegetation cover period while LE decreased (Fig. 2).This significant decrease of the LE component is possibly due to senescent vegetation accumulated on the ground since straw creates an insulating layer, partially interrupting the soil evaporation process.The soil heat flux, which is usually responsible for lower energy consumption, was lower than 8% for all crops.
Figure 3 shows how Rn was split between the LE and H fluxes in all analyzed dates (the graphs show H + LE).For most dates, LE is higher than H, consistent with the pattern observed in the mean values.However, for a few days, especially between the 97 and 145 days at the end of the summer crops of 2009 and 2010, H was higher than LE.Variations were observed only in the Rn partition proportions, a result of the variable weather conditions and soil water availability.Also, it is noteworthy that the LE and H values are close in the initial phase of winter crops.

Partition of the EB components -the OSEB and SEBAL models
The experimental energy balance components were used as a reference for understanding how Rn partition varies over time, and the factors influencing it.The objective of this study was to determine reliable methods, mathematical models, to estimate EB from satellite images to allow drawing up maps showing the spatiotemporal distribution of EB. Figure 4 shows the scatter plots comparing the estimated and experimental data of the EB components for the OSEB and SEBAL models.The Rn estimated by Eq. 2 behaved similarly for both models and very close to the experimental values, with the dispersion points approaching a straight line 1:1 and RMSE equal or below 50 W•m -2 (Table 3) for all three vegetation covers.Overall, Rn is the easiest parameter to be estimated and has the lowest error in EB studies from the images.Tang et al. (2013) reported similar magnitude errors (31 W• m -2 ) for the OSEB model in areas with maize (summer) and wheat (winter) crops in northwestern China.Also, Timmermans et al. (2007) reported errors close to 44 W•m -2 for the SEBAL model in pasture areas in the semi-arid and sub-humid climates of Arizona and Oklahoma, respectively, in the United States.
The estimated and experimental H values showed the highest dispersion for both models.The OSEB model, based on the differential temperature between the air and surface, and aerodynamic resistance (Allen et al. 1998;Kustas et al. 2004;Tang et al. 2013), provided low H values and close to a straight line 1:1 for the summer and winter crops, with RMSE less than 51 W•m -2 (Table 3).However, for partial vegetation cover, the estimated H had a higher dispersion and RMSE.
In the winter crops, H was estimated as 20% and 40% of Rn for the OSEB and SEBAL models, respectively.For the partial cover period, experimental H measured in the station represented only 34% of Rn and 65% of Rn in the SEBAL model, exceeding the LE values.Likewise, Timmermans et al. (2007) compared the SEBAL and TSEB models and reported that the greatest deviations in the SEBAL model occurred in areas with bare soil.
Both models estimated low G values, similar to the experimental measurements obtained in the meteorological station in Cruz Alta.
LE is obtained as a residual term in the EB Eq. 1 in both models.Therefore, its estimate is directly linked to the performance of the H estimate, which is responsible for the second largest partition of energy consumption at the surface level, and whenever the H value is underestimated, the LE value is overestimated.This was observed in the OSEB model in the summer and winter images with MBE about -40 W•m -2 .The SEBAL model overestimated H in The H and LE values estimated by the OSEB (Fig. 5) and SEBAL models (Fig. 6) behaved similar to the experimental data (Fig. 3).The OSEB (Fig. 5) results show that, generally, LE is higher than H; however, this ratio was inverted especially at the end of the summer crops.This inversion, the H higher than LE, has also been observed during a few days in the early winter crop, but it was not observed in the micrometeorological station.At the station, we observed that LE was greater than H right after the winter crops had been implemented.This difference probably occurred because the scheduled plantings of winter crops in the region and the station footprint are not representative of the MODIS pixel area.
Unlikely, the SEBAL model estimated H higher than LE at various times over the period analyzed (Fig. 6).Monteiro et al. (2014) used Landsat images and reported overestimated H in four of the six images analyzed in Rio Grande do Sul, while H was higher than LE in one of the images.
In the SEBAL model, the LE and H distribution is highly dependent on the correct definition of cold (LE maximum and H zero) and hot (LE zero and H maximum) pixels.Therefore, three hypotheses have been proposed to justify the uncertainties of the SEBAL model when determining the LE and H proportions during partial vegetation cover periods and beginning of winter crops.The first hypothesis is that the humid climate of the study area hinders the occurrence of dry areas to the point that evapotranspiration does not occur (LE zero), especially when the images are of a low energy availability period with lower evaporative atmospheric demand.The second one is that the low spatial variability of temperature in the winter does not allow the correct determination of the extreme water conditions proposed by the model.Finally, the third hypothesis is that the 1,000 m spatial resolution of the MODIS sensor for the thermal bands homogenizes the temperature patterns, making it difficult to locate the small areas with appropriate water conditions to determine the hot and cold pixels.
The spatial resolution hypothesis has also been addressed by Kustas et al. (2004).These authors analyzed how image spatial resolution affects the different EB components and found that lower spatial resolution has more influence in the H estimation.The study showed that a particular area, with 120 m pixel, has coefficient of variation 0.42, while the same area imaged with 960 m pixel (approaching the 1,000 m of this study) lowers the coefficient of variation to only 0.26, increasing the average H value of the area.Similarly, Roerink et al. (2000) addressed the deficiency of using NOAA or MODIS images with large pixels, which reflect a combination of different surfaces, hindering the occurrence of pixels with extreme water conditions.Figure 7 shows the temperature histograms of MODIS images for two of the studied days, 129/2009 and 290/2011, which are examples of the effects related to the homogenization of surface temperature patterns observed in one 1-km pixel and low winter temperature range that hinder the correct partition between LE and H in the SEBAL model.In both cases, the temperature variation between the cold (vertical blue line) and hot (vertical red line) pixels of the image was approximately 10 °C.The temperature of the reference pixels (vertical black line) is close to the maximum temperature of the image.Given that the LE and H proportions are distributed between the hot and cold limits of the image, the proximity of the station portion with the maximum temperature causes the SEBAL model to overestimate H.
The authors suggest further studies to test whether models such as the TSEB, which are independent of determining the hot and cold pixels, could provide good estimation of the EB parameters for the humid climate of the region.
Brazil.The dataset was obtained from a micrometeorological station (Eddy Covariance) and MODIS (Moderate Resolution Imaging Spectroradiometer) products during 84 observation days between 2009 and 2011.The OSEB and SEBAL models were used to estimate the partition of net radiation (Rn) into the components latent heat flux (LE), sensible heat flux (H), and ground heat flux (G) estimated from the MODIS images while the experimental data measured in situ AGROMETEOROLOGY -Article

Figure 1 .
Figure 1.Spatial and temporal scales of the study; (a) site location in Rio Grande do Sul, Brazil; (b) micrometeorological station; (c) dates of the images selected for analysis over the three-year period (gray lines) and NDVI in Cruz Alta experimental site.

CSAT3Figure 2 .
Figure 2. Experimental average daily cycle of the components of the energy balance in Cruz Alta, Rio Grande do Sul, Brazil; (a) summer crops; (b) winter crops and (c) partial vegetation coverage.

Figure 3 .
Figure 3. Partition of the Energy Balance into the LE (latent heat flux) and H (sensible heat flux) components, measured with the micrometeorological station.Sum of H (red), LE (green) and NDVI (Normalized Difference Vegetation Index) profiles during the years: (a) 2009; (b) 2010; (c) 2011, in Cruz Alta, Rio Grande do Sul, Brazil.

Figure 4 .
Figure 4. Scatter plots of experimental versus mathematical simulation of the energy balance components for the different vegetal covers and the 84-day period.LE = latent heat flux; H = sensible heat flux; Rn = net radiation; and G = ground heat flux.Simulations using: (a), (c) and (e), the OSEB model; (b), (d) and (f), the SEBAL model.The experimental values correspond to the measurements conducted in the meteorological station in Cruz Alta, Rio Grande do Sul, Brazil.

Figure 5 .Figure 6 .Figure 7 .
Figure 5. Partition of the Energy Balance into the LE (latent heat flux) and H (sensible heat flux) components, estimated using the OSEB.Sum of H (red), LE (green) and NDVI (Normalized Difference Vegetation Index) profiles during the years: (a) 2009, (b) 2010; (c) 2011, in the experimental site in Cruz Alta, Rio Grande do Sul, Brazil.
used by the SEBAL model (LE = RN -G and H = 0) Hot pixel temperature used by the SEBAL model (LE = 0 and H = RN -G) Temperature of the pixels in the reference portion in Cruz Alta, Rio Grande do Sul, Brazil(a) (b)

Table 1 .
Specifications of sensor BE components installed at the micrometeorological station in Cruz Alta.

Table 2 .
(Tang et al. 2013rage values of the EB components for different vegetation covers obtained experimentally in Cruz Alta, Rio Grande do Sul, Brazil.The EB closure was achieved by maintaining the Bowen Ratio(Tang et al. 2013).

Table 3 .
Energy balance components estimated using the OSEB and SEBAL models, extracted from the images in the coordinates of the meteorological station in Cruz Alta, Rio Grande do Sul, Brazil.