Spatial variability of coffee plant water consumption based on the SEBAL algorithm

Awareness of evapotranspiration (ET) and crop coefficient (Kc) is necessary for irrigation management in coffee crops. ET and Kc spatial variabilities are disregarded in traditional methods. Methods based on radiometric measurements have potential to obtain these spatialized variables. The Kc curve and spatial variability of actual evapotranspiration (ETa) were determined using images from Landsat 8 satellite. We used images of young and adult coffee plantations from OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor) sensors over a two-year period. Evapotranspiration was estimated using the Surface Energy Balance Algorithm for Land (SEBAL). Moreover, the reference evapotranspiration (ETo) was estimated through the Penman-Monteith method. We obtained the values for the evapotranspiration fraction (ETf), analogous to Kc, according to ET and ETo values. The study was conducted in Buritis, Minas Gerais State, Brazil, in areas cropped with Coffea arabica irrigated by central pivots. A comparative analysis was made using different statistical indices. Average ETa was 2.17 mm d –1 for young coffee plantations, , and the Kc mean value was 0.6. For adult coffee plantations, average ETa was 3.95 mm d–1, , and the Kc mean value was 0.85. The ETc and Kc data obtained based on the SEBAL algorithm displayed similar values to studies that used traditional methods. This model has huge potential to estimate ET of different stages of coffee plantation for the region studied.


Introduction
Crop evapotranspiration (ET c ) in non-deficit cultivation areas is determined by the crop coefficient (K c ) level, which correlates ET c and the referential evapotranspiration (ET o ).However, evapotranspiration (ET) is a complex function of soil properties, atmospheric conditions, soil use, vegetation and topography, and is influenced by these parameters in space and time (Navarro et al., 2016).
In the field, ET c can be measured along a homogeneous surface, using conventional techniques such as the Bowen ratio, Eddy covariance, soil water balance and lysimetric procedure.However, in situ measurements are limited in generating area estimates in terms of cost and accuracy, because of natural heterogeneity (Allen et al., 2011;Kiptala et al., 2013).Thus, estimating or measuring the variable on a regional scale is difficult.
With the advent of satellites for Earth observation (Navarro et al., 2016), numerous ET-based remote sensing algorithms have been developed and validated (Paul et al., 2013).One of these methods is the SEBAL (Surface Energy Balance Algorithm for Land) (Bastiaanssen et al., 1998a, b), developed in the 1990s.The main advantage of remote sensing is that ET can be calculated on a large scale, on a pixel-to-pixels basis, in a consistent set of equations (Teixeira et al., 2009).
This method has demonstrated great potential in estimating ET of several crops in large areas and in different regions worldwide, using limited meteorological data (Bala et al., 2016;Bastiaanssen et al., 2005;Kiptala et al., 2013;Li et al., 2008;Mkhwanazi et al., 2015;Paul et al., 2013;Rawat et al., 2017;Ruhoff et al., 2012;Sari et al., 2013;Teixeira et al., 2009).Some studies on ET were conducted in coffee growing areas (Flumignan et al., 2011;Marin et al., 2005); however, estimation of water consumption using SEBAL algorithm in coffee cultivation should be better studied, since remote sensing can improve irrigation management considering the spatial variability of cultivated areas.
This study aimed to estimate spatial distribution of ET c and K c curve of coffee plantation, using Landsat 8 sensor and SEBAL algorithm.It also aimed validate spatial estimates of ET c using FAO guidelines and other reference studies for K c of coffee plantations in Brazil.

Location and characterization of the study site
The study was carried out in coffee crops from the northwestern regions of the state of Minas Gerais (MG), Brazil (Brazilian Cerrado).The farms are located in the municipality of Buritis and neighboring municipalities (15°48'S lat., 46°30'W long.and elevation of 600 m), totaling approximately two thousand irrigated hectares of Coffea arabica (Figure 1).The irrigation systems used in the area are central pivot-type, equipped with LEPA (Low Energy Precise Application) emitters.
The local weather is type Aw (Alvares et al., 2013), tropical with dry season in the winter.The average annual temperature is 24.3 °C and the average annual rainfall is 1,815 mm.The terrain is flat and the predominant soils were classified as Oxisol (Typic Ustox).

