A MODIS-Based Energy Balance to Estimate Evapotranspiration for Clear-Sky Days in Brazilian Tropical Savannas

Abstract

Evapotranspiration (ET) plays an important role in global climate dynamics and in primary production of terrestrial ecosystems; it represents the mass and energy transfer from the land to atmosphere. Limitations to measuring ET at large scales using ground-based methods have motivated the development of satellite remote sensing techniques. The purpose of this work is to evaluate the accuracy of the SEBAL algorithm for estimating surface turbulent heat fluxes at regional scale, using 28 images from MODIS. SEBAL estimates are compared with eddy-covariance (EC) measurements and results from the hydrological model MGB-IPH. SEBAL instantaneous estimates of latent heat flux (LE) yielded r²= 0.64 and r²= 0.62 over sugarcane croplands and savannas when compared against in situ EC estimates. At the same sites, daily aggregated estimates of LE were r² = 0.76 and r²= 0.66, respectively. Energy balance closure showed that turbulent fluxes over sugarcane croplands were underestimated by 7% and 9% over savannas. Average daily ET from SEBAL is in close agreement with estimates from the hydrological model for an overlay of 38,100 km² (r² = 0.88). Inputs to which the algorithm is most sensitive are vegetation index (NDVI), gradient of temperature (dT) to compute sensible heat flux (H) and net radiation (R_n). It was verified that SEBAL has a tendency to overestimate results both at local and regional scales probably because of low sensitivity to soil moisture and water stress. Nevertheless the results confirm the potential of the SEBAL algorithm, when used with MODIS images for estimating instantaneous LE and daily ET from large areas.

1. Introduction

Evapotranspiration (ET) is one of the most important regulating factors of climate, at both local and global scales, linking energy, climate and hydrology. ET is traditionally calculated as a difference of water-balance terms, whilst it is measured by instruments such as lysimeters or eddy-correlation (EC) systems. However such methods are limited in that they provide values of ET at specific sites and not at a regional or larger scale. Because the costs of installing complex measurement systems are high, and human and economic resources scarce, methods have therefore been developed in recent years which provide rapidly available hydrological data over large areas based on remote sensing.

To estimate energy fluxes between surface and atmosphere, a number of algorithms have been developed, including METRIC (Mapping EvapoTranspiration with Internalized Calibration) [1], Alexi (Atmosphere-Land Exchange Inverse) and DisAlexi (Disaggregated Atmosphere-Land Exchange Inverse) [2], SEBS(Surface Energy Balance System) [3], S-SEBI(Simplified Surface Energy Balance Index) [4], SEBAL (Surface Energy Balance Algorithm for Land) [5,6],T-SEB(Two-Source Energy Balance) [7], R-SEB(Aerodynamic Resistance Surface Energy Balance) [8], among others. The aim of such algorithms is to estimate latent heat fluxes (LE) and ET at both local and regional scales, using remote sensing with complementary meteorological data. Specifically, results obtained by SEBAL are regarded as consistent with field measurements [9–15]. Other research has reported results that are less clear [16,17].

SEBAL was developed conceptually for local and regional scale applications using Landsat images [11] with spatial resolution of 30 m and a temporal resolution of 16-day. However following the launch of other sensor systems, especially the MODIS sensor with 500 m of spatial resolution and a twice-daily temporal resolution, monitoring biogeochemical processes at the earth’s surface has gained a new boost. The availability of remote sensing data at different temporal, spatial and spectral resolutions allows new methodologies to be proposed for the purpose of monitoring and understanding the characteristics of complex environmental systems. The great advantage of MODIS data is the temporal resolution, since they can be used to estimate energy fluxes at regional, continental and global scales at daily time intervals [18], which is not possible with sensors such as Landsat TM and ETM, with broad images and less frequent revisits [19]. Much research has reported the influence of the remote sensing scale factor on ET estimates, where higher spatial resolution can contribute significantly to increased estimate accuracy [20–22], whilst lower resolutions can result in loss of important information for environmental applications [23,24]. Compared to Landsat TM or ETM spatial resolution, the moderate spatial resolution from MODIS sensor when used in surface energy balance algorithms at regional scales can contribute significantly to increase errors in ET estimates due to scale factor. However the high radiometric sensitivity of MODIS’s 36 spectral bands associated with geometric and atmospheric corrections can compensate for lower spatial resolution. ET estimates at continental and global scales based on moderate spatial resolutions have been shown to be accurate [25–32].

