Evaluation of MOD16 Algorithm over Irrigated Rice Paddy Using Flux Tower Measurements in Southern Brazil

Abstract

Evapotranspiration (ET) is an important component of the hydrological cycle. Understanding the ET process has become of fundamental importance given the scenario of global change and increasing water use, especially in the agricultural sector. Determining ET over large agricultural areas is a limiting factor due to observational data availability. In this regard, remote sensing data has been used to estimate ET. In this study, we evaluated the Moderate-Resolution Imaging Spectroradiometer (MODIS) land surface ET product estimates (hereafter MOD16 ET – MODIS Global Terrestrial Evapotranspiration Product) over two rice paddy areas in Southern Brazil, through the ET measured using the eddy covariance technique (hereafter EC). The energy balance components were evaluated during fallow and flooded seasons showing latent heat flux dominates in both seasons. The results showed that MOD16 ET underestimated EC measurements. Overall, the RMSE (root mean square error) ranged between 13.40 and 16.35 mm 8-day⁻¹ and percent bias (PBIAS) ranged between −33.7% and −38.7%. We also assessed the ET (measured and estimated) main drivers, with EC yielding higher correlation against observed net radiation (R_n) and global radiation (R_g), followed by air temperature (Temp) and vapor pressure deficit (VPD), whilst MOD16 ET estimates yielded higher correlation against leaf area index (LAI) and fraction of photosynthetically active radiation (fPAR). The MOD16 algorithm was forced with meteorological measurements but the results did not improve as expected, suggesting a low sensitivity to meteorological inputs. Our results indicated when a water layer was present over the soil surface without vegetation (LAI around zero), the largest differences between EC measurements and MOD16 ET were found. In this period, the expected domain of soil evaporation was not observed in MOD16 ET physical processes partition, indicating the algorithm was not able to detect areas with high soil moisture. In general, the MOD16 ET product presented low accuracy when compared against experimental measurements over flooded rice paddy, suggesting more studies are necessary, in order to reduce uncertainties associated to the land cover conditions.

1. Introduction

The increasing competition for water associated to its scarcity is becoming a threat to sustainable development [1]. In a climate changing scenario associated with an increasing global water demand for cropland production, that currently accounts for approximately 70% of the world water consumption [2], water footprint quantification is fundamental to improve water management and mitigate future freshwater scarcity [3].

One of the most challenging issues in modern global agricultural management is the proper control of the amount of water made available for plantation. To efficiently determine how much water to apply in the crop requires the knowledge of the amount of water lost. The major natural processes involving water losses include evaporation from the surface and plant transpiration, also defined as evapotranspiration (ET). Therefore, a correct estimation of ET over irrigation-based areas is critical to improve productivity while reducing water use as well as to mitigate the impact of extreme events in agriculture, such as droughts [4]. Nevertheless, determining ET over large agricultural areas is a challenging task due to limited observational data availability. The eddy covariance (EC) is the state of art technique to estimate actual ET, but it is expensive and represents a region usually ranging from 100 to 1000 m of radius [5]. To overcome footprint limitations, remote sensing techniques have been largely used to estimate ET over large areas with high spatial and temporal resolution and high accuracy [6].

Currently, several models to estimate ET using remote sensing data are available, from local and regional scales (SEBAL [7], SEBS [8], METRIC [9], ALEXI [10]) to continental and global scales (MOD16 [11], ALEXI/DisALEXI [10], JPL-PT [12], GLEAM [13,14]), over a wide range of temporal scales. According to Biggs et al. [4] models can be classified in three different categories: (i) physical models based on vegetation indices and meteorological data, (ii) energy balance models based on land surface temperature and (iii) empirical/statistical models. Several studies were conducted to evaluate the accuracy of ET derived from remote sensing using different methods [15,16,17,18,19,20], showing an overall uncertainty between 10% to 30%, when compared to ET obtained by eddy covariance measurements or to water balance. A detailed discussion of remote sensing-based ET models can be found in Gowda et al. [21], Li et al. [22] and Biggs et al. [4].

A widely used algorithm to estimate ET at continental and global scale is MOD16 (MODIS Global Terrestrial Evapotranspiration Product). The algorithm was developed by Mu et al. [11] and Mu et al. [23] based on the Penman–Monteith equation, to monitor temporal and spatial variability of ET over global land surface areas. The algorithm combines MODIS (Moderate-Resolution Imaging Spectroradiometer) 8-day and 16-day remotely sensed inputs for vegetation representation and daily Modern-Era Retrospective analysis for Research and Applications (MERRA) meteorological reanalysis data. Most of the experimental data used for calibration and validation were obtained from sites over North America (AmeriFlux and FluxNet) between 2000 and 2010 [23]. Independent validation of MOD16 ET estimation, often performed using ET obtained by eddy covariance measurements [24], have been carried out in several land cover and climatological conditions around the world. High correlations were yielded for grassland in the United Kingdom [25]; for forests in arid and polar regions in Asia [26]; for cropland and grassland in Europe [18]. Nevertheless, several authors have found moderate to low correlations for forests [18,27,28], woodland savannas [15,18], arid and semi-arid regions in the United States [29]. Only a few studies have assessed the accuracy of MOD16 ET algorithm over irrigated agricultural areas, usually focusing on rainfed cropland [4,30]. Tang et al. [28] perform a study on maize, wheat and cotton over an irrigated area, finding an underestimation of MOD16 for several crops. Hu et al. [18] assess the MOD16 algorithm performance on wheat and barley croplands, finding reasonable results for Europe. Kim et al. [26] validate the MOD16 over two (irrigated) rice paddy sites, concluding that at one site ET was underestimated whereas at the other site ET was overestimated. Biggs et al. [4] highlighted that, despite being recommended mainly for rainfed systems and areas where vegetation cover is high, MOD16 ET algorithm needs further assessment in irrigated cropland areas with high soil evaporation.

