Modeling bulk surface resistance and evaluation of evapotranspiration using remote sensing and MATLAB

In developing countries, computation of Actual Evapotranspiration (AET) is challenging due to the lack of ground-based ﬂ ux measurement data. The estimation of AET is crucial for water resources management involving the allocation of water for different land cover (LULC) classes. The study ’ s novelty was mapping pixel-by-pixel spatial variations of bulk surface resistance and evaluating the derived actual evapotranspiration in a sub-humid tropical river basin where ﬂ ux tower data was lacking for validation. This study aimed to map bulk surface resistance and evaluate the estimated AET by global evapotranspiration data product (MOD16A2). MODIS data products, including Land Surface Re ﬂ ectance (LSR), Land Surface Temperature (LST) and Leaf Area Index (LAI) data, were used as input in the MATLAB for mapping pixel-wise variations to analyze the seasonal variations in bulk surface resistance (r s ) and actual evapotranspiration in pre-monsoon and post-monsoon seasons during the years 2019 and 2012. The years 2019 and 2012 were selected because 2019 experienced a relatively wet pre-monsoon and post-monsoon, whereas 2012 experienced the opposite conditions, which proved useful when interpreting variations that are in ﬂ uenced by wetness conditions. Overall the results indicated signi ﬁ cant variability in the r s and AET for different LULC classes. MOD16A2 AET was determined to be slightly higher than the LULC classes estimated AET. This study ’ s Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data products provided information on surface characteristics at a reasonable resolution. This permitted the identi ﬁ cation of differences in LULC classes and changes in surface characteristics by season and wetness conditions, which are helpful when estimating AET. words: actual evapotranspiration, bulk surface resistance, MATLAB, remote sensing HIGHLIGHT (cid:129) The study ’ s novelty was mapping pixel-by-pixel spatial variations of bulk surface resistance and evaluating the derived actual evapotranspiration in a sub-humid tropical river basin where ﬂ ux tower data was de ﬁ cient for validation.


INTRODUCTION
Characterizing the spatial variabilities of Actual Evapotranspiration from heterogeneous landscapes is essential in hydrology, climate variability, agriculture, irrigation, water resources engineering and management as environmental impact assessments (Droogers et al. 2010;Minacapilli et al. 2016;Khaniya et al. 2020). AET combines two separate processes where moisture is lost from the soil surface by direct evaporation and via vegetation by transpiration. Furthermore, AET is a significant water balance component, sometimes accounting for over 60 percent of the annual water balance. Therefore, enhanced methods are required for the accurate computation of AET in a field to improve the efficient use of water resources and productivity (Hobbins et al. 2001;Howell 2001;Allen et al. 2007;Taylor et al. 2013).
Modeling AET is complex and the direct measurement of AET in large basins can be challenging. The MODIS-based remote sensing technology, which can also provide estimates in ungauged areas where ground-based flux observations are unavailable, is one viable approach for computing regional AET. Numerous literature reviews have indicated that MODIS data products provide LULC characteristics information to estimate AET conveniently and cost-effectively. (e.g., Cleugh et al. 2007;Mu et al. 2007;Leuning et al. 2008;Zhang et al. 2009).
AET can be computed using the Penman-Monteith (PM) equation. Several researchers have proposed using remotely sensed data to estimate the daily AET using the PM approach at regional and global scales (Nagler et al. 2005;Mu et al. Water Supply Vol 00 No 0, 2 Uncorrected Proof 2007;Leuning et al. 2008;Zhang et al. 2008;Guerschman et al. 2009). However, bulk surface resistance is a complex parameter influenced by numerous variables; estimating it when using the PM method is stimulating. Bulk surface resistance modeling has been the subject of intense research to develop a PM model. The global evapotranspiration data product was developed by Mu et al. (2007) based on the Cleugh et al. (2007) surface conductance model.
A study of the literature on evaluating the spatial pixel-wise AET estimates of global evapotranspiration data product (MOD16A2) was conducted and the results are summarised below. Ruhoff et al. (2013) used data from flow towers to estimate AET, then validated using the MGB-IPH hydrological model and MOD16A2 AET data. Cherif et al. (2015) compared the results to MOD16A2 with a remotely sensed SEBAL model to improve the accuracy of AET estimates. Finally, Ke et al. (2017) utilized a machine learning model to compute AET-based spatial estimates and compared to MOD16A2, the study found that the Landsat data fusion provided high accuracy for estimated AET.
It is evident from the literature review that several studies exist and have modeled surface resistance and AET using different data inputs, methodologies, algorithms, satellites, spatial scales and validation techniques in other climatic conditions (Li et al. 2016;Dimitriadou & Nikolakopoulos 2021;Mobilia & Longobardi 2021). However, most of the algorithms use satellite data combined with ground-based meteorological and flux data. Since due to the lack of flux data for validation of the proposed model, the AET estimation was undertaken using the surface conductance model of Penman-Monteith integrated with Leuning et al. (2008) (PML), which requires satellite data with limited meteorological data inputs. Hence the objective of this study was to use satellite imagery with meteorological data inputs and MATLAB to map bulk surface resistance and evaluate the estimated AET by the MOD16A2.