The purpose of this study was to estimate energy fluxes between surface and atmosphere for clear-sky days (days that are not affected by cloud cover at any time of the day) using SEBAL algorithm and Terra MODIS remote sensing datasets in tropical biomes in Brazil. Results were validated at local scale using EC data and at the spatial (regional) scale using hydrological modelling. The SEBAL algorithm [5,6,11], primarily driven by remote sensing data, was used to estimate ET at the regional scale, whilst the MGB-IPH (Large Basins Hydrological Model) hydrological model [33] was used to estimate ET at basin scale.

2. Methodology

2.1. Surface Energy Balance Algorithm for Land

2.1.1. Model Description

SEBAL uses multispectral remote sensing data associated with complementary meteorological data to estimate instantaneous surface energy balance components. The requirement for multispectral and thermal information means that the algorithm can only be used with images acquired on cloud-free days. Instantaneous LE (W·m⁻²) was estimated as the residual in the energy-balance equation (Equation (1)), as the difference between net radiation (R_n: W·m⁻²), soil heat flux (G: W·m⁻²) and sensible heat flux (H: W·m⁻²). A full description of the algorithm is given by Bastiaanssen et al. [5,6] and Allen et al. [11].

$$\text{LE}={\text{R}}_{\text{n}}-\text{G}-\text{H}$$

(1)

Surface albedo (α_s) was estimated using MODIS daily surface reflectance [34] (Equation (2)) while surface emissivity (ε_s) was estimated using MODIS 16-days NDVI composition [35] (Equation (3)).

$${\alpha }_{\text{s}}={\sum }_{\text{b}=1}^7\left[{\rho }_{\text{s},\text{b}}{\omega }_{\text{b}}\right]$$

(2)

$${\varepsilon }_{\text{s}}=1.009+0.047\hspace{0.17em}\text{ln}(\text{NDVI})$$

(3)

in which ρ_s,b is the at-surface reflectance for band “n” and ω_b is the weighting coefficient representing the fraction of at-surface solar radiation occurring within the spectral range represented by a specific band [34].

From the residual in the instantaneous energy-balance equation and the evaporative fraction (EF) (Equation (4)) the daily ET (ET_24h: mm·day⁻¹) were estimated (Equation (5)). EF has an important characteristic which is its regularity and constancy in cloud-free days [36]. Thus, its instantaneous value can be taken as the daily mean value, so that the spatial variability in daily ET can be predicted over large scales.

$$\text{EF}=\frac{\text{LE}}{{\text{R}}_{\text{n}}-\text{G}}$$

(4)

$${\text{ET}}_{24\text{h}}=\frac{86400\hspace{0.17em}\text{EF}\hspace{0.17em}{\text{R}}_{\text{n},24\text{h}}}{\lambda }$$

(5)

in which λ is the latent heat of evaporation (J·kg⁻¹) and R_n,24h (W·m⁻²) is the average daily net radiation balance (Equation (6)) estimated using a sinusoidal function [37], assuming that during the night R_n and the average daily G is zero [36,38].

$${\text{R}}_{\text{n},24\text{h}}=\frac{2{\text{R}}_{\text{n}}}{\pi \text{sin}\left[\left(\frac{{\text{t}}_{\text{overpass}}-{\text{t}}_{\text{rise}}}{{\text{t}}_{\text{set}}-{\text{t}}_{\text{rise}}}\right)\pi \right]}$$

(6)

where t_overpass (h) is the time at which the image was acquired and t_set (h) and t_rise (h) are the times of sunset and sunrise, respectively. The daytime hours were calculated based on latitude and the day of the year.

2.1.2. Remote Sensing Input Data

To run SEBAL we used Terra MODIS remote sensing datasets. The datasets used to estimate the surface-energy balance can be classified into: (a) daily surface reflectance (MOD09 GHK C5 [39]), (b) daily land surface temperature (MOD11 C5 [40]) and (c) 16-day vegetation indices (MOD13 C5 [41]). MODIS is operating onboard Terra and Aqua satellites. A ±55° scanning pattern at 705 km altitude achieves a 2,330 km swath that provides global coverage every one to two days. Terra has a 10:30 am equator over-passing time. MOD09GHK gives estimates of spectral reflectance from the surface in seven short-wave spectral bands. MOD11A1 corresponds to land surface temperature (T_s) on cloud-free days. Pixels are selected using a cloud mask (MOD35L2) which excludes pixels affected by aerosols and clouds. The accuracy of the T_s data is better than 1K, and is 0.5K in most cases [40]. However in pixels contaminated by cloud cover and heavy aerosols the accuracy may vary from 4 to 11 °C, since the cloud-cover mask does not discard all affected pixels, particularly those near cloud edges [42]. MOD13A1 gives a normalized difference vegetation index (NDVI), which provide information on vegetation spatial and temporal variations, as well as conditions of photosynthetic activity and biophysical and phenological changes [43]. SEBAL was used on 28 cloud-free days between February and November 2001, 19 of which were in the dry season (usually lasting from April to September) and 9 in the wet season (usually from October to March). Besides the cloud cover which made selection of MODIS images difficult in the wet season, there were also many EC data missing, limiting the application of the model for a larger number of days.