The accuracy of MOD16 ET estimates depends mainly on factors related to the model structure and its parameterization and accuracy of input data [31,32]. The Penman–Monteith equation is known for its robust framework for ET estimation because (i) it involves the main drivers of evapotranspiration; (ii) it provides an energy constraint on the evaporation rate and (iii) estimated ET fluxes are not exceedingly sensitive to any of the inputs [33]. Despite its strong theoretical basis, the Penman–Monteith equation application at larger spatial scales is hindered by micrometeorological data requirements, which are difficult to obtain over large areas [33], and by issues of scale and footprint [34].

To ensure water security and to enhance the ability to simulate hydrological and climate feedbacks in a climate change scenario, there is a need to represent ET with high precision and accuracy, improving differentiation of water use and water stress among different crops, species and ecosystems [6]. Furthermore, Baldocchi et al. [35] highlight that long-term research of energy and water fluxes over flooded rice paddy, using eddy covariance measurements, are needed to understand ET patterns and to improve water management practices. In the context of this research, the main goal was to address the question of what drivers control the ET process, performing an assessment of surface observed data and remotely sensed data, over flooded areas. To perform this evaluation, we conducted an analysis using ET obtained by EC data: (1) evaluation of the energy partitioning between latent (LE) and sensible (H) heat fluxes; (2) assessment of the seasonal evaporative fraction (Λ) pattern during fallow and flooded rice paddy cycle; and (3) assessment of the main control drivers of ET. We also sought to assess whether global ET models, specifically the MOD16 algorithm, can physically reproduce this pattern and estimate accurately ET over flooded areas. Our analysis focused on: (4) validation of MOD16 ET estimates using EC measurements during flooded and non-flooded periods; (5) evaluation of control drivers in MOD16 ET estimation; (6) comparison between MOD16 ET estimates against observed meteorological variables to address whether the algorithm can reproduce the physical process of ET; and (7) analysis of residual errors attributed to meteorological reanalysis inputs.

2. Materials and Methods

2.1. Site Description

Brazil is one of the top ten worldwide rice producers, with approximately 12.5 million tons per year and an average production of 5424 kg ha⁻¹ and Rio Grande do Sul State, located in Southern Brazil, is responsible for more than 77% of Brazilian rice production [36]. In Rio Grande do Sul, rice production occurs generally in natural flooded regions due to its large water consumption [37]. The rice paddies are cultivated once a year, during spring and summer seasons, with soil remaining in fallow for the rest of the year. Atmospheric datasets collected from two rice paddy sites located in Rio Grande do Sul State, Brazil, were used in this research (Figure 1). An experimental site was located in Paraíso do Sul city (PRS), at 29°44′39.6″ S and 53°8′59.8″ W, altitude of 108 m (hereafter PRS site), in an area rice paddy of approximately 50 hectares (ha), with wind direction predominantly from the southeast. An experimental site was located in Cachoeira do Sul city (CAS), at 30°16′37.59″ S and 53°8′52.25″ W, altitude of 40.5 m (hereafter CAS site) in an area of approximately 1000 hectares, with wind direction predominantly from the east. In both sites, rice paddies are cultivated over large lowland areas (floodplains), using a flooding irrigation system. After rice harvesting, soil remains with spontaneous vegetation (fallow period), without crop rotation. The experimental sites are around 60 km apart.

The climate in the region is classified, according to Köppen climate classification, as humid subtropical climate (Cfa), with hot summer and without dry season (daily mean temperature greater than 22 °C) [38]. The climatological precipitation for Cachoeira do Sul is 1477.1 mm year⁻¹ (1961–1990 according to Brazilian National Meteorological Service (INMET) in http://www.inmet.gov.br/portal/index.php?r=clima/normaisClimatologicas). The precipitation is distributed throughout the year between 83 and 157 mm by month. There is no registration of climatological precipitation for Paraíso do Sul. The land use and land cover in the study areas are dominated by grassland and evergreen broadleaf forest, whilst cropland areas (cultivated with flooded rice paddy) are mosaicked with natural vegetation. The soil type at PRS is silt loam and in CAS is sandy clay loam, both typical flooded soil where rice is grown in Rio Grande do Sul state. More details about physical properties of the soil at the experimental sites are described in Timm et al. [37] for PRS and in Diaz et al. [39] for CAS.

The period of analysis for the PRS experimental site was from July 2003 to July 2004. Rice was cultivated at the end of November 2003 and harvested at the beginning of April 2004. The measurements for the CAS experimental site span from January 2011 to December 2014. Rice was cultivated between October and November and harvested between March and April from the next year. The field was flooded until harvest, to a depth of 0.1 m. The PRS site was flooded around 30 days after seeding, while the CAS site was flooded around 30 days before seeding. During the other periods of the year, the areas were maintained in fallow and the soil remained near saturation throughout the time due to high precipitation rates.

2.2. Meteorological and Flux Tower Measurements

During the period of this study, the PRS and CAS experimental sites were equipped with flux towers measuring atmospheric variables and surface fluxes. Air temperature, wind speed components and H₂O concentration were measured in high frequency to compute the vertical fluxes obtained by the EC technique. The PRS and CAS flux towers were equipped the sensors described in

Table 1. Sensors installed at the flux towers used in this study.

Variable	Sensor Model and Manufacturer	Position (m)—Sites
Wind speed components	3D sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, Utah)	10—PRS 3—CAS
H2O concentration	Open gas analyzer (LI-7500, LI-COR Inc., Lincoln, Nebraska)	10—PRS 3—CAS
Air temperature and relative humidity	Thermohygrometer (HMP45C-L; Campbell Scientific Inc., Logan, Utah)	10—PRS
Precipitation	Pluviometer (TB4 Rain Gauge; Campbell Scientific Inc., Logan, Utah)	10—PRS 6—CAS
Net radiation	Net Radiometer (NR LITE, Kipp&Zonen, Delft, The Netherlands)	8—PRS 3.4—CAS
Global Radiation	Pyranometer (CMB6, Kipp&Zonen, Delft, The Netherlands)	8—PRS
Soil temperature	Thermometer (STP01 Hukseflux Thermal Sensors B.V., Delft, The Netherlands)	−0.02; 0.05—PRS
Ground heat flux	Heat flux (HFP01SC-L; Hukseflux Thermal Sensors B.V., Delft, The Netherlands)	−0.07—PRS
Air temperature	3D sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, Utah)	3—CAS
Relative humidity	Combined data from CSAT3 and LI-7500	10—PRS