Satellite imaging
Images of young and adult coffee crops, from OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor), from the Landsat 8 satellite, were used during a two-year period (2014-2015) (Table 1).Regarding selection, images without clouds were prioritized and regularly distributed throughout the study period for a better evaluation of the variables (evapotranspiration and crop coefficient) according to the different phenological stages of coffee plants.An image of the ASTER GEDEM V2 satellite was also applied, referring to the digital terrain model (DTM).
Landsat 8 satellite images were re-routed into the Universal Transverse Mercator (UTM) coordinate system, Zone 23 South and Datum Sirgas 2000, using the EPSG 31983 code.In addition, the images were converted from digital numbers to reflectance values at the top of the atmosphere.These data preprocessing procedures were performed in the geographic information systems (GIS) Gdal (Geospatial Data Abstraction Library, version 2.0) and QGIS (QGIS Development Team, version 2.8.3), the latter using the GRASS GIS program (GRASS Development Team, version 7.0) plugin, which enables the interaction between both GISs.
To estimate ET at different temporal and spatial scales, SEBAL algorithm uses the residual approaches of surface energy balance.The net energy from the sun and atmosphere in the form of short and long-wave radiation is transformed and used for: (i) heating the soil (soil heat flux into the ground); (ii) heating the surface environment (sensible heat flux to the atmosphere), and; (iii) transforming water into vapor (latent heat flux from the crop/soil surfaces).All the energy involved in the soil-plant-atmosphere system can be given as the energy balance equation: where: R n is the net radiation, (W m −2 ); G is the soil hear flux, (W m −2 ); H is the sensible hear flux, (W m −2 ), and; λ ET is the latent head flux, (W m −2 ).Therefore, adjusting the energy balance equation and considering the λ ET as the "residual" energy, the ET is estimated as: The latent heat is expressed as hourly ET (mm) by dividing LE by the latent heat of vaporization and water density, multiplying by 3600 s h −1 .
Net radiation (R n ) represents the balance between input and output flux radiation, expressed as: where: R S ↓ is incoming shortwave radiation, (W m −2 ); α is surface albedo, (decimal); R L ↓ is incoming longwave radiation, (W m −2 ); ε0 is broadband surface thermal  The iteration of internal calibration is carried out to assign two evapotranspiration extreme conditions.One when H is equal to zero called "cold" pixel, and another for the latent heat flux (λET) equal to zero called "hot" pixel.
The "cold" pixel is selected on well-irrigated agricultural surfaces covered by vegetation, representing the maximum quantity of energy available that was consumed by the evapotranspiration process.The "hot" pixel is selected in fields with dry bare soil, where evapotranspiration is assumed zero (λET = 0) or H = Rn − G.The selection of the pixels is described in Allen et al. (2007).In this study, "hot" and "cold" pixels were selected manually.
The last stage of the SEBAL algorithm is the evapotranspiration estimate.Once calculating the energy balance fluxes, the instantaneous evapotranspiration (ET i ) at satellite overpass time is calculated as follow: where: ET i is the instantaneous evapotranspiration, (mm h −1 ); 3600 is the number of seconds at one hour, and; λ is the latent heat of vaporization, (J kg −1 ) stated by Allen et al. (2007) as: Therefore, daily (ET o ) and hourly (ET oi ) reference evapotranspirations were calculated through the Penman-Monteith method (Allen et al.,1998).After estimating the values of ET i , from SEBAL, and ET oi , the ratio between these variables provides a coefficient called reference evapotranspiration fraction (ET of ).When obtained on agricultural surfaces, this coefficient is similar to crop coefficient (K c ) (Tasumi et al., 2005).
ET of at hour of satellite overpass is relatively constant during the day, as highlighted by Allen et al. (2007).The daily or actual evapotranspiration (ET a ) is calculated as: where: ET a is the actual evapotranspiration, (mm d −1 ); ET of reference evapotranspiration fraction, (decimal), and; ET o is the reference of daily evapotranspiration (mm d −1 ).The calculation steps of the SEBAL algorithm were performed using a Python script SEBAL GRASS (Wolff, 2016), at GRASS GIS (GRASS Development Team, version 7.0).