2.2. The Hydrological Model MGB-IPH

2.2.1. Model Description

The MGB-IPH is a distributed hydrological model developed for large basins with drainage areas in excess of 10,000 km² [33]. This model calculates the complete water balance at daily or monthly time intervals. The model runs using a regular grid with spatial resolution ranging from 5 to 20 km. The grid-cells are connected by channels representing the drainage network [44]. Each grid-cell is divided into classes that combine soil type and vegetation, which are called hydrological response units (HRU) [45,46] to account for the fractional contributions from different physical characteristics within each grid-cell. In the MGB-IPH model, evaporation and transpiration are calculated separately based on Penman-Monteith approach [47] (Equation (7)).

$${\text{ET}}_{24\text{h}}=\frac{\mathrm{\Delta }\left({\text{R}}_{\text{n}}-\text{G}\right){+\rho }_{\text{a}}{\text{C}}_{\rho }\left({\text{e}}_{\text{s}}-{\text{e}}_{\text{a}}\right)/{\text{r}}_{\text{a}}}{\mathrm{\Delta }+\gamma \left(1+\raisebox{1ex}{${\text{r}}_{\text{s}}$}\!\left/ \!\raisebox{-1ex}{${\text{r}}_{\text{a}}$}\right.\right)}$$

(7)

where ET_24h (mm·day⁻¹) is the daily evapotranspiration, Δ (kPa·K⁻¹) is the gradient of saturated vapour pressure, R_n (MJ·m⁻²·day⁻¹) is the daily net radiation, G (MJ·m⁻²·day⁻¹) is the daily soil heat flux, ρ_A (kg·m⁻³) is the air density, C_p (MJ·kg⁻¹·K⁻¹) is the specific heat of air at constant pressure, e_s (kPa) and e_a (kPa) are the saturated vapour pressure and actual vapour pressure, respectively, γ (0.066 kPa·K⁻¹) is the psychometric constant, whilst r_s (s·m⁻¹) and r_a (s·m⁻¹) are the surface and aerodynamic resistance, respectively.

It is assumed that evaporation of water intercepted by the canopy leads soil evaporation and plant transpiration. Evaporation of soil water occurs subsequently. If there is still a demand for evaporation, water is evaporated from a second soil layer. The maximum depth of canopy intercepted water (EI_max: mm·day⁻¹) is determined for each HRU as a function of leaf area index (LAI). After estimating evaporation, the remaining evaporative demand fraction (f_DE) (Equation (8)) is met by plant transpiration for each vegetation cover type, which is used as a correction factor in calculating the plant transpiration (Equation (9)).

$${\text{f}}_{\text{DE}}=\frac{\text{EIP}-\text{EI}}{\text{EIP}}$$

(8)

$${\text{ET}}_{24\text{h}}={\text{f}}_{\text{DE}}*\frac{\mathrm{\Delta }\left({\text{R}}_{\text{n}}-\text{G}\right)+{\rho }_{\text{a}}*{\text{C}}_{\rho }\left({\text{e}}_{\text{s}}-{\text{e}}_{\text{a}}\right)/{\text{r}}_{\text{a}}}{\mathrm{\Delta }+\gamma \left(1+\raisebox{1ex}{${\text{r}}_{\text{s}}$}\!\left/ \!\raisebox{-1ex}{${\text{r}}_{\text{a}}$}\right.\right)}*\frac{1}{\lambda *\rho \text{w}}$$

(9)

where EIP (mm·day⁻¹) is the potential evaporation of intercepted water, EI (mm·day⁻¹) is the intercepted water by the canopy and ρ_w (kg·m⁻³) is the density of water. In this model, r_a depends only on the canopy height (m) and wind speed (m·s⁻¹) whilst r_s is characteristic of each vegetation type but which vary according to the restriction of soil moisture (W: mm·day⁻¹) [48]. It is assumed that soil conditions do not restrict ET if W is above 50% of soil water capacity (W_m: mm·day⁻¹) [49]. In this case, r_s is regarded as a minimum value typical of vegetation unaffected by soil moisture conditions. If soil water storage lies between wilting point (W_PM: mm·day⁻¹) and the beginning of the effect on r_s (W_L: mm·day⁻¹), r_s increases (Equation (10)). If W is less than the value W_PM, the restriction is at its greatest and ET is zero.

$${\text{r}}_{\text{s}}={\text{r}}_{\text{sm}}*\frac{{\text{W}}_{\text{L}}-{\text{W}}_{\text{PM}}}{\text{W}-{\text{W}}_{\text{PM}}}$$

