Spatial Disaggregation of Latent Heat Flux Using Contextual Models over India

Estimation of latent heat flux at the agricultural field scale is required for proper water management. The current generation thermal sensors except Landsat-8 provide data on the order of 1000 m. The aim of this study is to test three approaches based on contextual models using only remote sensing datasets for the disaggregation of latent heat flux over India. The first two approaches are, respectively, based on the estimation of the evaporative fraction (EF) and solar radiation ratio at coarser resolution and disaggregating them to yield the latent heat flux at a finer resolution. The third approach is based on disaggregation of the thermal data and estimating a finer resolution latent heat flux. The three approaches were tested using MODIS datasets and the validation was done using the Bowen Ratio energy balance observations at five sites across India. From the validation, it was observed that the first two approaches performed similarly and better than the third approach at all five sites. The third approach, based on the disaggregation of the thermal data, yielded larger errors. In addition to better performance, the second approach based on the disaggregation of solar radiation ratio was simpler and required lesser data processing than the other approaches. In addition, the first two approaches captured the spatial pattern of latent heat flux without introducing any artefacts in the final output.


Introduction
Latent heat flux (λE, W m −2 or evapotranspiration, ET, expressed in mm) is a widely required input for agricultural water management.Proper crop water monitoring and irrigation scheduling require the estimation of λE at the field scale and at regular time intervals.This need for repeated monitoring of λE over a large number of agricultural fields in a region have led to the development of different remote sensing (RS)-based models based on the surface energy balance equation (e.g., [1,2]).Furthermore, the RS-based models have been operationalized (e.g., [3]) over some regions of the globe to estimate the field scale λE using high spatial resolution data (~100 m) from satellites, like Landsat-5, Landsat-7, and Landsat-8.
The primary input to the surface energy balance models is the land surface temperature (T Rad , K) observed by thermal infrared (TIR) sensors.Among the currently-operational satellites, most of the TIR sensors, except the one aboard Landsat-8, have spatial resolution on the order of 1000 m or more.Landsat-8 has a nominal revisit cycle of 16 days and hence, the number of high-resolution T Rad images available from the single Landsat satellite is not sufficient for periodic monitoring of land surface fluxes λE [4].Due to this reason, sensors with frequent revisits, but coarse spatial resolution, like MODIS, are often used for repeated estimation of λE.The coarse spatial resolution of TIR data from sensors, such as MODIS, necessitated the development of models to disaggregate T Rad [5][6][7][8] and use the disaggregated T Rad for estimating λE with improved spatial resolution.In addition to this, studies have also been carried out to disaggregate evaporative fraction (EF, defined as the ratio of latent heat flux to the sum of sensible and latent heat fluxes) [9] or ET itself [10].
The need for such spatial disaggregation models is further increased over India, where agricultural fields with dimensions of one or two hundred of meters are more common.A single MODIS pixel at 1000 m resolution may include several tens of plots and, hence, it becomes highly difficult to monitor λE at the field scale.However, multiple studies have been carried out for the estimation of λE over India [11][12][13], and they all analyzed data at a coarse spatial resolution (≥ 1000 m).A previous study [8] tested the utility of disaggregated T Rad for the estimation of fine resolution λE over an area in south India, however, the disaggregated λE was not validated due to the lack of ground measurements.Another study [9] developed a model called disaggregation of evaporative fraction (DEFrac), which takes in coarse-resolution EF and disaggregates it to a finer resolution using the EF-NDVI relationship.The DEFrac model, though found to work with good accuracy over Indian sites, the disaggregated EF should be converted into λE by multiplying with net available energy at the surface ((R n − G), R n is the net radiation and G is the soil heat flux, W m −2 ) for use in various applications.The estimation of (R n − G) at fine spatial resolution require several other inputs that are either available at the target resolution (e.g., NDVI, surface albedo etc.) or that can be disaggregated from coarser resolution data (e.g., T Rad ).The complexities involved in estimating disaggregated (R n − G) can be avoided if EF can be replaced with some other ratio, which requires minimal data to be converted into λE.Previous studies, e.g., [14][15][16], suggested that the global solar radiation ratio (R g , defined as the ratio of latent heat flux to the incoming global solar radiation) can be utilized for the estimation of λE.The only variable required to convert R g into λE is the global incoming solar radiation at the surface (R sd ) which can be observed or modelled relatively easier than (R n − G).To the best knowledge of the authors, there are hardly any studies available in the literature that compare different approaches for estimating λE at finer resolutions over Indian ecosystems.Based on these, the objectives of this study are formulated as follows: 1.
Disaggregate satellite estimated EF and R g from 1000 m into 250 m and subsequently convert the disaggregated EF and R g into daytime integrated latent heat flux (λE day ) at 250 m resolution; and

2.
Compare the λE day estimated using the first two approaches with that estimated from the disaggregation of T Rad using the DisTrad thermal sharpening model [9,17].