. Further details are described in References [37,39,40,41].

We filled the missing values of air temperature and global radiation. At PRS, data gaps were shorter than three successive hours and a simple linear interpolation was used to fill missing data. CAS site gap-filling was performed using a linear regression scheme based on meteorological measurements at a nearby station (World Meteorological Organization, WMO 86977), from the INMET. The station is located approximately 80 km from the CAS site. Climatological precipitation for both sites was obtained from INMET estimates (http://www.inmet.gov.br/).

In addition to meteorological and surface flux measurements, leaf area index (LAI) was also measured at the CAS site. LAI was measured at approximately every two weeks, for the 2010–2011 and 2011–2012 growing season using LI-3000A (LI-COR Inc., Lincoln, Nebraska) and 2014–2015 growing season using LAI-2200C Plant Canopy Analyzer (LI-COR Inc., Lincoln, Nebraska). All instruments were factory calibrated and measurements were conducted at 9:00 a.m. LAI measurements were used to evaluate and validate MOD15 estimation. The MOD15 product was used as forcing data to run the MOD16 algorithm. MOD15 validation results are presented in Figure S1.

2.3. ET Experimental Data Processing

The experimental latent heat flux (LE) was computed from EC measurements. The EC technique estimated surface fluxes using the covariance between vertical wind speed measurements and H₂O (or air temperature) concentration in the atmosphere to compute LE fluxes (or sensible heat flux (H)). The EddyPro Advanced software (version 5.1, LI-COR) was used to process the flux data over half-hourly time scales. Further details of the configurations used are described in References [39,42].

The footprint of the EC was calculated using the software EddyPro Advanced according to the methodology proposed by Kljun et al. [43]. The fetch was approximately 470 m for PRS and 200 m for CAS in all directions. The flux towers were located in flat terrains with homogeneous surface for both sites, however the footprint analyses were not necessary to filter any flux data. LE and H were filtered in periods that presented physically inconsistent data [44], i.e., LE < −50 W m⁻² and LE > 650 W m⁻², H < −70 W m⁻² and H > 300 W m⁻² were removed from the analysis. During the night period, LE > 100 W m⁻² was also removed. During precipitation periods, LE and H values were excluded. After data filtering, which also included malfunctioning periods, total LE data gaps accounted for 30.1% and 29.15%, for PRS and CAS sites, respectively. These percentages of gaps were similar to those found in others studies conducted in rice paddies [37,45,46]. Gap filling for H and LE missing data was performed through the REddyProc package available in RStudio software, using the following entry data: LE, H, air temperature (Temp), global radiation ($R_g$). This package was developed using methods proposed by Reichstein et al. [47], further details are also described in Wutzler et al. [48].

Several studies show the imbalance between turbulent fluxes (H + LE) and available energy ($R_n$ − G; where $R_n$ is defined as net radiation and G is the soil heat flux) [35,37,42,49,50]. Many authors suggest the necessity of adjustment in the values of H and LE [51,52,53,54]. The energy balance closure can be forced using the Bowen ratio ($\beta$) technique, defined as $\beta =H/LE$. The residual available energy (RAE = $R_n-H-LE-G$) is then distributed between H and LE to close the energy balance, using the following equations (according to Twine et al. [52]):

$$H=\frac{\beta \left(R_n-G\right)}{1+ \beta },$$

(1)

$$LE=\frac{R_n-G}{1+ \beta }.$$

(2)

For the PRS site, G was estimated through the sum of ground heat flux (Fg) and heat storage in the soil and water layer (ΔG). In this work, ΔG was estimated according to Timm et al. [37] using:

$$\Delta G= \frac{c_s\Delta T_{5}d{z}_s}{dt}+\frac{c_a\Delta T_{2}d{z}_a}{dt}$$

(3)

where c_s is the volumetric heat capacity of soil (3.03 × 10⁶ J m⁻³ °C⁻¹) and c_a is the volumetric heat capacity of water (4.186 × 10⁶ J m⁻³ °C⁻¹), ∆T₅ and ∆T₂ are the soil temperature variation at depth of 0.05 m and 0.02 m, respectively, dt is the time step; dz_s is the soil layer above the soil temperature sensor and dz_a is the water layer. The heat storage in the water layer was calculated only for the flooded period. The values of Fg and ΔG represent 1.15% and 0.19% of the total $R_n$, respectively, for the PRS site. For the CAS site, Fg and ΔG were not available. However, we estimated G for CAS, using a linear regression between G and R_n based on data from the PRS site, according to the following equation:

$$G=0.19\times R_n-17.$$

(4)

To evaluate the energy partitioning between LE and H during flooded and fallow (non-flooded) periods, we also computed the evaporative fraction (Λ), defined as the amount of available energy used for evapotranspiration [55], according to the equation:

$$\mathsf{\Lambda }=\frac{LE}{R_n-G}.$$

(5)

2.4. MOD16 Global Evapotranspiration Algorithm

2.4.1. MOD16 Algorithm Structure

Several techniques using remote images have been developed to estimate ET rate over large areas. MOD16 algorithm was developed to estimate land surface ET [23,56,57], based on MODIS remotely sensed compositions and MERRA meteorological reanalysis data. This methodology was adapted from Cleugh et al. [33] using a Penman–Monteith approach to estimate LE. ET is given by the sum of soil evaporation ($\lambda$E_s), canopy evaporation ($\lambda$E_c) and canopy transpiration ($\lambda$T_c), according to the equation:

$$\mathrm{ET}=\lambda E_S+\lambda E_C+\lambda T_C.$$

(6)

Estimation of $\lambda$E_s, $\lambda$E_c and $\lambda$T_c is given by Equations (7)–(9), respectively [23]:

$$\lambda E_S=f_w\frac{\Delta \times A_s+\frac{\left(1-f_c\right)\times {\rho }_a\times C_P\times \left(e_s-e_a\right)}{r_a^s}}{\Delta +\gamma \times \frac{r_s^s}{r_a^s}}+R{H}^{\frac{\left(e_s-e_a\right)}{{\beta }_{sm}}}\left(1-f_w\right)\frac{\Delta \times A_s+\frac{\left(1-f_c\right)\times {\rho }_a\times C_P\times \left(e_s-e_a\right)}{r_a^s}}{\Delta +\gamma \times \frac{r_s^s}{r_a^s}},$$

(7)

$$\lambda E_c=f_w\frac{\Delta \times A_c+\frac{f_c\times {\rho }_a\times C_P\times \left(e_s-e_a\right)}{r_a^{wc}}}{\Delta +\gamma \times \frac{r_s^{wc}}{r_a^{wc}}},$$

(8)

$$\lambda T_c=\left(1-f_w\right)\frac{\Delta \times A_c+\frac{f_c\times {\rho }_a\times C_P\times \left(e_s-e_a\right)}{r_a^t}}{\Delta +\gamma \times \left(1+\frac{r_s^t}{r_a^t}\right)},$$

(9)

where $\Delta$ (Pa K⁻¹) is the gradient of the saturation vapor pressure–temperature; f_c is the vegetation cover fraction; $A_s$ and $A_c$ (W m⁻²) are the available energy to the soil and canopy, respectively (as $A_c=\left(f_c\times R_n\right)$ and $A_s=\left(1-f_c\right)\times \left(R_n-G\right)$); $C_P$ (J Kg⁻¹ K⁻¹) is the specific heat of the air, ($e_s-e_a)$ is the vapor pressure deficit (VPD, in Pa); $\gamma$ (Pa K⁻¹) is the psychrometric constant; ${\beta }_{sm}$ is a parameter related to the soil moisture constraint (${\beta }_{sm}$ = 200 [23,58]); $r_s^s$ and $r_a^s$ (s m⁻¹) are the surface and aerodynamic resistance for the soil surface; $r_s^{wc}$ and $r_a^{wc}$ (s m⁻¹) are the surface and aerodynamic resistance for the wet canopy evaporation and $r_s^t$ and $r_a^t$ (s m⁻¹) are the surface and aerodynamic resistance for the canopy transpiration; f_w is wet surface fraction, calculated as a function of e_a (when relative humidity (RH, in %) is higher than 70%, with $f_w=R{H}^4$ [23]), representing the fraction of canopy and soil covered by water.

The MOD16 final products provided averages of 8-day, monthly and annual potential and actual LE and ET, with spatial resolution of 1 km. MOD16 ET estimates are freely available at http://www.ntsg.umt.edu/project/mod16 and at https://search.earthdata.nasa.gov/.

2.4.2. MODIS Remotely-Sensed Inputs

MODIS 8-day and 16-day compositions, 500 m and 1 km spatial resolution, were used as MOD16 inputs: MOD15A2 LAI and fraction of photosynthetically active radiation (fPAR) [59] and combined Aqua & Terra Products: MCD43B2/B3 albedo [60] and MCD12Q1 land cover [61]. MOD15A2 [62] consists of an 8-day composite, where LAI and fPAR values are retrieved by an algorithm based on the radiative transfer process (RT). MCD43B2/B3 consists of a 16-day composition and it is used to estimate the reflected solar radiation. For both LAI/fPAR and albedo, unreliable pixels were replaced by the closest reliable pixels, using a linear interpolation technique proposed by Zhao et al. [63]. Although MOD16 relies on LAI/fPAR and albedo derived from MODIS compositions to estimate land surface ET, the evaluated algorithm is heavily dependent on land cover parameterization through the use of a Biome-Properties Look-up Table (BPLUT) related to a 14 classes land cover classification scheme (Land Cover Type 2 in the MCD12Q1 data). MCD12Q1 has an accuracy around 75%, but the range in class-specific accuracies is large (between 22.6% and 99.3% when considered errors of producer (commission) and user (omission) errors) [61]. MODIS subsets (MOD15, MCD43 and MCD12) are freely available at https://modis.ornl.gov/ and at https://search.earthdata.nasa.gov/.

2.4.3. Global Modeling and Assimilation Office (GMAO) MERRA Meteorological Reanalysis Inputs

As meteorological inputs, the current MOD16 version uses Global Modeling and Assimilation Office (GMAO) Modern-Era Retrospective analysis for Research and Applications (MERRA), with spatial resolution of 0.5° × 0.6°. Meteorological reanalysis inputs include daily total downward shortwave radiation (also defined as global radiation) (R_s), daily average air temperature (T_avg), daytime average air temperature (T_day), daily minimum air temperature (T_min) and specific humidity (q). To combine meteorological and remote sensing inputs, meteorological reanalysis data were interpolated from spatial resolution of 0.5° × 0.6° to 1 × 1 km to match the spatial resolution of MODIS pixels and to remove abrupt changes from the meteorological reanalysis data, using a fourth power cosine function proposed by Zhao et al. [63]. The VPD based on MERRA reanalysis data (VPD_MERRA) was computed as described by Mu et al. [56] and Running et al. [57].

2.5. Statistical Analysis

To evaluate and compare MOD16 ET estimates and EC measurements, we aggregated the 30-min and daily data to 8-day averages, corresponding to the MOD16 intervals. The following statistical metrics were used:

Pearson’s correlation

$$r=\frac{{{\displaystyle \sum }}^{\text{}}({\mathrm{EC}}_{\mathrm{i}}-\overline{\mathrm{EC}})\left({\mathrm{ET}}_{\mathrm{MODIS}}{}_{\mathrm{i}}-\overline{{\mathrm{ET}}_{\mathrm{MODIS}}}\right)}{\sqrt{{\left({{\displaystyle \sum }}^{\text{}}({\mathrm{EC}}_{\mathrm{i}}-\overline{\mathrm{EC}}\right)}^2)(}{{\displaystyle \sum }}^{\text{}}{\left({\mathrm{ET}}_{\mathrm{MODIS}}{}_{\mathrm{i}}-\overline{{\mathrm{ET}}_{\mathrm{MODIS}}}\right)}^2)},$$

(10)

root mean square error

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}^{\text{}}{\left(\mathrm{EC}-{\mathrm{ET}}_{\mathrm{MODIS}}\right)}^2}{\mathrm{n}}},$$