Reference evapotranspiration (ET o ) and crop coefficient (K c )
The reference evapotranspiration for the satellite days of passage was estimated using the Penman-Monteith method, according to the methodology emissivity, (decimal), and; R L ↑ is outgoing longwave radiation, (W m −2 ).The term αR S ↓ corresponds to fraction of incoming longwave radiation reflected from the surface.
Net shortwave radiation available on Earth surface depends on R S ↓ and α, which was calculated from spectral radiance for each infrared and visible band, following the mathematical expressions for spectral integration and atmospheric correction given by Silva et al. (2016).R S ↓ was estimated from the theoretical radiation on top of atmosphere (R TOA ) and from atmosphere transmissivity (τ SW ), both expressions are found in Allen et al. (1998).R L ↑ was computed according to the Stefan-Boltzman equation (Boltzmann, 1884), which expresses a theoretical quantity of radiation emitted by the surface with the surface temperature on the fourth power.After the atmospheric correction and using the thermal band 10, the surface temperature (T s ) was assigned by Weng et al. (2004).In turn, ε0 was estimated through the Soil-Adjusted Vegetation Index (SAVI) and Leaf Area Index (LAI) according to Allen et al. (2007).
As noted for R L ↑, R L ↓ is estimated by the Stefan-Boltzman equation Boltzmann (1884), which is related to the empirical emissivity of the atmosphere (ε A ) and the fourth power of air temperature, as demonstrated by boundary conditions given by Allen et al. (2007).
R n was obtained after the incoming and outgoing radiation fluxes of soil-plant-atmosphere system were calculated.Soil heat flux (G) was calculated through the ratio G/R n adopted by Bastiaanssen (2000) and Allen et al. (2007), expressed as: where: T S is surface temperature, (K); NDVI is Normalized Difference Vegetation Index, (dimensionless).For water surfaces, where the NDVI value is negative, the ratio G/ R n is defined as 0.5.Sensible heat flux (H) is estimated from the near surface temperature gradient given by the aerodynamic equation: where: ρ is the air density, (kg m −3 ); c p is the specific heat of air at constant pressure, (≈ 1004 J kg −1 K −1 ); r ah is the aerodynamic resistance to head transfer, (s m −1 ), and; a and b are empirical coefficients determined through an internal calibration for each satellite image.
The term (a + bT S ) is an equation represented by the near surface temperature gradient and air temperature, for heights 0.1 m and 2 m above surface.Determining rah requires iterative calculus for stability correction of atmosphere using the Monin-Obukhov length (L) (Bastiaanssen, 2000;Koloskov et al., 2007).described by Allen et al. (1998) The meteorological data used in ET o estimates were obtained from the automatic meteorological station in the study site.This station was installed at 2 m above ground and recorded the time data, storing them in a data acquisition system.After obtaining ET o , the K c value was obtained from the ratio between actual and reference evapotranspiration.
The term crop coefficient (K c ) was applied throughout the study representing K c multiplied by the soil water stress coefficient (K s ) (Allen et al., 1998).If K c includes K s and K c , it influences the comparison with theoretical K c values of FAO 56.In other words, the soil was considered to be at field capacity on the assessed days.
The K c values obtained from the images were only considered for the following study stages in adult and young coffee crops.The mean K c for young coffee plants was calculated using values from areas of four central pivots that cover an area of 2,250,000 m 2 , corresponding to 2,491 pixels (30 m × 30 m).The mean K c for adult coffee plants, on the other hand, was calculated using values from the areas of 13 central pivots, which cover an area of 12,840,000 m 2 , corresponding to 14,267 pixels (30 m × 30 m).
Coffee crops were considered young with less than three years and adult with three years or more.The polygon delimitation of young and adult coffee crops considered edge pixels that entered neighboring fields.However, the ET c values of these edge pixels were only considered to calculate the averages when more than 50 % of their coverage area was within a certain pivot.