Study Sites and In Situ Data
The three disaggregation approaches were tested at five sites across India.Each of the five sites is equipped with a Bowen Ratio energy balance (BREB) micrometeorological tower and the location of all the five sites are indicated in Figure 1.Site 1 is in a young pine forest amidst an urban town in Northern India.All other sites are in farmlands with cropping predominantly carried out in two seasons known as kharif (June-September, summer monsoon season) and rabi (November-January, winter cropping season).Additional details about the five tower sites can be obtained from [13,17].For the proper determination of the various context spaces (e.g., T Rad -NDVI) for the models used in this study, spatial grids of approximately 100 km × 100 km were delineated enclosing the tower sites.
Each of the BREB towers are 10 m in height.Air temperature and relative humidity are measured at 2 m, 4 m, and 6 m above ground level.Wind speed and wind direction are measured at 2.5 m, 5 m, and 10 m heights.The different components of surface radiation are measured using a Kipp and Zonen (Delft, The Netherlands) four-component net radiometer (3 m above the ground).Two soil heat flux plates are installed at 0.15 m and 0.3 m below the surface and atmospheric pressure is measured at 2 m above ground level using a transducer.Apart from these, the tower has a tipping bucket rain gauge (at 1 m above ground) for measuring rainfall and a shaded pyranometer (at 4 m above ground) for diffuse radiation measurements.Three soil thermometers and tensiometers to measure, respectively, the soil temperature and the soil suction head are installed at 0.05 m, 0.15, m and 0.3 m below the surface.All the measured variables are sampled every five minutes and the five-minute samples are averaged over half-hourly intervals.The half-hourly averages are stored in a data logger and transmitted to the server of the Space Applications Centre located at Ahmedabad in India.
above ground) for diffuse radiation measurements.Three soil thermometers and tensiometers to measure, respectively, the soil temperature and the soil suction head are installed at 0.05 m, 0.15, m and 0.3 m below the surface.All the measured variables are sampled every five minutes and the fiveminute samples are averaged over half-hourly intervals.The half-hourly averages are stored in a data logger and transmitted to the server of the Space Applications Centre located at Ahmedabad in India.

Satellite Data
The satellite data required for the study were obtained from the MODIS sensor onboard the Terra and Aqua satellites.All of the data products used are daily datasets acquired during the daytime overpass of the satellites.TRad, Land surface emissivity (εs) and satellite view time were obtained from the MOD11A1/MYD11A1 dataset at 1000 m spatial resolution.The seven-band surface reflectance at 500 m and solar zenith angle at satellite view time (1000 m) were obtained from the MOD09GA/MYD09GA product.In addition to this, for disaggregation purposes, the surface reflectance at red and near-infrared bands at 250 m were also obtained from the MOD09GQ/MYD09GQ data product.Furthermore, air temperature (Ta) data was obtained from the MODIS atmospheric profile level 2 product at a 5000 m pixel resolution.A total of 248 datasets were downloaded from the two satellites for all five sites together (59, 54, 35, 42, and 58 datasets for sites 1-5, respectively).

Brief Overview of the Disaggregation Approaches
In this study, λEday was estimated either as a product of EF and (Rn − G)day or as a product of Rg and (Rsd)day (the subscript 'day' refers to the daytime integrated value).To disaggregate λEday from MODIS thermal data at 1000 m to 250 m spatial resolution, three different approaches were adopted.In Approach 1, EF was estimated, first, at 1000 m resolution and disaggregated to 250 m using the DEFrac model [9].Simultaneously, the MODIS TRad data was disaggregated from 1000 m to 250 m using the DisTrad thermal sharpening model [17] and subsequently used in the estimation of (Rn − G)day at 250 m resolution.These disaggregated EF and (Rn − G)day were multiplied to get λEday at 250 m

Satellite Data
The satellite data required for the study were obtained from the MODIS sensor onboard the Terra and Aqua satellites.All of the data products used are daily datasets acquired during the daytime overpass of the satellites.T Rad , Land surface emissivity (ε s ) and satellite view time were obtained from the MOD11A1/MYD11A1 dataset at 1000 m spatial resolution.The seven-band surface reflectance at 500 m and solar zenith angle at satellite view time (1000 m) were obtained from the MOD09GA/MYD09GA product.In addition to this, for disaggregation purposes, the surface reflectance at red and near-infrared bands at 250 m were also obtained from the MOD09GQ/MYD09GQ data product.Furthermore, air temperature (T a ) data was obtained from the MODIS atmospheric profile level 2 product at a 5000 m pixel resolution.A total of 248 datasets were downloaded from the two satellites for all five sites together (59, 54, 35, 42, and 58 datasets for sites 1-5, respectively).