STUDY AREA
On the Kaveri River's northern bank (Figure 1), the Hemavathi River is among the most important tributaries. It rises in the Western Ghats near Ballalara-Yanadurga in the Mudigere taluk of Chikmagalur district. The river basin is located between the North latitudes of 13°22 0 30″ to 12°35 0 15″ and the East longitudes of 75°31 0 30″ to 76°39 0 45″. The river joins the Cauvery River in the Krishnarajasagar Reservoir in Akkihebbal after a 245-kilometer. The Hemavathi River basin covers 5,427 square kilometers. From March through May, the pre-monsoon season begins. From June until October, the rainy season lasts. During the rainy season, extreme rainstorms are common. November through February are the post-monsoon months where the weather is quite cold. The investigation's catchment is hilly, with steep to moderate slopes. The slope is steeper at higher elevations and eventually decreases. Figure 1 shows the overall elevation of the basin, which ranges from 748 to 1,853 metres above sea level. Plantations and agriculture are the backbone of the basin's economy. The Western Ghat section to the west is exposed to the full effect of southwest monsoon winds. Temperatures rise gradually from January through April in most parts of the basin. As a result, winds rise over the Western Ghats, bringing rain to the region. The passing of depressions is associated with heavy to highly severe rainfall. From June through September, the wind blows from the southwest to the northeast, causing the southwest monsoon to storm. Except for the Palghat gap, through which the southwest winds blow, the Western Ghats form a continuous barrier. As a result, the southwest monsoon shifts its direction from North-East to South-West around the end of September. The maximum quantity of cloudiness persists during the receding monsoon season. The cloudiest month is July, while the minimum cloudiest month is March. The Western Ghats partially covers the study area. As a result, there is a wide range of vegetative cover in this western region.

MODIS satellite data products
In this study, the MODIS satellite data products mounted on the Terra platform were used. The data can be freely downloaded from the website. The products used in this study included the MOD11A2 product of land surface temperature (LST), the MOD09A1 product of land surface reflectance (LSR), the MOD15A2 product of leaf area index (LAI) and the MOD16A2 product for global terrestrial evapotranspiration. MODIS data products were chosen based on post-monsoon seasons (November) and the pre-monsoon (April) of years 2012 and 2019.
Rainfall (Table 1) data were collected from the Indian Meteorological Department. Between 2012 and 2019, yearly rainfall ranged from a high of 2,165 mm to a low of 1,042 mm, with a 1,530 mm average. Compared to the year 2019, which had 2,165 mm of rainfall and a lower Land Surface Temperature (LST), the year 2012 had 1,042 mm of rainfall, which was less than the normal annual rainfall and had a high Land Surface Temperature (LST). The climate in the catchment  Uncorrected Proof region is typical monsoon climate. According to rainfall records, 2012 received significantly less rainfall both in pre-monsoon (dry condition) and post-monsoon (wet condition) than the annual average rainfall, while 2019 received significantly more rainfall both in pre-monsoon (dry condition) and post-monsoon (wet condition) than the annual average. As a result, these two years were chosen with care to study seasonal changes in LULC classes due to dry and wet moisture levels in pre-monsoon and post-monsoon seasons.

