Applied Different Pixel Selection in METRIC Model for Estimating Spatial Daily Evapotranspiration of Oil Palm in East Kalimantan Province, Indonesia

Abstract


INTRODUCTION
Elaeis guineensis, so-called oil palm, is an agricultural commodity with high economic value.In the past decades, the supply-demand of oil palm products has increased drastically and the palm oil industry, which has a continual prospect, led to the expansion of oil palm plantations chiefly in Indonesia [38].Indonesia has become one of the largest countries with major land use for oil palm plantations, especially in Sumatra and Kalimantan islands.In 2019, East Kalimantan Province had an oil palm plantation that covers around 1,254.20 ha of the total oil palm area in Indonesia (14,456.60 ha).The oil palm coverage is expanding every year, which in 2020 increased by about 4.74% from 2019 and in 2021 increased by around 3.99% from 2020 [2].Therefore, along with increasing the demand for oil palm products to maintain and enhance the oil palm productivity every year, water crop requirements for oil palm become a vital study to be monitored and evaluated, which can be achieved by quantifying and estimating the evapotranspiration as a loss of water amount [11,20].
Estimating evapotranspiration as a residual surface energy balance can be carried out with remote-based sensing to monitor the agricultural area [7].There are several energy balances models that have been developed for use in estimating evapotranspiration, such as Surface Energy Balance Algorithm for Land (SEBAL) [32,35], Two Sources Energy Balance (TSEB) [9,36], and Mapping Evapotranspiration at High Resolution with Internalized Calibration (METRIC) [21].In particular, METRIC has advantages for mapping evapotranspiration which considers the separated energy balance calculation and pixel selections, which are hot and cold pixels, for more accurate results [8,23].
SEBAL model has been conducted in oil palm plantation by [24], which depicted latent heat fluxes as the energy used to evaporate water from the surface and plant, had a high value in mature oil palm than in young oil palm.[27] denoted the appropriate spatial daily evapotranspiration using the SEBAL model, which was validated by the aerodynamics method from eddy covariance flux tower observation.However, TSEB and METRIC models have not been conducted to estimate the evapotranspiration in oil palm plantations.Thus, we concern with (i) estimating and depicting the spatial daily evapotranspiration of oil palms at different ages and (ii) analyzing the statistical value between different treatment of pixel selections (with and without pixel selection), in particular of the utilization of the METRIC model in PT Teladan Prima Agro, East Kalimantan Province.

METHOD 2.1 Area of Study
The study was conducted at the Muara Bengkal Estate, PT.Teladan Prima Agro (See figure 1), which located at latitudes of 0˚ 28' -0˚ 36' N and longitudes of 116˚ 44' -116˚ 50' E in East Kalimantan Province Province, Indonesia.

Plotting Design of Oil Palm Area
We performed our research on 4 sample parts of the plantation area, so-called afdeling: afdeling 01, 02, 03, and 04 and seven blocks, which were blocks G, J, M, N, O, P, and T (see Figure 2).Afdeling 01 was consisted of Block G, M and N. Block G had 138.13 ha of the planted area from 143.71 ha surface area, which had a flat land topography, consisting of oil palm (planted in July and August 2007), road, pond, and marsh.Block M had 91.95 ha of the planted area from 91.95 ha of surface area, which had a wavy land topography, consisting of oil palm (planted in September and November 2007), road and marsh.Block N had 14.39 ha of the planted area from 14.39 ha of surface area, which had a wavy land topography, consisting of oil palm (planted in January 2007), road, forest and marsh.
Afdeling 02 was consisted of block J. Block J had 39.19 ha of the planted area from 40.95 ha of surface area, which hds a flat land topography, consisting of oil palm (planted in October 2009) and road.Moreover, Afdeling 03 was consisted of Block P and T. Block P had 104.34 ha of the planted area from 153.79 ha of surface area, which had a hills land topography, consisting of oil palm (planted in January 2008), road, marsh, steep road, emplacement, and pond.Block T had 97.36 ha of the planted area from 117.34 ha of surface area, which had a hills land topography, consisting of oil palm (planted in March and June 2008), road, steep road, and forest.Subsequently, Afdeling 04 was consisted of Block O. Block O had 36.43 ha of the planted area from 69.18 ha of surface area, which had a hills land topography, consisting of oil palm (planted in August and September 2011), road, marsh, steep road and forest.