(10)

where r_sm (s·m⁻¹) is surface resistance for a given W. Full details of soil water balance in the MGB-IPH model are given by Collischonn et al. [33]. The MGB-IPH hydrological model was applied successfully in several South American large-scale river basins [50–53].

2.2.2. Input Data and Model Validation

The hydrological model calibration and validation were performed using data from 273 rain gauges distributed over the basin were used for the analyses. The data were obtained from the Brazilian Water Agency [54] and from the Integrated System for Water Resource Management of the State of São Paulo [55]. Meteorological data were obtained from 14 data collection platforms provided by the National Institute for Spatial Research in Brazil [56]. To characterize the HRU, we used soil type information from the Radam Brasil [57] as well as land use data derived from classification of LANDSAT 7 ETM+ images. LANDSAT images were classified into four distinct groups: (i) open water; (ii) forest and reforested areas; (iii) agricultural areas; and (iv) pasture and bare soil. Soil types were regrouped into three basic classes according to water storage capacity: (i) high; (ii) medium; or (iii) low. After combining information on soil type and land use, six HRU’s were defined (10 × 10 km cell):

• pasture, grassland and cropland areas with soils of medium storage capacity (HRU 1);
• cropland areas with soils of high storage capacity (HRU 2);
• soils with low storage capacity (HRU 3);
• forest and reforested areas on soils with medium storage capacity (HRU 4);
• pasture, grassland and bare soil on soils with high storage capacity (HRU 5);
• open water surfaces (HRU 6).

The model was run using a daily time interval. The period from 1970 to 1980 was used to calibrate the model, and the period from 1981 to 2001 was used for model validation [58]. Nash-Sutcliffe coefficients of model fit were greater than 0.7 in all sub-basins in both stream flow calibration and validation periods, whilst errors in volume were less than 0.5% at calibration and 7% at validation [59]. The calibration and validation procedures and results are fully presented in Collischonn et al. [58] and Tucci et al. [59].

2.3. Site Description

The flux sites were set up in 2000 installed at a savanna ecosystem (PDG site) and a sugarcane cropland (USE site). Turbulent fluxes (H and LE) were measured at heights of 7 m (USE) and 21 m (PDG). For the 7 m height at USE, the maximum contributing area has a footprint of between 40 and 60 m, whilst for the 21 m height at PDG, it was between 120 and 170 m [58–60]. The accumulated flux for neutral conditions shows that more than 80% of measured flux at PDG site originates within 500 m, while at USE site originates within 1,500 m [58]. Turbulent flux measurements of H and LE at the USE site had calibration problems between August and December 2001 while some of the equipment o the tower installed at the PDG site had missing calibrations during January and September 2001 [58]. Ancillary meteorological data measurements needed for the analysis included R_n and G. Detailed descriptions of the equipment used, measurement procedures and results at both PDG and USE sites can be found in Juárez [60], Cabral et al. [61] and Rocha et al. [62].

2.4. Study Area

The study area lies between the Brazilian States of São Paulo and Minas Gerais. The Rio Grande basin is located between latitudes 19°15´S and 23°00′S and longitudes 43°30′W and 51°00´W. Most of the basin is covered by a savanna vegetation, locally called Cerrado. This natural vegetation was extensively replaced during the last centuries by intensive production of sugarcane, soya and by pastures. Savanna ET has a strong seasonality, varying between 6 mm·day⁻¹ in the wet season to 1 mm·day⁻¹in the dry season. In agricultural areas, ET varies according to the type of the cropland and its cultivation cycle, with significant variations due to vegetation structure that may range from bare soil to fully developed canopy cover.

The MGB-IPH was applied to the Rio Grande basin (145,000 km²) while SEBAL was applied over an area of 2° × 2° (45,000 km²) selected to cover the flux tower sites. The overlay between MGB-IPH and SEBAL is about 38,100 km² (Figure 1).

2.5. Data Analysis

Measured data (provided at half-hour) were averaged to 1-hour and to daily periods to compare with model-predicted data. We did not gap-fill with a model and we excluded any given daily average if fewer than 75% of data time steps were available. The performance of the models were evaluated based on values of root mean square error (RMSE), mean absolute error (MAE) and coefficient of determination (r²) for model-predicted versus measured data. The RMSE is the overall error in the predictions relative to the actual measured value while MAE measures the average magnitude of the errors in a set of predictions, without considering their direction. The r² assesses how well the model reflects the variations in the data.

3. Results and Discussions

3.1. Validation of SEBAL’s Instantaneous Energy Fluxes