(11)

Bias

$$\mathrm{bias}=\frac{{{\displaystyle \sum }}^{\text{}}\left({\mathrm{ET}}_{\mathrm{MODIS}}-\mathrm{EC}\right)}{\mathrm{n}},$$

(12)

Percent bias

$$\mathrm{PBIAS}=\frac{\mathrm{bias}}{\left(\frac{1}{\mathrm{n}}\right){{\displaystyle \sum }}^{\text{}}\mathrm{EC}}\times 100,$$

(13)

where EC is the ET measurement obtained by the eddy covariance technique, ET_MODIS is the MOD16 ET estimation and n is the number of samples. For Pearson’s correlation, we considered r ≤ 0.4 a low correlation, 0.4 < r ≤ 0.7 a moderate correlation and r > 0.7 a high correlation.

3. Results and Discussion

3.1. Energy Balance Analysis

Since the EC technique provides a direct measurement of LE and there is a well-documented missing energy or an energy balance closure “error”, with the turbulent fluxes (H + LE) being typically less than the 10%–30% of the available energy ($R_n-G$) [51,64], we evaluated the average energy balance components ($R_n=H+LE+F_g+\Delta G$, where $G=F_g+\Delta G$) during fallow (harvest to the subsequent planting) and flooded seasons (planting to harvest) (Figure 2). At the PRS site we used measured R_n, H, LE, F_g and calculated ∆G (Equation (3)), whilst at the CAS site we used measured R_n, H and LE and estimated G (according to Equation (4)). During both flooded and fallow seasons, the energy balance closure yielded an r of 0.95 at the PRS site. On the other hand, at the CAS site, the energy balance closure yielded an r of 0.91 during the fallow season and 0.93 during the flooded season. Ikawa et al. [65] find an r value for the interannual energy balance closure between 0.85 and 0.94 over 13 years of EC measurements in irrigated crops. According to Timm et al. [37] when ∆G is not included in the energy balance analysis, lower r coefficients are typically found, suggesting that the lower energy balance closure is due to uncertainties in the estimate of heat storage capacity, especially during the irrigation season and changes in water temperature measured in a single depth [66].