Brief Overview of the Disaggregation Approaches
In this study, λE day was estimated either as a product of EF and (R n − G) day or as a product of R g and (R sd ) day (the subscript 'day' refers to the daytime integrated value).To disaggregate λE day from MODIS thermal data at 1000 m to 250 m spatial resolution, three different approaches were adopted.In Approach 1, EF was estimated, first, at 1000 m resolution and disaggregated to 250 m using the DEFrac model [9].Simultaneously, the MODIS T Rad data was disaggregated from 1000 m to 250 m using the DisTrad thermal sharpening model [17] and subsequently used in the estimation of (R n − G) day at 250 m resolution.These disaggregated EF and (R n − G) day were multiplied to get λE day at 250 m resolution.Approach 3 was similar to Approach 1, in the estimation of λE day ; however, the EF was estimated directly at 250 m using the T Rad disaggregated to 250 m.In Approach 2, the EF estimated at 1000 m was converted into R g and this was disaggregated to 250 m using the DiSoRa disaggregation model (Section 3.3).This disaggregated R g was multiplied with (R sd ) day to get λE day at 250 m.The schematic of the three approaches are presented in Figure 2 and the various models used are explained in the following sections.All the models used in this work were forced with inputs that are obtainable only from RS datasets without any ground measurements.In Figure 2, FVC and NDWI indicate, respectively, the fraction vegetation cover and normalized difference water index (sometimes referred as normalized difference moisture index), which was used in the disaggregation of T Rad [17].
resolution.Approach 3 was similar to Approach 1, in the estimation of λEday; however, the EF was estimated directly at 250 m using the TRad disaggregated to 250 m.In Approach 2, the EF estimated at 1000 m was converted into Rg and this was disaggregated to 250 m using the DiSoRa disaggregation model (Section 3.3).This disaggregated Rg was multiplied with (Rsd)day to get λEday at 250 m.The schematic of the three approaches are presented in Figure 2 and the various models used are explained in the following sections.All the models used in this work were forced with inputs that are obtainable only from RS datasets without any ground measurements.In Figure 2, FVC and NDWI indicate, respectively, the fraction vegetation cover and normalized difference water index (sometimes referred as normalized difference moisture index), which was used in the disaggregation of TRad [17].

Estimation of EF and Rg from Satellite Data
The EF was estimated using the triangle model [18], which is based on the TRad-NDVI relationship.Several studies have successfully utilized the triangle model for the estimation of EF or λE (e.g., [19][20][21][22][23][24]) at various spatial scales.Over India, the triangle model for EF estimation was tested by [9,13] and a detailed review of the various studies that utilized the TRad-NDVI relationship was provided by [25].In this study, for the triangle model, the dry edge was determined using the automated algorithm proposed by [21] and the wet edge was determined as the minimum value of TRad at maximum NDVI [18] with NDVI being calculated from the first two bands of the MODIS surface reflectance data product.The TRad-NDVI context space was determined using all the MODIS pixels (at both 1000 m and 250 m resolutions) with in the spatial grids enclosing the tower sites.The EF from the triangle model was converted into Rg as given by: where the subscript 'inst' indicates the estimation of energy terms at the instant of satellite overpass.(Rn − G)inst was estimated as [11]:

Estimation of EF and R g from Satellite Data
The EF was estimated using the triangle model [18], which is based on the T Rad -NDVI relationship.Several studies have successfully utilized the triangle model for the estimation of EF or λE (e.g., [19][20][21][22][23][24]) at various spatial scales.Over India, the triangle model for EF estimation was tested by [9,13] and a detailed review of the various studies that utilized the T Rad -NDVI relationship was provided by [25].In this study, for the triangle model, the dry edge was determined using the automated algorithm proposed by [21] and the wet edge was determined as the minimum value of T Rad at maximum NDVI [18] with NDVI being calculated from the first two bands of the MODIS surface reflectance data product.The T Rad -NDVI context space was determined using all the MODIS pixels (at both 1000 m and 250 m resolutions) with in the spatial grids enclosing the tower sites.The EF from the triangle model was converted into R g as given by: where the subscript 'inst' indicates the estimation of energy terms at the instant of satellite overpass.
(R n − G) inst was estimated as [11]: where α is the surface albedo, ε s and ε a are surface and air emissivity, respectively, σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −1 ), and T a and T Rad are the air and land surface temperature (K), respectively.α was estimated from the MODIS seven-band surface reflectance data following [26] and the ε s was estimated from the emissivity of bands 31 and 32 in MODIS following [27].
In Equations ( 2) and ( 3), (R sd ) inst [11] and ε a [28] were estimated as: where a (=0.75) and b (=1.28) are empirical coefficients, S 0 is the solar constant (= 1367 W m −2 ), f is the Sun-Earth distance correction factor in astronomical units, and θ is the solar zenith angle (radians).The detailed explanation of the models used to estimate EF and (R n − G) is provided in [13].