SEBAL estimations were spatially averaged over a 3 × 3 1-km pixels window centred on each flux tower to achieve spatial representativeness of the measured data for each site. Instantaneous values of R_n derived from SEBAL near mid-day (Terra MODIS has a ∼11:30 am local over-passing time) significantly explained over 80% of the variation of measured values at the USE (r² = 0.88, p < 0.05) and at the PDG sites (r² = 0.83, p < 0.05) (Figure 2(a)). There was a very marked seasonality, associated with the variation in solar radiation at the earth’s surface. RMSE was 37 W·m⁻² at USE and 47 W·m⁻² at PDG. Largest errors in the estimation of R_n occurred during the wet season, when days may have been partially cloudy or when clouds were present but not detected in MOD11A1. Instantaneous values of G at USE varied between 12 and 58 W·m⁻² in dry and wet seasons, respectively, and at PDG between 8 and 24 W·m⁻². Estimated values were moderately correlated with observed data (Figure 2(b)), varying between 40 and 120 W·m⁻² at USE, and between 29 and 70 W·m⁻² at PDG, in dry and wet seasons respectively. In absolute values, G was overestimated on average by 35 and 41 W·m⁻². Estimates of H showed greater uncertainties relative to observed data. The mean instantaneous flux of H estimated by SEBAL at USE site was 217 ± 68 W·m⁻², whilst the mean value given by the EC system was 175 ± 55 W·m⁻². At the PDG site, the estimated mean was 170 ± 71 W·m⁻², compared with the observed EC mean of 179 ± 51 W·m⁻². The RMSE was 36 and 38 W m⁻² at USE and PDG sites, respectively. We found a r² of 0.60 (p < 0.05) and 0.52 (p < 0.05) at USE and PDG sites when comparing modelled H values with the in situ EC measurements (Figure 2(c)). Discrepancies in estimates of H can be explained by the combination of a series of factors related to simplifications and uncertainties introduced in algorithms to predict ET [17,63]. However these discrepancies are similar to results found by other studies [17,63–68] when validating different models for turbulent flux estimation under different conditions of land use and land cover. Although LE fluxes estimated using SEBAL showed the correct seasonal pattern, an overestimation was found. At the USE site the mean overestimation was 23 W·m⁻² with r² = 0.64 (p < 0.05), whilst at PDG overestimation was approximately 70 W·m⁻² with r² = 0.62 (p < 0.05), when comparing modelled LE values with the in situ EC measurements (Figure 2(d)). The mean estimated LE at USE was 271 ± 100 W·m⁻², whilst for the EC system the mean was 287 ± 87 W·m⁻², with RMSE 34 W·m⁻². At PDG, mean LE was estimated as 386 ± 93 W·m⁻² compared with a mean of 285 ± 122 W·m⁻² from the in situ EC data, with RMSE of 81 W·m⁻².

The energy balance closure, given by comparing the turbulent fluxes (H+LE) estimated by SEBAL with the observed available energy (R_n–G), shows an underestimation of about 7% in turbulent flux at the USE site and 9% at PDG, with r² = 0.88 (p < 0.05) and r² = 0.85 (p < 0.05) respectively (Figure 3(a)). Since the measurements and estimates of R_n, G, H and LE are found by different methods, the energy balance closure shows whether or not there is agreement between these values. The mean value of H+LE estimated by SEBAL for the USE site was 488 ± 78 W·m⁻² whilst the mean of measured R_n–G was 516 ± 68 W·m⁻². At PDG, the mean H+LE given by SEBAL was 536 ± 89 W·m⁻² compared with a mean of measured R_n–G, 589 ± 109 W·m⁻². From an analysis of energy balance closure, Juarez [60] found underestimations of 24% and 23% in observed turbulent fluxes at the sites USE and PDG respectively, suggesting that there may be limitations inherent in the EC measurement system. It is not always possible to reach a satisfactory result because of measurement errors, particularly in the EC technique, which normally underestimates turbulent fluxes because of the difference between the levels at which R_n sensors are mounted (some meters above the surface) compared with sensors for G (some centimeters below the surface), and because of surface heterogeneity [69–71].

3.2. Control of SEBAL’s Instantaneous Latent Heat Flux

To understand what drives estimations of ET using SEBAL algorithm, coefficients of determination (r²) were calculated between instantaneous LE (on MODIS overpass) and the main inputs to the model (albedo, land surface temperature (T_s) and NDVI) and intermediate variables (surface emissivity (ε_s), roughness length (Z_om), gradient of temperature (dT), wind friction velocity (u_*), aerodynamic resistance (r_a) and net radiation (R_n)). Figure 4 shows: (1) the importance of NDVI for estimating LE (disregarding ε_s and Z_om since these are estimated as functions of NDVI, which explains the high correlation between these variables and LE), (2) the importance of precision in the choice of “anchor” pixels to determine dT and (3) the influence of R_n in controlling LE.