Data Analysis
Our research was performed using several applications, viz.ArcGIS 10.8.1, Python 3.6.8,ENVI Software 5.3, and Microsoft Office 2019 to visualize, calculate, and analyze the data.We used Landsat 8 OLI/TIRS with path/raw 125/61, acquisition date: 06/04/2015 and 03/04/2020.We also used ERA-5 Reanalysis data, which included an hour of air temperature, dew-point temperature, air pressure, and wind speed.The ERA-5 Reanalysis also provided incoming and outgoing short and longwave radiation.We used elevation data obtained from Google Earth Pro.
METRIC model is the model-based energy balance and remote sensing where calculated evapotranspiration from latent heat fluxes as the residual energy balance components.The calculation will show in equation 1 below: ) where, RN = net radiation (Wm -2 ), G = ground heat flux (Wm -2 ), H = sensible heat flux (Wm -2 ), and LE = latent heat flux (Wm -2 ).

Radiometric Calibration
Landsat 8 Imagery OLI/TIRS mostly in the form of digital number format.Radiometric calibration is necessary to alter from the digital number to spectrum radiance (equation 2), so the imagery can be calculated further.The alteration uses radiometric calibration that can utilize gain and offset parameters, which is available on image metadata [30].Subsequently, the cloud covers are removed in the Landsat 8 imagery.

Calculation of basic parameters in pixel selection
Hot and cold pixels selected by 2 parameters: Normalized Different Vegetation Index (NDVI) and Land Surface Temperature (LST).NDVI illustrates the Proportion of Vegetation (PV), which measure the difference between Near-Infrared (Red) electromagnetic spectrum reflected (absorbed) by vegetation.Reference [14], calculation of NDVI in Landsat 8 uses band 5 (NIR) and band 4 (Red), which is the equation will show in the following: where  = Reflectance of Near-Infrared band (band 5 in Landsat 8), = Reflectance of Red band (band 4 in Landsat 8) LST can be calculated by converse the brightness temperature [4,14], where the conversion conducted in order to change radiance spectral to temperature value in Kelvin or Celcius [30].