Disaggregation of R g and EF
For the disaggregation of R g , a model called DiSoRa (DIsaggregation of SOlar radiation RAtio) is proposed based on the observed inverse triangular relationship between R g and NDVI.The DiSoRa is a contextual model which utilizes the variation of R g with different values of NDVI in a given satellite scene.The schematic of the model is given in Figure 3.
Remote Sens. 2017, 9, 949 5 of 13 where α is the surface albedo, εs and εa are surface and air emissivity, respectively, σ is the Stefan-Boltzmann constant (5.67 × 10 -8 W m -2 K -1 ), and Ta and TRad are the air and land surface temperature (K), respectively.α was estimated from the MODIS seven-band surface reflectance data following [26] and the εs was estimated from the emissivity of bands 31 and 32 in MODIS following [27].In Equations ( 2) and ( 3), (Rsd)inst [11] and εa [28] were estimated as: 9.2 10 where a (=0.75) and b (=1.28) are empirical coefficients, S0 is the solar constant (= 1367 W m -2 ), f is the Sun-Earth distance correction factor in astronomical units, and θ is the solar zenith angle (radians).The detailed explanation of the models used to estimate EF and (Rn − G) is provided in [13].

Disaggregation of Rg and EF
For the disaggregation of Rg, a model called DiSoRa (DIsaggregation of SOlar radiation RAtio) is proposed based on the observed inverse triangular relationship between Rg and NDVI.The DiSoRa is a contextual model which utilizes the variation of Rg with different values of NDVI in a given satellite scene.The schematic of the model is given in Figure 3.In Figure 3, the minimum value of Rg for a given value of NDVI increases with increasing NDVI (indicated by hypotenuse AB).The scatter above the minimum Rg line is due to soil moisture variation and the effect of soil moisture is higher (indicated by larger scatter) under lower canopy cover and it decreases with the increase in NDVI (indicating an increase in canopy cover).This kind of inverse relationship between EF and NDVI was observed by [11] and was utilized in the development of the DEFrac model for the disaggregation of EF [9].The DiSoRa model is a direct extension of the EF-NDVI relationship to the Rg-NDVI relationship.
In the application of the DiSoRa model, it is assumed that, for a given satellite scene, the minimum and maximum bounds (AB and BC, respectively, in Figure 3) can be established between Rg and NDVI, and these bounds are scale independent.For the MODIS scale data at 1000 m resolution, a linear regression was performed between minimum vales of Rg and NDVI to determine the slope (m) and intercept (c) as given below: In Figure 3, the minimum value of R g for a given value of NDVI increases with increasing NDVI (indicated by hypotenuse AB).The scatter above the minimum R g line is due to soil moisture variation and the effect of soil moisture is higher (indicated by larger scatter) under lower canopy cover and it decreases with the increase in NDVI (indicating an increase in canopy cover).This kind of inverse relationship between EF and NDVI was observed by [11] and was utilized in the development of the DEFrac model for the disaggregation of EF [9].The DiSoRa model is a direct extension of the EF-NDVI relationship to the R g -NDVI relationship.
In the application of the DiSoRa model, it is assumed that, for a given satellite scene, the minimum and maximum bounds (AB and BC, respectively, in Figure 3) can be established between R g and NDVI, and these bounds are scale independent.For the MODIS scale data at 1000 m resolution, a linear regression was performed between minimum vales of R g and NDVI to determine the slope (m) and intercept (c) as given below: R min g,1000 = (m The regression was carried out in an automated way following the edge detection algorithm in [21] with suitable modifications.Then, for any given pixel i in the 1000 m image, the fraction (d) by which the R g of that pixel (R g,1000 ) deviates from the minimum bound was estimated as: where the minimum R g was estimated using NDVI of that pixel in Equation ( 6) and R max g,1000 is the maximum value of R g found in the scene.Then, for any subpixel j at 250 m resolution with in the given 1000 m pixel i (j will vary between 1 and 16), the disaggregated R g at 250 m ( Rg, 250 ) resolution was estimated as: where, R min g, 250 [j] was estimated using the NDVI at 250 m resolution in Equation ( 6).From Equation (8), it can be observed that the fraction d is kept constant for all the 16 subpixels within a 1000 m pixel.Further, it should also be noted that, though the fraction d is kept constant for all the subpixels, the actual soil moisture term (part of equation within the curly braces {}) will vary based on the NDVI of that subpixel j.For a pixel with higher value of NDVI, the added soil moisture control term will be less and vice versa.Since the disaggregation model is contextual in nature, for the successful implementation of the model, the satellite scene should have a wide range in NDVI and pixels evaporating at potential and minimal rates should be found within the scene.The NDVI required at 250 m was estimated from the MODIS two-band surface reflectance product available at 250 m.For the disaggregation of EF to 250 m, the DEFrac model was used by replacing R g with EF in Equations ( 6)- (8).Both the R g -NDVI context space and EF-NDVI context space were established over the same spatial grids (~100 km × 100 km) which are used in the triangle model enclosing the tower site.

Disaggregation of T Rad
For the estimation of (R n − G) at 250 m (for Approaches 1 and 3) and for the estimation of EF at 250 m (for Approach 3), it was necessary to disaggregate the MODIS T Rad from 1000 m to 250 m.For this purpose, the DisTrad model [5,6] was utilized.The DisTrad thermal sharpening model is based on the relationship between T Rad and a vegetation index.Originally, NDVI and FVC were proposed and used in the model, however, a recent study [17] showed that the choice of vegetation index in the DisTrad model depends on the dryness/wetness of the study area.Only during wetter conditions, NDVI/FVC yielded more accurate results in the DisTrad model.For drier conditions, NDWI was found to be the better-performing vegetation index.In this study, too, the choice of index was decided by the temperature vegetation dryness index (TVDI), as shown in [17].For the sake of clarity, the DisTrad thermal sharpening model is not explained here and the readers are directed to refer to [17] for the complete information about the DisTrad model and its application.

Estimation of λE day at 250 m Resolution
After the estimation of EF and R g at 250 m using the three approaches, λE day was estimated from them at 250 m resolution.The EF was converted into λE day using (Approaches 1 and 3): Here the subscripts '250'and 'day' indicate the target resolution of 250 m and daytime integrated values of the energy fluxes respectively.The (R n − G) inst estimated from Equations ( 2) and (3) was temporally extrapolated to daytime using a simple sinusoidal extrapolation equation [29]: where t pass , t rise , and t set are, respectively, the time of the satellite overpass, sunrise, and sunset.The sunrise and sunset times were estimated based on standard equations using the latitude and longitude of the satellite pixel.For the estimation of (R n − G) inst at 250 m resolution, all the necessary inputs were converted into 250 m as given in Table 1.In Approach 2, λE day at 250 m was estimated from the disaggregated R g .The (R sd ) day was modelled using Equation (10) after substituting (R n − G) inst with (R sd ) inst (Equation ( 4)).
Table 1.Getting the input variables for estimating (R n − G) day at 250 m resolution.

Variables Data Source with Available Resolution Method Adopted to Get the Variable at 250 m Resolution
Solar zenith angle (θ) Available from MODIS 7-band surface reflectance product at 1000 m resolution.
Assumed the same value for all the 250 m pixels within a single 1000 m pixel.
Available from MODIS at 250 m Surface emissivity (ε s ) Available at 1000 m resolution from MODIS LST product.But not used for disaggregation Estimated using NDVI at 250 m following [30] Air temperature (T a ) Available at 5000 m resolution from MODIS atmospheric profile product.
Assumed the same value for all the 250 m pixels within a single 5000 m pixel.
T Rad Available at 1000 m resolution from MODIS LST product.
Disaggregated to 250 m using the DisTrad thermal sharpening model [17] Satellite overpass time (t pass ) Available at 1000 m resolution from MODIS LST product.
Assumed the same value for all the 250 m pixels within a single 1000 m pixel.

In Situ Data Processing
The estimation of λE day from ground observation was carried out using the BREB method.Consistency checks were performed to the micrometeorological observations and any anomalous and missing data points were excluded.Air temperature and relative humidity observed at two different levels above the ground were used to estimate the Bowen ratio (β).Observations made at 2 m and 4 m heights were used when the crop height in the catch of the tower was less than 2 m whereas, when the crop height was greater than 2 m, observations made at 4 m and 6 m were used.
β was estimated every half hour and subjected to the checks prescribed by [31].This half-hourly β was further averaged during the entire daytime to get daytime averaged Bowen Ratio (β day ) and from this, daytime averaged EF was calculated as 1/(1 + β day ).In this study, daytime was defined as the time during which the net shortwave radiation was positive.Using the measurements from the four component net radiometer and soil heat flux plates, (R n − G) was estimated at half-hourly intervals and averaged over the entire daytime to get (R n − G) day and similarly, the observations of incoming shortwave radiation were averaged over daytime to obtain (R sd ) day .λE day was estimated by multiplying the daytime averaged EF with (R n − G) day and finally the λE day was divided by (R sd ) day to get daytime averaged solar radiation ratio R g .These daytime averaged EF, R g and λE day were used to validate the satellite-based estimates.

Validation of Disaggregated EF and R g at 250 m
The R g -NDVI relationship that was established for the satellite data during the disaggregation of R g was visually examined first to check for the proper determination of the minimum and maximum bounds.It was observed that R g and NDVI exhibited an inverse triangle relationship with the minimum bound clearly increasing with NDVI.As examples, the R g -NDVI scatter observed on selected days during the rabi (winter) cropping season over Sites 3 (6 January 2011), 4 (26 January 2011), and 5 (12 December 2012) are presented in Figure 4.This kind of clear scatter was observed for over 90% of the analyzed MODIS datasets.This suggested that the R g -NDVI relation, based on which the DiSoRa model is formulated holds well.

Validation of Disaggregated EF and Rg at 250 m
The Rg-NDVI relationship that was established for the satellite data during the disaggregation of Rg was visually examined first to check for the proper determination of the minimum and maximum bounds.After the visual examination, the disaggregated Rg and EF from the three approaches were validated with corresponding ground measurements.Though, the disaggregation of EF using Approaches 1 and 3 was already validated by [9] over four of the five sites in used in this study, it is revalidated here by repeating the entire exercise using additional satellite datasets and over one new site.This revalidation was to ensure the consistent performance of the DEFrac disaggregation model.Further, the DisTrad thermal sharpening model was improved from the one used in [9] by choosing NDWI for drier conditions and FVC for wetter conditions.The results of validation of disaggregated EF and Rg were analyzed in-terms of the root mean square error (RMSE), which was defined as follows: where, 'Modelled' denotes the satellite estimated land surface variable, 'Observed' indicates the same variable from ground observations, and n is the total number of datasets used.For all the five sites taken together, the RMSE in the satellite estimated Rg at 1000 m and 250 m (disaggregated using the DiSoRa model) resolutions was estimated to be 0.07 (17% of the mean observed Rg) and 0.06 (15% of the mean observed Rg), respectively.The coefficient of determination (R 2 ) also marginally improved from 0.82 to 0.83 for Rg after disaggregation.In the case of disaggregated EF, the RMSE was 0.10 (16% of the mean observed EF) for both EF at 1000 m (R 2 = 0.75) and 250 m (R 2 = 0.76) resolutions when the disaggregation was carried out using the DEFrac model (Approach 1).However, when Approach 3 was used for disaggregating the EF, the RMSE increased to 0.17 (28% of the mean observed EF) and R 2 decreased to 0.32.The comparison between EF and Rg estimated from ground measurements and the satellite data is presented in Figure 5 and the site-wise error values are listed in Table 2.After the visual examination, the disaggregated R g and EF from the three approaches were validated with corresponding ground measurements.Though, the disaggregation of EF using Approaches 1 and 3 was already validated by [9] over four of the five sites in used in this study, it is revalidated here by repeating the entire exercise using additional satellite datasets and over one new site.This revalidation was to ensure the consistent performance of the DEFrac disaggregation model.Further, the DisTrad thermal sharpening model was improved from the one used in [9] by choosing NDWI for drier conditions and FVC for wetter conditions.The results of validation of disaggregated EF and R g were analyzed in-terms of the root mean square error (RMSE), which was defined as follows: (11) where, 'Modelled' denotes the satellite estimated land surface variable, 'Observed' indicates the same variable from ground observations, and n is the total number of datasets used.For all the five sites taken together, the RMSE in the satellite estimated R g at 1000 m and 250 m (disaggregated using the DiSoRa model) resolutions was estimated to be 0.07 (17% of the mean observed R g ) and 0.06 (15% of the mean observed R g ), respectively.The coefficient of determination (R 2 ) also marginally improved from 0.82 to 0.83 for R g after disaggregation.In the case of disaggregated EF, the RMSE was 0.10 (16% of the mean observed EF) for both EF at 1000 m (R 2 = 0.75) and 250 m (R 2 = 0.76) resolutions when the disaggregation was carried out using the DEFrac model (Approach 1).However, when Approach 3 was used for disaggregating the EF, the RMSE increased to 0.17 (28% of the mean observed EF) and R 2 decreased to 0.32.The comparison between EF and R g estimated from ground measurements and the satellite data is presented in Figure 5 and the site-wise error values are listed in Table 2.  From Table 2, it can be observed that the disaggregation of EF using the DisTrad thermal sharpening model resulted in higher errors than the EF at 1000 m resolution.On the other hand, the disaggregation of EF and Rg using the DEFrac and the DiSoRa models, respectively, resulted in marginal improvements in RMSE (which is not seen in Table 2 due to rounding off) at all the five sites.The RMSE in both the disaggregated EF and Rg was within 15% of the mean ground observed values at all but one site (Site 2), where, the error was within 20% of the mean observed value.A relatively higher RMSE in EF was observed at Site 2 due to the mismatch between the landcover conditions of the catch of the flux tower and the MODIS satellite pixel [13].The disaggregation using the DisTrad thermal sharpening model (Approach 3) has introduced more errors in the EF estimated using the triangle model.A similar result was obtained in an earlier study [9] where it was mentioned that the relative positions of TRad pixels gets modified in the TRad-NDVI context space when disaggregated thermal data is used in the triangle model.Though the DisTrad model was improved by choosing a suitable vegetation index, the disaggregated TRad was not accurate enough to be used in the contextual models such as the triangle model for the estimation of disaggregated EF.

Validation of Disaggregated λEday at 250 m
The disaggregated EF and Rg were converted into λEday, which was compared against ground observations.For all five sites considered together, the RMSE in the 1000 m λEday was found to be 34 W m -2 .Similarly, the RMSE of the disaggregated λEday was 32 W m -2 , 30 W m -2 , and 48 W m -2 for Approaches 1, 2, and 3, respectively.The R 2 was determined to be 0.69, 0.71, and 0.37 for the three approaches taken in order, respectively.The RMSE values for individual sites are presented in Table 3.  From Table 2, it can be observed that the disaggregation of EF using the DisTrad thermal sharpening model resulted in higher errors than the EF at 1000 m resolution.On the other hand, the disaggregation of EF and R g using the DEFrac and the DiSoRa models, respectively, resulted in marginal improvements in RMSE (which is not seen in Table 2 due to rounding off) at all the five sites.The RMSE in both the disaggregated EF and R g was within 15% of the mean ground observed values at all but one site (Site 2), where, the error was within 20% of the mean observed value.A relatively higher RMSE in EF was observed at Site 2 due to the mismatch between the landcover conditions of the catch of the flux tower and the MODIS satellite pixel [13].The disaggregation using the DisTrad thermal sharpening model (Approach 3) has introduced more errors in the EF estimated using the triangle model.A similar result was obtained in an earlier study [9] where it was mentioned that the relative positions of T Rad pixels gets modified in the T Rad -NDVI context space when disaggregated thermal data is used in the triangle model.Though the DisTrad model was improved by choosing a suitable vegetation index, the disaggregated T Rad was not accurate enough to be used in the contextual models such as the triangle model for the estimation of disaggregated EF.

Validation of Disaggregated λE day at 250 m
The disaggregated EF and R g were converted into λE day , which was compared against ground observations.For all five sites considered together, the RMSE in the 1000 m λE day was found to be 34 W m −2 .Similarly, the RMSE of the disaggregated λE day was 32 W m −2 , 30 W m −2 , and 48 W m −2 for Approaches 1, 2, and 3, respectively.The R 2 was determined to be 0.69, 0.71, and 0.37 for the three approaches taken in order, respectively.The RMSE values for individual sites are presented in Table 3. From Table 3, it can be observed that the disaggregation of λE day to 250 m using Approaches 1 and 2 has marginally improved the RMSE at all five sites in comparison to λE day at 1000 m, and the performance of the first two approaches were similar at all five sites.However, disaggregation using Approach 3 had resulted in higher errors than 1000 m at all sites.Further analyses revealed that the RMSE in the satellite-estimated (R n − G) day was estimated to be 29 W m −2 and 27 W m −2 at 1000 m and 250 m resolutions, respectively, for all five sites considered together.This shows that the higher error in disaggregated λE day using Approach 3 was predominantly due to the error in the disaggregation of EF rather than disaggregation of (R n − G) day .
The spatial patterns of the λE day observed over Site 5 on a selected day (14 November 2012) in the rabi (winter) cropping season at 1000 m and 250 m is presented in Figure 6.From the Figure, it can be observed that spatial disaggregation to 250 m using Approaches 1 (Figure 6b) and 2 (Figure 6c) had improved the detail of the λE day image at 1000 m (Figure 6a) and at the same time preserved the spatial pattern observed in the 1000 m image.On the other hand, disaggregation using Approach 3 (Figure 6d) under-predicted λE day values especially in the west-southwest parts of the spatial grid.From Table 3, it can be observed that the disaggregation of λEday to 250 m using Approaches 1 and 2 has marginally improved the RMSE at all five sites in comparison to λEday at 1000 m, and the performance of the first two approaches were similar at all five sites.However, disaggregation using Approach 3 had resulted in higher errors than 1000 m at all sites.Further analyses revealed that the RMSE in the satellite-estimated (Rn − G)day was estimated to be 29 W m -2 and 27 W m -2 at 1000 m and 250 m resolutions, respectively, for all five sites considered together.This shows that the higher error in disaggregated λEday using Approach 3 was predominantly due to the error in the disaggregation of EF rather than disaggregation of (Rn − G)day.
The spatial patterns of the λEday observed over Site 5 on a selected day (14 November 2012) in the rabi (winter) cropping season at 1000 m and 250 m is presented in Figure 6.From the Figure, it can be observed that spatial disaggregation to 250 m using Approaches 1 (Figure 6b) and 2 (Figure 6c) had improved the detail of the λEday image at 1000 m (Figure 6a) and at the same time preserved the spatial pattern observed in the 1000 m image.On the other hand, disaggregation using Approach 3 (Figure 6d) under-predicted λEday values especially in the west-southwest parts of the spatial grid.These results suggest that both Approaches 1 and 2 can be used for spatial disaggregation of λE day .However, it can also be observed that the disaggregation using the DiSoRa model (i.e., Approach 2) is relatively simpler and has fewer steps than disaggregation using the other two approaches.Especially, the additional disaggregation of T Rad to estimate (R n − G) can be avoided if disaggregation based on R g (DiSoRa model) is carried out.Furthermore, the results from the R g -based disaggregation can be improved by either calibrating Equation (4) to the specific region of interest, or using other reliable data/models to obtain R sd instead of using the empirical equation.
The major limitation of the proposed disaggregation Approaches 1 and 2 is the accuracy of the disaggregated λE day is limited by the accuracy of EF/R g at 1000 m resolution.Hence, it becomes necessary to test a model for estimation of λE day at coarse spatial resolution over the region of interest and then use this disaggregation approach.For example, the EF and radiation models used in this study were chosen after comparison with other models over all the five sites [13].Further, in this study, data from only one Bowen ratio tower was available in each of five spatial grids.Ideally, data from multiple spatially-distributed towers within each grid would be required for performing a robust validation of the disaggregation approaches.This was not possible now as there are only a handful of flux towers available across India.However, with increasing number of hydrological observatories in India, in the near future, it should be possible to validate these approaches with data from multiple flux towers within a given spatial grid or river watershed.

Conclusions
This paper compared three simple approaches for the disaggregation of latent heat flux (λE).The first two approaches are, respectively, based on the disaggregation of the evaporative fraction (EF, DEFrac model) and the solar radiation ratio (R g , DiSoRa model), and the third approach is based on disaggregation of thermal data (T Rad ) for estimation of finer resolution latent heat flux.The three approaches were tested using MODIS datasets over five sites in India and the validation was carried out using data from Bowen Ratio energy balance (BREB) micrometeorological towers.The validation suggested that utilizing the disaggregated T Rad for the estimation of disaggregated λE (RMSE = 48 W m −2 ) resulted in larger errors than the original estimates at 1000 m (RMSE = 34 W m −2 ).The results further revealed that the higher error in disaggregated λE from this approach is due to errors in the disaggregated EF rather than the errors in disaggregated (R n − G).In this study, the triangle model was used for the estimation of EF, and the accuracy of the estimated EF from the triangle model is sensitive to the relative position of T Rad of a pixel with respect to the wet and the dry edges in the T Rad -NDVI context space.Using the disaggregated T Rad in the triangle model had perturbed this relative position of pixels in the T Rad -NDVI context space resulting in higher errors.
The approaches based on the disaggregation of EF and R g produced equally better results at all five sites, and marginally improved the error values of the disaggregated λE in comparison with 1000 m estimates.The RMSE of the disaggregated λE was 32 W m −2 and 30 W m −2 , respectively, for EF and R g -based disaggregation.Both these approaches have performed in a similar manner as they share the same physical basis.However, by converting EF into R g , the need to estimate (R n − G) at finer spatial resolution is avoided, thereby making Approach 2 simpler with fewer data processing steps.Future studies should focus on utilizing the DiSoRa and the DEFrac models for the disaggregation of λE estimated from models like the SEBAL, two source model, etc. and re-validating these approaches at multiple sites across the globe.

Figure 1 .
Figure 1.Location of the five BREB tower sites.

Figure 1 .
Figure 1.Location of the five BREB tower sites.

Figure 2 .
Figure 2. Flowchart indicating the steps involved in estimating λEday at 250 m using the three approaches.

Figure 2 .
Figure 2. Flowchart indicating the steps involved in estimating λE day at 250 m using the three approaches.
It was observed that Rg and NDVI exhibited an inverse triangle relationship with the minimum bound clearly increasing with NDVI.As examples, the Rg-NDVI scatter observed on selected days during the rabi (winter) cropping season over Sites 3 (6 January 2011), 4 (26 January 2011), and 5 (12 December 2012) are presented in Figure 4.This kind of clear scatter was observed for over 90% of the analyzed MODIS datasets.This suggested that the Rg-NDVI relation, based on which the DiSoRa model is formulated holds well.

Figure 5 .
Figure 5.Comparison between ground observed and disaggregated (a) EF using the DEFrac model, (b) Rg using the DiSoRa model, and (c) EF using the DisTrad thermal sharpening model over the five sites.

Figure 5 .
Figure 5.Comparison between ground observed and disaggregated (a) EF using the DEFrac model, (b) R g using the DiSoRa model, and (c) EF using the DisTrad thermal sharpening model over the five sites.

Figure 6 .
Figure 6.Spatial pattern of λEday observed on 14 November 2012 over Site 5. (a) λEday at 1000 m resolution and (b-d) λEday disaggregated to 250 m using approaches 1, 2 and 3 respectively.Figure 6. Spatial pattern of λE day observed on 14 November 2012 over Site 5. (a) λE day at 1000 m resolution and (b-d) λE day disaggregated to 250 m using approaches 1, 2 and 3 respectively.

Figure 6 .
Figure 6.Spatial pattern of λEday observed on 14 November 2012 over Site 5. (a) λEday at 1000 m resolution and (b-d) λEday disaggregated to 250 m using approaches 1, 2 and 3 respectively.Figure 6. Spatial pattern of λE day observed on 14 November 2012 over Site 5. (a) λE day at 1000 m resolution and (b-d) λE day disaggregated to 250 m using approaches 1, 2 and 3 respectively.

Table 2 .
Results of comparison between ground observations, disaggregated EF, and disaggregatedRg.

Table 2 .
Results of comparison between ground observations, disaggregated EF, and disaggregated R g .

Table 3 .
Results of the comparison between ground observations and satellite-obtained λE day .

Table 3 .
Results of the comparison between ground observations and satellite-obtained λEday.