The NDVI uses spectral information at near infrared (NIR) wavelengths (between 0.7 and 1.0 μm) to measure terrestrial photosynthetic activity and structural variations in the canopy [72]. Variations in LE are quantitatively associated with spatial and temporal changes in the vegetation canopy, since under normal conditions of water and nutrient availability there is a strong negative correlation between T_s and NDVI. Although NDVI and LE are highly correlated, two factors are limiting: (i) the saturation of MODIS NDVI in areas where vegetation is quite dense [73]; and (ii) the low correlation between NDVI and water stress. NDVI commonly shows an asymptotic saturation in areas with high biomass indices [74]. These saturation effects can mask changes in vegetation dynamics and biophysical processes. In addition, the ability of NDVI to convey information about vegetation water content is limited, since the index utilizes wavelength related to internal leaf structure and to photosynthetic processes, and does not give a measure directly related to the quantity of water in the leaf and to vegetation water stress [75]. The quantity of water available in a leaf’s internal structure largely controls infrared short-wave (SWIR) wavelengths (between 1.2 and 2.6 μm) [76]. This limitation of NDVI in relation to water stress could be partly overcome by using vegetation indices that are more sensitive to canopy water availability and leaf water content or by incorporating the relation between soil moisture and surface resistance. One alternative is to use the normalized difference water index (NDWI) [77]. The combination of NIR with SWIR removes variations introduced by internal leaf structure and quantities of dry material, improving accuracy in monitoring water stress and giving a measure that is more indicative of water availability than NDVI [78–80]. Predicting turbulent fluxes of H constitutes the main step in SEBAL’s partition of energy, incorporating many of the uncertainties generated by the algorithm. The subjective choice of hot and cold pixels determines the accuracy of H and consequently LE, as shown by the high correlation between dT and LE. It is difficult to avoid this subjectivity in the selection of reference pixels. Allen et al. [1,81] suggest that the choice of warm pixels should be determined not from T_s alone, but also with reference to its relationships with other parameters, particularly vegetation indices like NDVI or LAI.

3.3. Scaling SEBAL’s Instantaneous Latent Heat Flux to Daily Evapotranspiration

Estimation of daily ET requires the integration of instantaneous fluxes R_n, G and LE predicted by SEBAL using EF and R_n,24h. As a starting point, the accuracy of estimated EF was assessed relative to the measured data. Next, the validation of the EF as a constant indicator of energy partitioning throughout the day (taken as the period between 7:00 h and 18:00 h) was examined. When the estimated instantaneous EF was compared with observed values, an overestimation by 12% and 18% was found at USE and PDG sites, respectively, most probably as a consequence of underestimation of turbulent fluxes measured by the EC system and in SEBAL’s energy balance closure, causing significant reduction in EF (Figure 3(b)). Figure 5 shows estimates of instantaneous energy fluxes found using SEBAL and the daily energy cycle observed at EC monitoring sites for two days during wet and dry seasons, together with diurnal variability in EF at those sites.

When instantaneous data were converted to daily ET, the estimated seasonal variation in R_n,24h (Figure 6(a)) and ET_24h (Figure 6(b)) agreed, respectively, with in situ measurements. Estimated values of R_n,24h significantly explained over 85% of the variation of measured data ate USE site (r² = 0.85, p < 0.05) and at the PDG site (r² = 0.86, p < 0.05). R_n,24h was overestimated by about 23 W·m⁻² and 26 W·m⁻² at USE and PDG sites, since this method (based on Equation (6)) does not correct for the radiation losses observed throughout the day as a function of cloud cover and atmospheric aerosols. Furthermore, this sinoidal model reduces R_n to zero during the night-time when effectively R_n becomes negative. As R_n goes to zero, we are assuming that ET goes to zero as well. According to Fisher et al. [82] without large wet surfaces, night-time ET is minimal because both photosynthesis and transpiration approaches zero as radiation diminishes. Using this simple scheme to predict daily R_n we are trying to overcome the need for ground information as model input.

With regard to ET_24h derived from SEBAL, estimated values significantly explained between 66% and 76% of the variation of measured data at the USE site (r² = 0.76, p < 0.05) and at the PDG site (r² = 0.66, p < 005), respectively. At USE site, estimated ET_24h varied between 5.9 mm·day⁻¹ in the wet season and 1.2 mm·day⁻¹ in the dry season, whilst observed values varied from 5.0 mm·day⁻¹ to 1.5 mm·day⁻¹ in the wet and dry seasons. At PDG site, estimated ET varied from 6.6 mm·day⁻¹ in the wet season to 2.3 mm·day⁻¹ in the dry season, with observed values varying from 6.5 mm·day⁻¹ to 1.2 mm·day⁻¹ in wet and dry seasons. Estimated and observed values had similar standard deviations, but ET was overestimated by 0.8 mm·day⁻¹ at USE, and by 1.15 mm·day⁻¹ at PDG, on average, reflecting the overestimation of EF and R_n,24h.