Data analyses
The monthly mean and standard deviation calculations were performed for each situation considering ET a and K c values obtained through the SEBAL method.ET monthly average was calculated based on data from 2014 and 2015.For young coffee crops, we used mean values of four pivots and for adult coffee plantations, 13 pivots were used to calculate mean values.
ET a values of young and adult coffee crops based on the SEBAL algorithm were compared to ET c (mm d -1 ) calculated from observed data.The measured ET c was obtained by multiplying the estimated ET o in the Penman-Monteith (PM) method by the K c of recommendations from the FAO-56 bulletin and reference studies in Brazil (Table 2).
Statistical analyses were performed to determine differences between ET c measured and estimated: Root Mean Square Error (RMSE), Mean Bias Error (MBE), Normalized Root Mean Square Error (NRMSE) and Index of Agreement (d) to know the significance between measured and estimated values.

Results and Discussion
Figure 2 shows some steps applied in the SEBAL algorithm to estimate actual evapotranspiration and coffee crop coefficient.For this satellite date of passage (28 July 2015), the NDVI (index of vegetation that varies in a scale from -1 to 1) negative values correspond to water surfaces and positive values represent surfaces, low values for less vegetated and high values for more vegetated) displayed a maximum value of 0.87 and a minimum value of -0.88.
The SAVI varied between -0.42 and 0.83, and the LAI was between 0 and 6.Surface temperature exhibited a maximum value of 39.05 °C and a minimum value of 5.5 °C (average maximum surface temperature in July in the study site is 29.1 °C).The surface albedo was between 0.03 and 0.94 and the radiation balance varied from 34.41 W m -2 to 598 W m -2 .
Figure 3 shows examples of the spatial distribution ET annual and the K c , estimated using the SEBAL algorithm.The annual evapotranspiration ranged from 0.0 to 960.30 mm, and vegetation areas displayed the highest values of ET annual.Areas with bare soil were showed the lowest values of ET annual and areas with rivers and lakes showed average values of the variable.Regarding the crop coefficient, a variation between 0.0 and 1.20 was observed and, as expected, the vegetation areas had the highest values, while areas with bare soil, had the lowest.The central pivot areas with adult coffee plants can be easily identified since they retain high K c values.The areas with adult coffee crops present a higher percentage of vegetated surface when compared to the areas with young coffee plantations planted with 3.5 m spacing between rows in the study site.Furthermore, the SEBAL method is calibrated based on the "hot" and "cold" pixels that are selected according to the vegetal cover percentage.Thus, ET c estimated by SEBAL is higher in areas with higher vegetated surface, if the pixel is chosen properly, resulting in higher K c values.
Table 3 shows the monthly average values of ET a and K c of young and adult coffee plantations based on the SEBAL algorithm.For young coffee crops, the ET a was 2.17 mm d -1 , on average, with a standard deviation 0.7 mm d -1 , on average.The mean value K c was 0.6 for coffee plants in the same phase, with a standard deviation 0.43, on average.
Regarding the adult coffee plantation, the ET a was 3.95 mm d -1 , on average, with a standard deviation 0.78 mm d -1 , on average.The low standard deviation values observed among the pixels are indicative of a uniform water supply by the central pivots, which enables homogeneous growth of coffee plants, resulting in a smaller difference between the K c values at different points of the plantation.The mean K c value was 0.85 for coffee plants in the same phase, with a mean standard deviation 0.48.On average, the K c for adult coffee plants was higher in the final phase and lower in the initial phase, with median value during the intermediate phase.
Figure 4 (A, B, C, D and E) show values of ET c of young coffee plantations based on the SEBAL algorithm and ET c calculated with observed data.Comparing ET c obtained by the SEBAL algorithm with ET c calculated and recommended by Lima and Silva (2008) using the variety Rubi MG 1192 (Figure 4A) and Acaiá Cerrado (Figure 4B), the Index of Agreement (d) described exhibited values below 0.44, regardless of the variety, indicating low accuracy of the model in these situations.The Root Mean Square Error (RMSE) was below 0.94, the Mean Bias Error (MBE) was below 0.82 and the Normalized Root Mean Square Error (NRMSE) was below 0.75, indicating a rather unacceptable precision of the ET c estimate, based on the SEBAL method.
This may have occurred because crop coefficients show little or no variation throughout the year, recommended by the reference studies.This issue differs from the other estimation studies, carried out using the SEBAL model, which obtained a different crop coefficient value for each month.Furthermore, sites of young coffee plantations present a higher rate of soil without vegetation cover, reducing accuracy of the estimation model.For future studies in these conditions, changes are suggested in the algorithm including correction factors for estimates in young coffee plantations.
In the validation of the ET c obtained by the SEBAL algorithm with the values of the ET c calculated by observed data and recommended by Santinato et al. (2008) using the crop density of 2500 (Figure 4C), 3300 (Figure 4D) and 6700 plants ha −1 (Figure 4E), the index d described exhibited values lower than 0.50, regardless of the crop density, indicating low accuracy of the model in these situations.The RMSE was below 0.79, the MBE was below 0.67 and the NRMSE was below 0.63, also indicating unacceptable precision of the ET c estimate, based on the SEBAL method.
The comparison of ET c values of the adult coffee plantations based on the SEBAL method with values of ET c calculated is shown in Figures 5A, B, C, D, E,  F, G, and H.In the validation of the ET c obtained by the SEBAL algorithm with the ET c calculated through observed data and recommended by Allen et al. (1998) with (Figure 5A) and without (Figure 5B) plant covering, index d presents values higher than 0.62, regardless of plant cover, indicating relatively high accuracy of the model in these situations.The RMSE was below 0.51, the MBE was below 0.43 and the NRMSE was below 0.50, also indicating reasonable precision of the ET c estimate, based on the SEBAL method.Rawat et al. (2017) estimated the ET c of wheat crop by the SEBAL and standardized FAO-PM methods in Bhiwani District of Haryana, India, and the results obtained were evaluated through statistical performance measure tests.The statistical parameters found (RMSE = 0.56, NRMSE = 0.20 and d = 0.87) showed a good agreement between SEBAL ET c and PM ET c .The findings of this work suggested that the SEBAL model shows a good potential to estimate spatial ET c for the region studied.
In the validation of the ET c obtained by the SEBAL algorithm with the values of the ET c calculated with observed data and recommended by Villa Nova et al. ( 2002) with (Figure 5C) and without (Figure 5D) plant covering and with crop density of 2500 (Figure 5E), 3300 (Figure 5F) and 4000 plants ha −1 (Figure 5G), the index d described exhibited values higher than 0.40, regardless of the plant covering and the crop density, indicating relatively high accuracy of the model in these situations.The RMSE was below 1.20, the MBE was below 1.12 and the NRMSE was below 1.17, also indicating reasonable precision of the ET c estimate, based on the SEBAL method.
In the comparison of the ET c obtained by the SEBAL algorithm with the values of the ET c calculated with observed data and recommended by Sato et al. (2007) (Figure 5H), the index d presents exhibited value 0.50, indicating low accuracy of the model in this situation.The RMSE was 0.77, the MBE was 0.69 and the NRMSE was 0.75, also indicating reasonable accuracy of the ET c estimate, based on the SEBAL method.Bala et al. (2016) also estimated the ET c of wheat crop by SEBAL and compared it to ET from lysimetric technique by using statistical parameters.The authors found RMSE values equal to 0.19 and NRMSE 0.21 mm d −1 .The errors or deviations found in the comparisons between estimated and measured values may come from the estimation model itself or from the way these data are compared.The comparison methodology using fixed K c values for long periods of the year can also result in unexpected values of certain statistical parameters.