The average daily values of the energy balance components (without energy balance closure) during fallow and flooded seasons are presented in

Table 2. Daily averages of the energy balance components (without energy balance closure) during fallow and flooded seasons at the PRS and CAS sites.

	PRS		CAS
	Fallow	Rice	Fallow	Rice
Rn (W m−2)	67.62	146.16	59.84	146.08
LE (W m−2)	53.44	91.01	50.14	118.97
H (W m−2)	10.01	10.08	10.79	2.96
G (W m−2)	−0.36	4.56	−5.63	10.75
RAE (W m−2)	4.53	34.53	4.54	13.40
Λ ± stdev	0.87 ± 0.22	0.92 ± 0.13	0.86 ± 0.21	0.98 ± 0.11

. We found R_n values of 146.16 and 146.08 W m⁻² for the PRS and CAS sites during flooded season, with LE values corresponding to 62% and 81% of R_n, respectively. During fallow season, we found R_n values of 67.62 and 59.84 W m⁻² for the PRS and CAS sites, with LE values close to 79% and 84% of R_n, respectively. H represented 6.9% to 2.0% of R_n during irrigation season and 14.8% to 18.0% of R_n during fallow season, for PRS and CAS, respectively. The ground heat flux ($F_g+\Delta G$) was negative for both sites during fallow season (almost zero at PRS), representing −0.5% and −9.4% of R_n at PRS and CAS, respectively, whilst during irrigation season the ground heat flux was positive, with 3.1% and 7.3% of R_n at PRS and CAS, respectively. In general, PRS and CAS present similar behavior on the energy balance components. A detailed discussion about the energy partitioning over rice paddy, using the PRS dataset, is described in Timm et al. [37] and using the CAS dataset, is described in Diaz et al. [39].

Several studies investigate the energy partitioning and the energy balance closure problem [37,50,51,67], despite considerable uncertainties in irrigated crops, such as rice [68,69,70]. RAE represented approximately 7% of R_n during fallow season (for both sites) and 23.6% and 9.2% of R_n during flooded season for PRS and CAS, respectively. Baldocchi et al. [35] also find a difference up to 15% of RAE for interannual measurements over flooded rice.

After the energy balance closing at daily timescale, we estimated daily Λ, that presented high values, with 0.92 and 0.98 at PRS and CAS sites, respectively, during irrigation season. It indicated that almost all energy was converted to LE, reflecting high water availability and rice canopies fully developed. During fallow season, Λ decreased to 0.87 and 0.86 at the PRS and CAS sites, respectively, with a higher standard-deviation (approximately 0.22, compared to 0.13 during the irrigation season). Kar and Kumar [71] find Λ values ranging from 0.6 to 1.0 for irrigated rice crops, indicating that as rice crop emerges, and until reaching full development, the main energy partitioning driver is LAI, causing a decrease in soil evaporation and an increase in plant transpiration. Bouman et al. [72] find an increase of Λ during crop-growing season and a decreasing trend during leaf senescence and crop maturation (before harvesting). The seasonal variation of daily Λ at PRS and CAS sites is shown in Figure 3, where it can be seen that higher values were concentrated in rice growing season. In fallow season, daily Λ values increased from harvest until new growing season. This fact can be attributed to the seasonality of climatic conditions.

3.2. Comparison between MOD16 Evapotranspiration Estimates and Eddy Covariance Measurements

In this section, we evaluate the accuracy of MOD16 ET estimates against EC measurements (here we included the RAE partitioned between H and LE measurements) at the two flux towers sites over rice paddy areas. The performance of MOD16 (ET_MODIS) was assessed at the pixel within the same location of the flux tower, with results presented in