3.4. Estimation of Daily Evapotranspiration by Hydrological Modelling

After calibration and validation the MGB-IPH hydrological model, daily ET estimations were compared with ET measured by EC at the USE and PDG sites. PDG site was compared with HRU 4 (forest and reforested areas on soils with medium storage capacity) and the USE site with HRU 2 (cropland areas with soils of high storage capacity), lying in the same area as the flux towers. At the PDG site the average ET estimated by the MGB-IPH model was 2.5 ± 1.4 mm·day⁻¹, whilst the measured ET was 2.6 ± 1.2 mm·day⁻¹, with RMSE equal to 0.72 mm·day⁻¹ (Figure 7(a)). Average ET estimated by the MGB-IPH model at the USE site was 2.7 ± 1.1 mm·day⁻¹, while the EC-measured ET at the same site was 2.6 ± 1.1 mm·day⁻¹, with RMSE 0.51 mm·day⁻¹ (Figure 7(b)). Compared to measured data, estimated ET_24h significantly explained over 64% of the variation of measured data at the PDG site (r² = 0.64, p < 0.05) and at the USE site (r² = 0.67, p < 0.05). We also found that the MGB-IPH model overestimated ET at the USE site by about 0.06 mm·day⁻¹, and underestimated it by about 0.13 mm·day⁻¹ at the PDG site.

3.5. Evapotranspiration at Basin Scale

A crucial point in spatial modelling of the energy balance is the capacity of thermal images to distinguish patterns of land use and land cover, since T_s is one of the most important parameters for estimating ET, followed by soil moisture and the fraction of vegetation in each pixel [83]. The great advantage in using MODIS data for estimating ET is the high accuracy of surface temperature images [84,85] associated with the spatial variability of this process at the regional scale. Two analyses were made to validate ET at the basin scale: (i) average ET_24h estimated by SEBAL was compared with ET_24h estimated by MGB-IPH, in the overlay area for which the two models were applied (an area of about 38,100 km²), and (ii) average ET_24h estimated by the two models were compared for different land use and land cover types, taking the hydrological response units (HRU) in the MGB-IPH as the analytical classes.

Modelled mean ET_24h for the area of 38,100 km² explain over 88% of the variation observed in the MGB-IPH model (r² = 0.88, p < 0.05) (Figure 8) with a RMSE of 0.48 mm day⁻¹. The SEBAL overestimated ET_24h by about 0.45 mm·day⁻¹ when compared with the model MGB-IPH, a difference of 14%. In the wet season the overestimation was about 0.16 mm·day⁻¹; while during the dry season the overestimation was about 0.60 mm·day⁻¹. In the wet season the mean ET_24h from the MGB-IPH model was 4.6 ± 0.6 mm day⁻¹ in comparison to an average of 4.8 ± 0.6 mm day⁻¹ from the SEBAL model. In the dry season the mean ET_24h from the MGB-IPH model was 2.1 ± 0.8 mm·day⁻¹ in comparison to an average of 2.6 ± 0.8 mm·day⁻¹ from the SEBAL model. In both seasons, estimations of ET_24h from SEBAL were higher than those given by MGB-IPH, although the difference was greatest in the dry season, related to SEBAL overestimations of instantaneous LE, EF and R_n24h.

Land use and land cover in the Rio Grande basin is a complex mosaic of agricultural and natural areas. Of the six HRU’s used in MGB-IPH, the following four can be found in the overlay area: (1) pasture, fields and agricultural areas in areas of medium infiltration capacity, covering 22% of the area; (2) agricultural areas in soils with high infiltration capacity, covering 36% of the area; (4) forest and reforested areas on soils with medium infiltration capacity, 16% of area; (5) pasture, fields and bare soil where soil has high infiltration capacity, 25% of area. Although SEBAL overestimated ET_24h relative to MGB-IPH, correlations between the two models in the four HRU’s were moderate to high, significantly explaining from 46% to 86% of the variation of measured values. Correlation was higher for agricultural areas, fields and pastures on soils with medium and high infiltration capacity; in forested areas, it was lower than for other units (

Table 1. Comparisons between daily evapotranspiration estimated by SEBAL and MGB-IPH. Labels: HRU 1 (pasture, grassland and cropland areas with soils of medium storage capacity); HRU 2 (cropland areas with soils of high storage capacity); HRU 4 (forest and reforested areas on soils with medium storage capacity); HRU 5 (pasture, grassland and bare soil on soils with high storage capacity).

		HRU 1	HRU 2	HRU 4	HRU 5
Area (Km2)		8,510	13,686	6,085	9,364
r2		0.86 (p < 0.05)	0.65 (p < 0.05)	0.46 (p < 0.05)	0.79 (p < 0.05)
MAE (mm day−1)		1.2	0.4	0.9	0.7
RMSE (mm·day−1)		1.2	0.7	1.1	0.9
Avg±StDev (mm day−1)	MGB-IPH	2.4 ± 1.7	2.6 ± 1.3	2.8 ± 1.2	2.5 ± 1.5
	SEBAL	3.7 ± 1.2	3.1 ± 1.2	3.7 ± 1.2	3.3 ± 1.2

).