Calibration of Sensible Heat Flux
Sensible heat flux (H) is an energy balance component with a complex algorithm.H depicts energy changes imposed by alteration in temperature at different high above the plant [19,21,32,33,34,35].There are several supplementary calculations to calibrate H; they are roughness length, wind speed at 200 meters, friction velocity, aerodynamics resistance and temperature gradient.
=  +     (12) where  0 = roughness length (m), ℎ = plant height (m),  1 and  1 = constant regression from ln( 0 ) to    [31] where   = surface albedo [22],  200 = wind speed at 200 meters (m/s),  = wind speed at 2 meters (m/s),  = Von Karman constant (0,41),  * = friction velocity (m/s),  ℎ = Aerodynamics resistance (s m -1 ) [29],  1 = height 1,  2 = height 2,  = temperature gradient,   = LST (high and low LST in Celcius),  and  constitute correlation coefficient between dT and   [21].METRIC model has advantages to calibrate the sensible heat fluxes internally using pixel selection (see equations 13 and 14), where the cold and hot pixel is defined as a potential condition of evapotranspiration and hardly any evapotranspiration at the surface, respectively.The cold and hot pixel candidature can be determined using Normalized Different Vegetation Index (NDVI) and Land Surface Temperature (LST).Both pixels can be found at 1% and 2% in the area of interest, respectively.Most of the candidates' cold pixels were found in the agricultural field, whereas bare soil represents the hot pixel [1,3,17].Generally, 1% of the isolated area (candidate of the cold pixel) can be obtained by selecting 5% of the coldest LST and 20% of the highest NDVI.Besides, the hot pixel is founded in 10% of the warmest LST and 20% of the lowest NDVI [17].However, due to many pixels in the area of interest, which will impact the model run time, so we adjust the criteria percentage of NDVI from 20% to 15% reference to [37].On the other hand, without pixel selection uses the average of LST from all the pixels in study area.
= ( − ) −   (14) where  ℎ = sensible heat flux at hot pixel (Wm -2 ),   = sensible heat flux at cold pixel (Wm -2 ),  ℎ = latent heat flux at hot pixel (Wm -2 ), and   = latent heat flux at cold pixel (Wm -2 ) [21]. ℎ is assumed to zero, while   is calculated by equations below: In actual conditions, the mass and heat transfers at a small-scale range take place vertically and horizontally.However, the earth's surface has horizontally inhomogeneities that impact radiation distribution and other energy budget components.The other components (sensible, latent, and ground heat fluxes) are responses to the radiative forcing.Subsequently, the partitioning of the ratio's G/RN, H/RN, and LE//RN is expected to depend on micrometeorological data.Therefore, the assumption of G that is proportional to RN through empirical regression restricts the transfer of heat and mass only occurring vertically [26].
In the preceded foundation of the METRIC model, so-called SEBAL, an evaporative fraction (EF) as the ratio of available energy (RN-G) was used to simplify the energy balance.EF is defined to be the same for the 24 hours period.However, increasing advection may occur during the day, resulting in increased evapotranspiration in the proportion of RN-G.Thus, METRIC uses Reference Evapotranspiration Fraction (ETrF), to be constant in a day, which is capturable any impact of increasing advection during the day [21].ETrF is calculated using daily evapotranspiration (accumulation) and hourly evapotranspiration when Landsat-8 passes (equation [17][18][19]. = 3,600       (19) where  24 = daily evapotranspiration from METRIC model (mm d -1 ), ETrF = reference evapotranspiration fraction (mm h -1 ),  −24 = accumulation of reference evapotranspiration (mm d -1 ),   = instantaneous evapotranspiration (mm h -1 ), ETr = reference evapotranspiration fraction (mm h -1 ),  = water density (~1,000 kg m -3 ) [21].
Based on [21], the METRIC model uses the standardized ASCE Penman-Monteith equation as reference evapotranspiration.However, due to different characteristics of the surface in oil palm plantations, the principle of physics method so-called Penman-Monteith ( 0 ) from Food and Agricultural Organization (FAO) 56 was used as reference evapotranspiration (ETr) [27], which appropriated to estimate the heat flux in oil palm plantation either at young or adult ages [15].

Aerodynamic Transport
The sensible heat flux rate drives the buoyancy forces within the boundary layer to influence the aerodynamics resistance (rah).In the METRIC model, an iterative solution is used to estimate the rah, where friction velocity (u*) was computed at the first iteration.The first iteration used logarithmic wind law for neutral atmospheric conditions [21], which is the tangential rate of air parcel movement impacted by mechanical turbulence so that the neutral condition is favorable to determine u* when buoyancy forces less significance [28].The u* was impacted by the surface roughness, thereby figure 3 illustrates the difference between aerodynamic resistance (rah) and friction velocity (u*) values at several wind speed measurement heights.The friction velocity witnessed a decline along with rising wind speed measurements height until it reached constant at 200 meters.Decreasing the u* impacted increasing of rah until reached constant values when the u* has stable rates.The u* elevates the rate of the eddy diffusion process of latent heat fluxes.In unstable (stable) conditions, eddy diffusion increases (decreases), which will impact more (less) water vapor transfer [29].Reference [21], wind speed measurement height at 200 meters is assumed to be constant over all the pixels in the satellite imagery, which is unaffected by the surface features.The constant values define the relation between temperature gradient (dT) and surface temperature (LST) that can extend the entire imagery.Therefore, the utilization of wind speed measurement height at 200 m to calculate u* in the METRIC model is appropriate to describe the neutral condition of rah and also necessary for determining of specific roughness length for each pixel.
The u* in 2015 (oil palm aged 4, 6, 7, and 8 years) is higher than in 2020 (oil palm aged 7, 9, 12, and 13 years).They are from 0.10 -0.13 ms -1 and 0.11 -0.14 ms -1 , respectively.Then, it impacts the lower rah of young oil palm than mature oil palm, which are 52.85 -66.54 sm -1 and 56.29 -70.39 sm -1 , respectively.According to [28], young oil palm had a very small roughness length, which is a smooth surface that did not have a large absorbed momentum area.Besides, mature oil palm absorbed more momentum due to rough surface characteristics and triggered bigger turbulence production.Therefore, it imposed greater shear stress in mature than young oil palm and indicate the increasing u*, which had a strong and positive correlation to wind speed.

Pixel Selection and Spatial Daily Evapotranspiration
The principle of NDVI is concerned to plant growth effectiveness with increased absorbed radiation until it reaches a particular point.NDVI uses the reflected NIR wavelength that can measure the photosynthetically plants on the surface and at the various structure of different vegetation covers [12].Reference [24], the vegetation structure of young oil palm has less vegetation cover and canopy compared to mature oil palm.Thereby, it will impact the high surface temperature on the plantation.
Table 1 shows the statistical parameter such as minimum, maximum, mean and standard deviation of NDVI.The NDVI value at the consecutive imagery is in the range 0.39 -0.95 (4, 6, 7 and 8 years-oil palm) and 0.48 -0.96 (9, 11, 12 and 13 years-oil palm).Average values show mature oil palms have high NDVI than young oil palms.Young oil palms at 4 years have the lowest NDVI, which is 0,838 ± 0,073 due to less dense vegetation and tend to be open land.LST also plays an important role in the METRIC model as a basic determination of cold and hot pixels, which is assumed to have maximum evapotranspiration (low LST) and hardly any evapotranspiration (high LST).However, the surface condition and cloud covers must be considered.Therefore, not all the pixels with high LST represent the hot pixel.On the contrary, LST on oil palm plantations can be the cold pixel candidate representing the high vegetation cover.
LST in Table 2 shows average values of oil palm in 2015 and 2020; 23.56 ℃ -25.66 ℃ and 20.59 ℃ -22.87 ℃, respectively.The lowest LST showed by oil palms aged 12 years at block P, which is 21.41± 0.15 ℃.On the other hand, the highest LST is shown by oil palms aged 4 and 6 years, which is 24.51 ± 0.34 ℃ dan 24.51 ± 0.32 ℃, Pixel selection in the METRIC model must determine the location of the cold and hot pixels.The cold and hot pixels can represent the surface condition and above.The utilization of METRIC model needs to determine the anchor pixel to restrict the candidatures both pixels in each satellite imagery.However, errors may occur from the users.Therefore, the applied criterion in the selection process using NDVI and LST can decrease user error [17].The candidatures of cold and hot pixels are shown in Figure 4. NDVI values of the hot pixels (cold pixels) candidate at the consecutive imagery are 0.35 -0.45 (0.93 -0.94) and 0.30 -0.40 (0.93 -0.95).LST values of the hot pixels (cold pixels) candidate are 301.00K-302.00K(296.75K-297.00K) and 298.00K -299.00K(294.00K-294.50K).The amount of daily evapotranspiration at several oil palm ages received from the METRIC model with and without pixel selection in the calibration process is given in Figure 5.The automated calibration in the METRIC model uses reference evapotranspiration hourly and daily.Both values are used for scaling the instantaneous evapotranspiration to daily evapotranspiration.The hourly reference evapotranspiration at consecutive image acquisition is 0.52 mm h -1 and 0.69 mm h -1 .The daily reference evapotranspiration or accumulation of evapotranspiration 24 hours period is 4.28 mm d -1 and 5.07 mm d -1 , respectively.

Figure 1 .
Figure 1.The geographic location of the research site in the Muara Bengkal Estate

Figure 2 .
Figure 2. The area of oil palms at the different planted years

Figure 3 .
Figure 3. Aerodynamics resistance and friction velocity of oil palms at different measurement heights

Figure 4 .
Figure 4. Pixel selection in the METRIC Model of (a) hot pixels candidate and (b) cold pixels candidate

Table 1 .
Statistical parameter of NDVI at different oil palm ages