Conclusion
The ET c and K c data obtained based on the SEBAL algorithm displayed similar values to studies that used traditional methods.The SEBAL algorithm is a more parsimonious alternative in K c estimation; therefore, the K c obtained by this method can be used directly in coffee crops.The method of ET a and spatialized K c estimation allowed understanding the variability of the study site, supporting the use of precision irrigation.Therefore, this model has high potential to estimate ET of different stages of coffee plantations for the region studied.This may also help understand crop water requirement and irrigation scheduling in the region.

Figure 1 -
Figure 1 -Location of the study site, with image obtained by the Operational Land Imager (OLI) sensor -Landsat 8 on 28 July 2015.Note the central pivots in the young and adult coffee crops and the weather station.

Figure 2 -
Figure 2 -Spatial distribution of the variables used in the Surface Energy Balance Algorithm for Land (SEBAL).Normalized Difference Vegetation Index (NDVI), Soil Adjusted Vegetation Index (SAVI), Leaf Area Index (LAI), Surface temperature (T s ), Surface albedo (a S ) and Radiation balance (R n ).Variables were calculated based on images obtained on 28 July 2015.

Figure 3 -
Figure 3 -Spatial distribution of annual evapotranspiration (ET annual) and crop coefficient (K c ) in the studied areas.K c was calculated based on images obtained on 28 July 2015 and ET annual was calculated with the average evapotranspiration in 2014 and 2015.