Land Use and land cover map
LULC data for the years 2012 and 2019 was taken from the Bhuvan website (https://bhuvan.nrsc.gov.in). The data utilised in this study was derived from Linear Imaging Self Scanning Sensor (LISS) -III data collected by the Resourcesat-1 satellite at a scale of 1:250,000. There are 19 classes (2nd level) in this data; however, only 14 land cover classes were found. The derived LULC map of the study area is shown in Figure 2. The existence of plantations may be observed in the basin's higher elevations (towards the West). However, agricultural crops and plantations were the most significant classes in the subhumid tropical river basin. Figure 3 depicts the flow chart of steps involved in using the PM model proposed in this study to derive estimates of surface resistance and AET using inputs from MODIS imagery and meteorological data. First, Pre-processing and post-processing of MODIS data products was undertaken using the MODIS Reprojection Tool (MRT), followed by the application of ArcGIS. Finally, all the necessary MODIS images were post-processed into TIFF formats for final processing in MATLAB. The flowchart depicts the step-by-step procedure for estimating the variables, then describes the parameters in the following sections involved in computing surface resistance and AET using the PM approach.  (1) and (2). (1) where τ ¼ exp (Àk A LAI) and k A ¼ coefficient of extinction for energy available, LAI¼ Leaf Area Index, Q 50 ¼visible radiation flux at which stomatal conductance is half its maximum value (MJ/m 2 /day), Q h ¼ visible radiation flux density at the top of canopy (MJ/m 2 /day), k Q ¼ extinction coefficient for shortwave radiation, D 50 ¼ vapour pressure deficit when stomatal conductance is half its maximum value (kPa), D a ¼ vapour pressure deficit of air (kPa), g sx ¼ maximum stomatal conductance (m/s), G c ¼ bulk canopy conductance (m/s) and G i ¼ climatological conductance (m/s) i.e., G i ¼ γ (Rn-G) /(ρacpDa).

METHODOLOGY
As per the recommendations of Leuning et al. (2008) k A ¼ k Q ¼ 0.6, Q 50 ¼ 30 W/m 2 and D 50 ¼ 0.7 kPa values can be kept constant values across different LULC classes and optimization is required for only g sx and f. According to Table 1 of Kelliher et al. (1995) the values of g sx are fixed for the different LULC classes. The optimum value of f was determined by running the model algorithm where the values of f was varied from 0 to 1 in increments of 0.1 until the minimum RMSE was obtained between actual evapotranspiration estimated.

The penman-monteith (PM) approach
The latent heat flux (lET) (W/m 2 ) is the energy used to evaporate the water which is estimated using the Penman-Monteith (PM) Equation (3) combined with a bulk surface resistance (r s ) which is estimated from leaf area index (LAI) as the inverse of Gs. The PM equation is: where D ¼slope of saturation vapor pressure curve at air temperature (kPa/°C) which was calculated using Richards (1971), e s and e a ¼ saturation and actual vapor pressure (kPa) of the air which was computed using Allen et al. (2007), γ, ρ a and c p ¼ psychrometric constant (kPa/°C), mean air density (kg/m 3 ) and specific heat of air (MJ/kg/°C) which was computed using per Allen et al. (2007), R n ¼ net radiation (W/m 2 ) which was estimated using the procedures given by Allen et al. (2007), Bastiaanssen (2000) and Tasumi et al. (2003), G ¼ soil heat flux (W/m 2 ) which was computed using Bastiaanssen (2000). r s and r a ¼ the (bulk) surface and aerodynamic resistances (s/m) which were computed by combining the procedures given by Leuning et al. (2008). The terms D, e s , e a , ρ a , r a and γ were calculated using Cleugh et al. (2007).

Estimating actual evapotranspiration (AET)
The instantaneous evaporative fraction (EF) Equation (4), (EF ¼ Rn-G) is a ratio of latent heat flow values to accessible energy flux values. It's used to show how energy is partitioned in order to derive daily energy balance information from satellite data. The daily AET may be computed using EF and daily average net radiation data, according to Morse et al. (2000).
Due to the use of instantaneous satellite sensors to estimate all parameters, 24 hours AET is computed using Equation (5) assuming that EF remains constant due to latent heat of flow (Sugita & Brutsaert 1991).

RESULTS AND DISCUSSION
As discussed earlier, the years 2012 and 2019 were selected for analysis as representative dry and wet years. These results with respect to antecedent rainfall conditions proved useful in the subsequent interpretation of differences in LULC characteristics and the influence of wetness conditions. Furthermore, implementing equations using satellite imagery yielded pixelwise values of r s and AET across the Hemavathi River basin on pre-monsoon and post-monsoon days. Figure 4 maps the bulk surface resistance variations derived form MATLAB with low average values observed in post-monsoon compared to pre-monsoon for all the LULC classes due to the high moisture availibility with a high leaf area index. It can be observed that high r s values are likely where the LAI is low and has less moisture availability, limiting the AET process. Low r s values are observed under high wetness and high LAI conditions, enhancing the AET. It can also be observed that the minimum values of r s were lower for all LULC classes during the post-monsoon and pre-monsoon of 2019 (wet year) in comparison to the post-monsoon and pre-monsoon of 2012 (dry year). The r s values also varied across the LULC classes with the lowest values calculated for the Plantations and evergreen forest, which lies in the higher elevation towards the left side. Table 2 shows the performance of the pixelwise AET values (5,427 pixels ¼ 5,427 km 2 ) computed by PML model compared to MOD16A2. For the LULC classes, the AET estimated by MOD16A2 was slightly higher than that estimated using the PML method. with an R 2 ¼ 0.78 to 0.84 with a root mean square error (RMSE) ¼ 0.39 to 0.48 mm/day and a bias (BIAS) ¼ 0.13 to 0.18 mm/day for the post-monsoon and pre-monsoon seasons of years 2012 and 2019.

Uncorrected Proof
Others (Ruhoff et al. 2013;Cherif et al. 2015;Autovino et al. 2016;Ke et al. 2017) found that MOD16A2 overestimated the AET (RMSE ¼ 0.4-0.9 mm/day) when comparing spatial pixelwise AET estimations with MOD16A2 . Similarly, the AET estimated using MOD16A2 was found to be slightly higher than that estimated by the PML method in this investigation. Figure 5 maps the daily evapotranspiration estimated using the PML approach in MATLAB, generated using AET pixelwise values. The variations in the AET were high and appeared to be related to wetness conditions and also the LULC class. As expected, the average AET values were low during the post-monsoon and high during the pre-monsoon. The ranking of the estimated AET was highest in Waterbodies, followed by Plantations, evergreen forest and grassland. It was found to be low in the gullied and other wastelands LULC classes. The topography, LULC, LST and moisture conditions influence the spatial distribution of AET variations across the basin. In both pre-monsoon and post-monsoon in 2012 and 2019, the AET at higher elevations of the basin was consistently higher than at lower elevations (see Figures 1 and 5). Due to high LST at the higher elevation, the average AET values are higher in pre-monsoon than post-monsoon in 2012 and 2019. For most classes, the AET in the post-monsoon of 2019 was slightly lower than in the post-monsoon of 2012 due to significantly less moisture in 2012 associated with reduced rainfall during the monsoon season.

CONCLUSIONS
The MODIS global evapotranspiration dataset was assessed and evaluated in the Hemavathi River Basin using the surface conductance model and MODIS satellite data products. The focus was on assessing the patterns of bulk surface resistance and evaluating the derived AET across different LULC classes present in the river basin. It was found that agricultural crops and plantations are the major classes in the basin. The analysis was undertaken in pre-monsoon and post-monsoon in 2012 (dry year) and 2019 (a wet year). The surface resistance and AET were computed for each pixel in the basin using MODIS-based satellite imagery products. Very high surface resistance values were observed under low wetness and low LAI conditions. The variations in AET appeared to be related to wetness conditions, LST, and the LULC class. Waterbodies had the highest AET losses estimated by the PML and MOD16A2 methods, followed by Plantations, evergreen forest, and grassland. For all the LULC classes, the AET estimated by MOD16A2 was slightly greater than that estimated using the PML method. It was concluded from the assessment that MODIS satellite data products provide information on earth surface characteristics at a resolution that permits identification of not only differences in LULC classes but also on changes in these characteristics as a function of season, LST and moisture conditions. This work's developed methodology provided reasonably accurate regional/catchment scale surface resistance, AET estimations using satellite images and limited data inputs. Different satellite images with finer spatial resolution like AVHRR and LANDSAT can be used to compare the performance of MODIS imagery. This will be critical in water balance investigations and hydrological models calibration.