Table 3. Statistical analysis between MOD16 evapotranspiration (ET) estimates against eddy covariance (EC) measurements in two rice paddy sites located in Southern Brazil.

Site	Period	r	RMSE (mm 8-Day−1)	PBIAS (%)	Bias (mm 8-Day−1)
PRS	2003–2004	0.65	13.40	−38.70	−8.77
CAS	2011–2014	0.50	16.35	−33.70	−8.26
CAS	2011	0.66	13.55	−20.70	−4.62
CAS	2012	0.61	14.22	−43.00	−9.83
CAS	2013	0.41	19.12	−34.80	−8.99
CAS	2014	0.38	17.83	−35.60	−9.58
PRS	Fallow	0.34	10.75	−44.20	−6.88
PRS	Rice	0.41	17.17	−34.20	−12.14
CAS	Fallow	0.44	9.49	−43.30	−6.93
CAS	Rice	−0.05	23.04	−27.6	−10.22

and Figure 4a,b. When computing the correlation between EC measurements and MOD16 ET estimates, we found r values of 0.65 and 0.50, respectively, at the PRS and CAS sites. At the CAS site, the highest correlations and the lowest errors were observed during 2011 and 2012 (with r close to 0.63 for both years, RMSE between 13.55 and 14.22 mm 8-day⁻¹). On the other hand, during 2013 and 2014, we found lower correlations and higher errors, with RMSE between 17.83 and 19.12 mm 8-day⁻¹. Bias and PBIAS showed negative values, resulting in an overall underestimation of MOD16 ET over flooded rice paddy, for all evaluated periods at both sites. At the PRS site, MOD16 ET was underestimated by 8.77 8-day⁻¹ resulting in PBIAS value of −38.7. Similar statistics were found at the CAS site for the 2011–2014 period and for each growing season.

The evaluation of MOD16 ET was also performed separately for fallow and flooded (rice growing) seasons (Table 3), with no correlation observed in rice growing season at the CAS site (r ≈ 0) and low correlations for the PRS site. Moderate correlation was found during the fallow season for both sites. At both sites, RMSE and bias increased from fallow to flooded season, on the other hand PBIAS decreased. During the beginning of the crop growing season, the pixels may present a spectral mixture related to land cover (vegetation) and large amount of water in soil background, compromising the accuracy of ET estimates obtained from remote sensing. This sub-pixel mixture may lead to inaccuracies in fPAR and albedo estimates, used as input parameters to MOD16 ET algorithm. This problem is reported by Kim et al. [26] when validating MOD16 estimates for different land cover classes in Asia.

The annual accumulated EC measurements and MOD16 ET estimates at the PRS and CAS sites are presented in

Table 4. Annual MOD16 ET, tower-driven MOD16 ET estimates, and EC measurements at the PRS and CAS sites.

Site	Year	ETMODIS (mm Year−1)	ETMODIS-Tower (mm Year−1)	EC (mm Year−1)
PRS	2003–2004	665.9	592.8	1079.6
CAS	2011	814.5	617.6	1027.1
CAS	2012	599.4	529.6	1051.7
CAS	2013	776.0	642.6	1189.5
CAS	2014	796.3	764.5	1237.2

. The MOD16 ET underestimation values ranged from 212.6 mm year⁻¹ to 452.3 mm year^-1, representing 20.7% and 43% of the EC. These results and the statistical analyses (Table 3) indicated a poor relationship between MOD16 ET estimates and EC measurements over rice paddy in Southern Brazil. The validation of MODIS ET estimates over the rice paddy sites, showed several discrepancies when compared to the EC measurements, mainly underestimating ET. Kim et al. [26] validate MOD16 in two rice paddy sites, reporting overestimation in the site located in South Korea, and underestimation in the site located in Japan, with RMSE of 6.12 and 10.69 mm 8-day⁻¹, respectively. Khan et al. [73] assess MOD16 algorithm performance in nine EC sites in Asia, finding RMSE between 6.71 and 10.16 mm 8-day⁻¹ over three irrigated rice sites. Overall, the r and RMSE values for the rice paddy in Southern Brazil were, in magnitude, similar or larger than those found in the reported studies. However, in this region, rice production occurs once a year, a timeline differing from the practice in Asia, where twice a year the rice is cultivated.

Several authors also evaluate MOD16 around the globe. Hu et al. [18] validated MOD16 ET product over Europe, finding a mean RMSE value of 0.76 mm day⁻¹ (≈ 6 mm 8-day⁻¹) and bias ranging from −0.93 to 1.11 mm day⁻¹ with an average value of 2 mm 8-day⁻¹. For arid regions, Ramoelo et al. [15] validate MOD16 ET in South Africa, yielding r values between 0.51 to 0.93 and RMSE lower than 9 mm 8-day⁻¹ when compared to EC measurements. Meanwhile, Velpuri et al. [29] find monthly r values ranging from 0.40 to 0.77 and RMSE of 29.5 mm month⁻¹ (representing an error of approximately 8 mm 8-day⁻¹). Therefore, the MOD16 ET for dry regions is better related to EC than in wet regions, represented here for flooded areas used for rice paddy.

Despite uncertainties related to the algorithm structure and inputs, there was also uncertainty related to spatial resolution. A single MODIS pixel can represent multiple land cover conditions, whilst reanalysis data with coarse spatial resolution was not able to represent micrometeorological processes at smaller scales [30,74]. The mismatch between eddy covariance measurements and coarse spatial resolution from remote sensing data could introduce large errors, reducing ET up to 15% at large scale and 50% at pixel scale when compared to higher resolution estimates [75] due to losses on surface heterogeneity to compute surface fluxes. Therefore, lower spatial resolution can result in an increased likelihood of discrepancies between estimated and measured fluxes [76].

3.3. Seasonality and Climatic Drivers of Evapotranspiration