Figure 4 -
Figure 4 -Crop evapotranspiration (ET c ) of young coffee plantations based on the SEBAL algorithm compared to ET c (mm d -1 ) calculated with observed data.ET c observed was obtained by multiplying the reference evapotranspiration (ET o ) calculated by the Penman-Monteith (PM) method by the crop coefficient (K c ) of reference studies in Brazil.Lima and Silva (2008) using the variety Rubi MG 1192 (A) and Acaiá Cerrado (B) and Santinato et al. (2008) with crop density of 2500 (C), 3300 (D) and 6700 plants ha -1 (E).

Figure 5 -
Figure 5 -Crop evapotranspiration (ET c ) of adult coffee plantations based on the SEBAL algorithm compared to ET c (mm d -1 ) calculated with observed data.ET c observed was obtained by multiplying the reference evapotranspiration (ET o ) calculated by the Penman-Monteith (PM) method by the crop coefficient (K c ) of recommendations from the FAO-56 bulletin and reference studies in Brazil.Allen et al. (1998) with (A) and without (B) plant covering, Villa Nova et al. (2002) with (C) and without (D) plant covering and with crop density of 2500 (E), 3300 (F) and 4000 plants ha -1 (G) and Sato et al. (2007) (H).

Table 1 -
Satellite transits date referring to Operational Land Imager (OLI)-Thermal Infrared Sensor (TIRS) imagery used in the study. Sci.

Agric. v.76, n.2, p.93-101, March/April 2019
, using the REF-ET software (Reference Evapotranspiration Calculation, version 2.01.14).Solar radiation, relative humidity, wind speed and precipitation are required to estimate ET o in REF-ET and a dataset of maximum and minimum temperatures.

Table 2 -
Conditions for the development of coffee crop coefficient (K c ) studies and values used in comparison of the coffee K c estimated by the SEBAL model with the K c recommended by the reference studies.

Table 3 -
Monthly average values of actual evapotranspiration (ET a ) and crop coefficient (K c ) of young and adult coffee plantations based on the SEBAL algorithm and the respective standard deviation (SD).
*Monthly average in 2014 and 2015; **Mean values in the four pivots analyzed; ***Mean values in the 13 pivots analyzed; ET a = actual evapotranspiration (mm d -1 ); K c = crop coefficient; SD = standard deviation.