Differences between correlations may be mainly the consequence of the complex system of land use and land cover and of the different input data and parameterizations used by the two models. In the MGB-IPH, the estimate of ET is based mainly on meteorological data obtained from field stations. Surface parameters are adjusted to agree with existing land use and land cover maps. Thus the variation in time and space of vegetation parameters such as albedo, LAI and r_s is determined by average values given in the literature, which probably do not accurately describe the best spatial variation in vegetation. For example, HRU 2, covering an area of about 13,600 km² and characterized by different agricultural cultivation, uses the same surface parameters for all of them. By contrast, the more important input variables to SEBAL, such as NDVI, albedo and T_s more faithfully describe temporal and spatial vegetation processes, giving less attention to micrometeorological information which is difficult to obtain at a regional scale.

4. Concluding Remarks

Because of the large quantities of surface information obtained by MODIS sensor and its available products, estimates of energy fluxes of the Rio Grande basin given by SEBAL algorithm can be regarded as consistent. Instantaneous estimates of R_n, H and LE were significantly related to measured values, yielded a r² varying from 0.52 to 0.88. We found, however, that the equation for determining G overestimated results, with a RMSE ranging from 35 and 41 W m⁻², since this equation was specifically calibrated for semi-arid areas and may not accurately represent the local conditions of our study region. Discrepancies between predicted and measured H and LE may be associated with the low spatial resolution of MODIS products and high heterogeneity in land use and land cover within the same pixel. Results given by energy balance closure showed a r² of 0.88, suggesting a good agreement between available energy and turbulent fluxes, although the partition of this energy into H and LE does not yet show the same consistency of results as the energy balance suggests. We also found that inputs to which SEBAL is most sensitive are vegetation index (NDVI), gradient of temperature (dT) and net radiation (Rn).

Estimates of ET_24h based on SEBAL algorithm were overestimated by 0.80 mm·day⁻¹ at sugarcane croplands and 1.15 mm·day⁻¹at savannas when compared to recorded data by the EC system, probably in consequence of factors such as:

• The omission of a soil moisture constraint for a region of known water limitations.
• The description of vegetation water stress by NDVI, which is limited by the asymptotic saturation level in areas where biomass index is high and information about vegetation water content is difficult to obtain.
• The subjective determination of the gradient of temperature (dT) to estimate sensible heat fluxes (H).
• Limitations from the compensatory or cumulative errors entered through the residual energy balance, partitioning turbulent fluxes and estimating EF and R_n,24h.
• The omission of night net radiation (R_n) when it becomes effectively negative or even assuming that average daily soil heat flux (G) is zero can lead to overestimations.

Regarding to the hydrological model MGB-IPH, when compared with measured data yielded a r² of 0.64 at savannas and 0.67 at sugarcane croplands. Both models, SEBAL and MGB-IPH, provided similar results for average estimated ET over an area of 38,100 km² with a r² of 0.88. However it was found that SEBAL overestimated ET by 14% relative to MGB-IPH.

Comparing SEBAL and MGB-IPH to estimate daily ET, both methods have different spatial and temporal capabilities. MGB-IPH hydrological model can be run with different time-steps and levels of data although the overall data requirements is high. SEBAL algorithm has the advantage that the data requirements are low and spatial resolutions is high. The disadvantage is that it only can be applied to clear-sky days. To estimate seasonal ET the necessity of analysing many images can be an expensive problem. Remote sensing data extraction under cloud cover still remains a challenging task considering surface energy balance models to predict ET.

To conclude, results obtained from SEBAL algorithm confirm the potential of modelling surface energy fluxes on clear-sky days. The technique is not restricted to irrigated areas (as originally developed) but can be applied to a broad range of biomes. Despite the overestimation of R_n,24h and ET_24h, we argue that the strength of the SEBAL algorithm should be evaluated, not only by how closely the estimated results are from measured data, but by the ability to provide spatial information over large areas.