To evaluate MOD16 ET control drivers, we evaluated the cross-correlation between EC measurements, MOD16 ET estimates, measured environmental variables (Temp, $R_g$ and VPD) and MOD16 ET algorithms inputs (LAI, fPAR, and albedo from MODIS datasets and air temperature (Temp_MERRA), global radiation ($R{g}_{MERRA}$) and vapor pressure deficit (VPD_MERRA) derived from MERRA reanalysis). The full cross-correlation is presented in

Table 5. Cross Correlation coefficients between EC measurements, MOD16 ET products and atmospheric variables for two rice paddy sites in Southern Brazil.

	ETMODIS	EC	Temp	${\mathit{R}}_{\mathit{g}}$	${\mathit{R}}_{\mathit{n}}$	VPD	LAI	fPAR	Albedo	TempMERRA	RgMERRA	VPDMERRA
ETMODIS	1	0.65	0.57	0.47	0.54	0.44	0.92	0.82	0.48	0.55	0.47	0.43
EC	0.5	1	0.77	0.91	0.95	0.81	0.69	0.50	0.76	0.77	0.89	0.75
Temp	0.52	0.82	1	0.68	0.78	0.82	0.56	0.48	0.77	0.99	0.72	0.77
${\mathit{R}}_{\mathit{g}}$	0.42	0.91	0.78	1	0.98	0.86	0.55	0.35	0.83	0.70	0.96	0.82
${\mathit{R}}_{\mathit{n}}$	0.43	0.96	0.83	0.95	1	0.87	0.59	0.40	0.82	0.79	0.95	0.83
VPD	0.34	0.65	0.72	0.74	0.67	1	0.54	0.40	0.85	0.84	0.86	0.92
LAI	0.88	0.53	0.55	0.50	0.50	0.38	1	0.92	0.52	0.57	0.54	0.56
fPAR	0.80	0.44	0.51	0.44	0.43	0.30	0.90	1	0.36	0.47	0.32	0.42
Albedo	0.71	0.29	0.43	0.35	0.31	0.23	0.72	0.86	1	0.81	0.88	0.84
TempMERRA	0.52	0.82	0.98	0.79	0.83	0.67	0.58	0.56	0.48	1	0.76	0.83
RgMERRA	0.41	0.90	0.79	0.94	0.82	0.73	0.49	0.40	0.29	0.80	1	0.86
VPDMERRA	0.24	0.80	0.77	0.83	0.82	0.61	0.49	0.47	0.34	0.84	0.84	1

* white cells present data for PRS, whilst light grey cells present data for the CAS site and dark grey are the 1:1 cells. Units of the variables: ET_MODIS and EC (mm 8-day⁻¹); temperature (Temp) and temperature from the Modern-Era Retrospective analysis for Research and Applications (MERRA) (Temp_MERRA) (°C); global radiation (R_g) and R_gMERRA (W m⁻²); vapor pressure deficit (VPD) and VPD_MERRA (hPa); and leaf area index (LAI) (m² m⁻²).

. Measured environmental variables ($R_g$, VPD and Temp) were well correlated with EC measurements at both sites, whilst remote sensing variables (LAI, albedo and fPAR) were less correlated with EC measurements. For instance, when compared to EC measurements, $R_g$ and $R_n$ yielded correlations higher than 0.91, Temp yielded correlations higher than 0.77 and VPD yielded correlations higher than 0.65. In contrast, remote sensing estimates showed moderate to low correlations, with LAI and albedo yielding correlations between 0.37 and 0.74, respectively. Similar results were found by Fang et al. [24] for different land cover types. Our results revealed that the main drivers of measured ET were radiation, followed by Temp and VPD, with LAI showing a secondary role in measured ET. On the other hand, for MOD16 ET estimates, the main drivers were LAI and fPAR, yielding correlations higher than 0.88 for LAI and 0.80 for fPAR, with meteorological reanalysis data playing a secondary role in MOD16 algorithm.

To assess MERRA reanalysis accuracy, we compared surface meteorological data used as MOD16 inputs against observed meteorological data. A full detail of MERRA reanalysis validation is presented in Tables S1 and S2. All three MERRA reanalysis (Temp, $R_g$ and VPD) variables were overestimated at daily time steps. Temp and $R_g$ presented high accuracy, with PBIAS lower than 6.7% and 16.2% and r values higher than 0.95 and 0.87, respectively. On the other hand, VPD showed higher PBIAS (ranging from 69% to 87%) and moderate r (with values lower than 0.78), possibly related to the spatial resolution of the MERRA gridded data, that do not adequately capture subgrid meteorological scale. The VPD values, in general, are lower in wet than in dry regions [77]. Moreover, these statistical differences can be related to the methodology to compute VPD in MOD16 algorithm using a function of terrain elevation only [56,57].

The mean seasonal pattern of EC measurements, for MOD16 estimation and MERRA environmental variables at the PRS and CAS sites, are presented in Figure 5 and Figure 6, respectively. Overall, it is shown that MOD16 ET estimates followed the seasonal patterns of LAI and fPAR, whilst ET measurements presented a seasonal pattern similar to radiation components and air temperature. The greatest differences between EC and MOD16 ET were found during the spring fallow season (from September to November), with EC measurements around 20 mm 8-day⁻¹ greater than MOD16 ET estimates. In this period, Λ values increased (Figure 3), but LAI values were still low. During the winter fallow season (from May to August), MOD16 ET products and EC measurements showed similar behaviors, but in this period, Λ and LAI values are low. In contrast, during the maximum plant development (maximum LAI), MOD16 ET estimates were higher than EC measurements at the CAS site, reaching differences near to 10 mm 8-day⁻¹, whilst at PRS EC measurements and MOD16 ET estimates showed similar values. Overall, our results were similar to those reported by Biggs et al. [19], with larger differences in the beginning and in the middle of growing season, coincident with low LAI values and high water availability (during flooding), suggesting that MOD16 underestimates soil evaporation when soil wetness is high.

3.4. Sensitivity of MOD16 Algorithm on the Accuracy of Meteorological Data

To investigate the sensitivity of meteorological inputs on MOD16 ET estimates accuracy, we ran the MOD16 algorithm with measured meteorological forcing data and land cover parameterization as cropland [23,56,57] for the PRS and CAS sites (called “tower-driven” results), using the pixel at the same location of the flux tower to extract MODIS remotely-sensed inputs, as reported in Khan et al. [73] and Yao et al. [78]. The annual integrated values are showed in Table 4. The values were smaller than MOD16 ET at approximately 13%, on average. Therefore, using tower-driven data did not improve the results of MOD16 ET over flooded area.

The 8-day tower-driven estimates were then compared to EC ET measurements, as shown in Figure 7a,b. Tower-driven ET estimates did not indicate a high accuracy improvement, with RMSE increasing from 13.40 to 13.42 mm 8-day⁻¹ at PRS and decreasing from 16.35 to 15.14 mm 8-day⁻¹ at CAS, when compared to MOD16 ET estimations using MERRA reanalysis. Both estimates (using MERRA and observed meteorology) yielded similar r values (increased from 0.65 to 0.76 at PRS and from 0.50 to 0.62 at CAS) and similar bias and PBIAS (for both sites, we found BIAS increasing at ≈2 mm 8-day⁻¹ and PBIAS increasing ≈7% for PRS and ≈10% for CAS). Other studies also suggest that tower-driven ET estimates indicate a low sensitivity of MOD16 algorithm to meteorological inputs, suggesting that biases are related to algorithm limitations, such as model structure to compute surface (r_s) and/or aerodynamic (r_a) resistance and canopy conductance (c_c) [31,79] or to land cover parameterization [80]. In this case, further assessments are needed to improve global ET algorithms for flooded rice paddy or water-saturated soils.

The ET partitioning between λEs, λEc and λTc is presented in Figure 7c,d. In the spring fallow season (from September to November), months of preparation and sowing, LAI is minimal (i.e., without vegetation influence on surface fluxes). Moreover, with a water layer over the soil surface, mainly at the end of this period and beginning of the growing season, we expected no contribution to the final sum of ET from transpiration (λTc) and canopy evaporation (λEc), soil evaporation (λEs) being the key parameter in this period. However, this was not observed in the estimation of partition by MOD16 ET, that represented a larger contribution of λTc and λEc. In this period, we found the largest differences between EC measurements and MOD16 ET. Therefore, the equation of λEs needs adjustment in its parameterization Equation (7) to accurately describe the influence of soil moisture in the ET process, as discussed in Section 3.3. He et al. [81] describes that MOD16 algorithm does not directly represent soil moisture water supply constraints to ET, increasing significantly errors in the model. Recently, Purdy et al. [82] suggest soil moisture incorporation into of λEs improves ET estimates thus reducing uncertainties and provide a rich dataset to evaluate land surface and climate models.

Diaz et al. [39], showed the transpiration in CAS in fallow season (April to November) almost vanish, while evaporation dominates the ET. This behavior was also represented by MOD16 algorithm (Figure 7c,d). However, the processes of λEc and λTc were overestimated when LAI was close to zero, suggesting more studies in these parametrizations (Equations (8) and (9), respectively) are necessary. As expected, during the growing season (October to March), the transpiration is the dominant process for estimating ET, followed by canopy evaporation. Hossen et al. [83] find that about 70% of seasonal ET comes from transpiration in the growing season, whilst Diaz et al. [39], analyzing ET partition for the CAS site, find that around 60% of ET is due to transpiration process. Here we found 80% for PRS and 50% for CAS. However, the contribution of the water layer in the λEs during growing season should probably be greater than the estimate of λEs by the model (8% for PRS and 14% for CAS). An alternative to decrease the error in MOD16 ET estimates over flooded rice paddy is to improve the flooded fraction remotely estimated [84].

4. Conclusions

The main goal of this study was to conduct a comprehensive assessment of ET over two flooded rice paddy sites located in Southern Brazil using measured and remote sensing data, to address the main control drivers of the ET process. ET in irrigated crops is the largest consumer of the available energy. In this irrigation system, the fallow season (non-flooded) also showed a large importance for the available energy partitioning. During the irrigation season, almost all available energy was converted to LE. Since both flux tower sites are located in lowland areas (floodplains) with high water availability, the environmental measured variables indicated that the main drivers of measured ET were radiation ($R_g$ and $R_n$), followed by Temp and VPD.

MOD16 ET estimations over flooded rice paddy at PRS and CAS were underestimated when compared to EC measurements, with largest differences during the rice growing season. MOD16 algorithm main control drivers were LAI and fPAR, with meteorological reanalysis playing a secondary role. This result was supported by the sensitivity analysis of the algorithm. When forced with meteorological measurements, tower-driven ET estimates did not improve accuracy as expected, yielding similar statistics when compared to the meteorological reanalysis ET estimates.

Discrepancies between MOD16 ET estimates and EC measurements may be related to several factors, including flux tower measurement errors, flux tower footprint vs. MERRA–MODIS spatial resolution and land cover parameterization of the Penman–Monteith equation, as well as algorithm limitations. The cause of underestimation observed in this study may be also associated to model structure to compute r_s and c_c and land cover parameterization. Since the algorithm showed a low sensitivity to VPD, possibly caused by the simple methodology used to its determination as a function of surface elevation only, further studies are necessary to improve resistance factors and ET estimates. Our findings support that vegetation-based models to estimate global ET underestimate soil evaporation and plant transpiration in saturated soils, possibly due to inaccurate land cover parameterization and resistance factors. Therefore, the MODIS algorithm to estimate land surface ET is not able to detect areas with high soil moisture, such as flooded agricultural areas, mainly because the MOD16 algorithm relies on vegetation indices as surface processes indicators and MERRA reanalysis data as an indicator of moisture stress.

Results obtained in this research can guide additional model refinements to increase the accuracy of ET estimates. Therefore, we would recommend a regional scale calibration of the MOD16 algorithm. Additionally, as an alternative to reduce the impacts of parameterization, land surface temperature or soil moisture could be included as an algorithm input, to improve ET estimates over irrigated crop. Improvements in remote-sensed ET models could bring global ET estimates closer to experimental measurements and, therefore, help to improve water management.