Seasonal Crop Water Balance Using Harmonized Landsat-8 and Sentinel-2 Time Series Data

: Efﬁcient water management in agriculture requires a precise estimate of evapotranspiration ( ET ). Although local measurements can be used to estimate surface energy balance components, these values cannot be extrapolated to large areas due to the heterogeneity and complexity of agriculture environment. This extrapolation can be done using satellite images that provide information in visible and thermal infrared region of the electromagnetic spectrum; however, most current satellite sensors do not provide this end, but they do include a set of spectral bands that allow the radiometric behavior of vegetation that is highly correlated with the ET . In this context, our working hypothesis states that it is possible to generate a strategy of integration and harmonization of the Normalized Difference Vegetation Index ( NDVI ) obtained from Landsat-8 ( L 8) and Sentinel-2 ( S 2) sensors in order to obtain an NDVI time series used to estimate ET through ﬁt equations speciﬁc to each crop type during an agricultural season (December 2017–March 2018). Based on the obtained results it was concluded that it is possible to estimate ET using an NDVI time series by integrating data from both sensors L 8 and S 2, which allowed to carry out an updated seasonal water balance over study site, improving the irrigation water management both at plot and water distribution system scale.


Introduction
Under the constant pressure of population growth, food consumption is increasing in nearly all regions of the world. The world population is expected to reach 9 billion by 2026 [1], which implies an increase in the worldwide area under irrigation of 300,000 km 2 , along with a 40% increase in water and energy demand over the next 20 years [2]. The sector with the greatest water use in Chile is agriculture, accounting for 78% of the total water availability [3], which supplies an irrigated area of 11,000 km 2 [4]. However, in various regions of the country, water use rights exceed the actual availability of water resources, which has led to numerous regions being declared depleted in terms of both surface and groundwater [5]. In this scenario, it is clear that one of the main challenges of the 21st century is to increase agricultural water productivity [6]. Thus, precise information on agricultural water demands is crucial for efficient management of water and crop productivity [7,8]. This high-water-demand scenario has triggered a search for solutions to narrow the gap between irrigation water demand and availability in terms of quantity and quality through the use of new technologies [9][10][11].
Water resource management strategies must rely on the estimation of crop water demand. One of the most widely used methods to determine water demand is the estimation of evapotranspiration (ET). Several studies have been done to evaluate the effect of water applied in crop yield [12][13][14][15][16] to optimize the water management in agriculture. Nonetheless, this process is difficult to correctly quantify when dealing with large areas, as there is great spatiotemporal variability due to the complex interactions between the soil, vegetation, and climate [17]. Currently, ET estimates are mainly based on observations from terrestrial weather stations. Several approaches have been proposed in the literature to estimate ET using weather station observations: (i) A two-step approach by multiplying the weather-based reference evapotranspiration ET r by crop coefficients (K c ) [18][19][20]. Crop coefficients are determined according to the in situ type of crop and the crop growth stage [19]. (ii) On the basis of the Penman-Monteith (P-M) equation [21], with crop to crop differences represented by the use of specific values of surface and aerodynamic resistances [22][23][24][25][26].
Presently it is possible to estimate ET for different crops, providing spatial and temporarily distributed information over a wide area, using information gathered from aircraft or satellite platforms. Two main strategies for ET estimation from remote sensing can be distinguished as follows: (i) Methods that use visible and near infrared sensors to extract a vegetation index (V I) and the radiative surface temperature to estimate the corresponding skin temperature [27,28]; (ii) Residual methods using the surface energy balance (SEB). These methods calculate ET by subtracting sensible heat and soil heat fluxes from net radiation [29]. Among the best-known methods that belong to this category are the surface energy balance index (SEBI) [30], the two-source model (TSM) [31], the evapotranspiration mapping algorithm (ETMA) [32], and the Mapping Evapotranspiration at high Resolution using Internalized Calibration (METRIC) model [33], which is based almost entirely on the Surface Energy Balance Algorithm for Land (SEBAL) model developed by Bastiaanssen et al. [34]. METRIC and SEBAL estimate crop ET by calculating the surface energy balance using spectral information from multispectral satellite images in the optical, near infrared, and thermal ranges. Only remote sensing imagery that provides spectral information in the thermal band may be used as input for these models. Unfortunately, most current satellite sensors do not provide this information, but they do include a set of spectral bands that allows the radiometric behavior of vegetation to be determined by focusing on the spectral contrast presented by plant cover in the red and near infrared bands [12]. Most V I are based on this principle; one which allows multitemporal data series to be constructed, thus providing essential information on water consumption patterns in various crop types and helping keep information on different agricultural covers up to date, as well as for monitoring of the biophysical properties of plants, such as plant cover, vigor, and growth dynamics [35]. In particular, a strong relationship between ET and the V I, known as the Normalized Vegetation Difference Index (NDV I), has been demonstrated [36]. In [6,37], the authors demonstrated that it is possible to estimate ET using NDV I using models trained previously with ET maps obtained by the METRIC model (considering visible, infrared, and thermal bands of Landsat 7 ETM+) and NDV I maps.
Images from the Landsat Program have generally been used for water consumption analysis and crop yield estimation. This is because of their high-spatial resolution in the visible spectrum (30 m) and the inclusion of the thermal bands (60 m), which allow for the quantification of the evapotranspiration in large areas at a temporal resolution of approximately 16 days in a precise and consistent way [38]. This acquisition frequency, however, can become a problem, especially in areas where climate conditions do not always allow good quality data to be available. As one might expect, a higher temporal frequency would be better. Sentinel-2 from European Space Agency (ESA) provides good spatial (20-60 m) and temporal (five days) resolutions for the visible spectrum, but lacks the Thermal Infrared (TIR) band [38]. Recently, methods have been developed with the aim of taking advantage of the characteristics of both products (Landsat-8 and Sentinel-2), with efforts to generate harmonized time series of surface reflectances for land monitoring applications being especially relevant [39][40][41]. The combination of the products of both programs allows for an effective increase in spatial and temporal coverage, providing a greater availability of data to users [42].
In this context, our working hypothesis states that it is possible to prepare a seasonal water balance in crops using serial data from Landsat-8 and Sentinel-2, obtained through the integration and harmonization between the NDV I of both sensors, known as NDV I . The expected result of this work is a NDV I time series used to estimate ET through specific adjustment equations for each type of crop, which allows one to continuously characterize the demand for water during an irrigation season.

Study Site
The study was carried out in one of the characteristic zones of the Central Valley of Chile. It has an area of 70 km 2 irrigated by the Convento Viejo reservoir, which has a water storage capacity of approximately 237 million m 3 . The area is located in the O'Higgins Region (252,626.83 E and 6,153,791.95 S, Zone 19, Datum WGS 84) along the South Lolol Canal (Figure 1). It is characterized by a temperate Mediterranean climate (with winter rain), a condition that favors the development of various crops such as forestry, orchards, cereal, and table grapes and vine, with 85.8% of the region's agroforestry area concentrated there [43]. The main agricultural land-use types identified in the study site were plums, olives, almonds, blueberries, table grapes and vines, industrial tomatoes, maize, wheat, cereals, and alfalfa. In this sense, for this work, the area was divided according to two types of cover: orchards (2.7 km 2 ) and annual crops (1.1 km 2 ) ( Figure 2).

Satellite Imagery
For this study, satellite images obtained from the Landsat-8 (L8) (L1T processing level, radiometrically calibrated and orthorectified) and Sentinel-2 (S2) satellites (1C processing level, with upper atmosphere reflectance values and orthorectified) were used. In Table 1, the characteristics of the spectral bands of each satellite are presented. For the image selection, the absence of clouds over the study area during the agricultural season (December 2017 to March 2018) was considered. The images selected are shown in Table 2.

Methods
The workflow used to achieve the research objective is shown in Figure 3. At the start, the spectral bands from the L8 and S2 satellites were atmospherically corrected using the Semi-Automatic Classification plugin in the open-source software QGIS 2.8.9. Then, the co-registration of the L8 and S2 images was carried out with the AutoSync tool in Erdas IMAGINE v.2011. This tool algebraically recognizes the coordinates of common points between the two satellite images. To carry out this process, the L8 satellite image was chosen as a reference, such that the other images were geometrically fitted to it. Subsequently, the subsampling of S2 from 10 m to 30 m pixel size was carried out using the linear interpolation method. To obtain the harmonization between the NDV I maps of the L8 and S2 sensors, the first step is to generate the NDV I maps during the season for both sensors, applying Equations (1) and (2), respectively; where NDV I L8 and NDV I S2 are the normalized difference vegetation indices for the L8 and S2 satellites, respectively.
For NDV I L8 and NDV I S2 harmonization, images from coinciding dates or with minimal temporal separation during the analyzed season were selected. For each pair of images, 618 sampling points in the annual crops and 638 in the orchard cover were selected. These points were chosen at random, and those out of the confidence interval were eliminated. The confidence interval is defined as absolute value difference between NDV I L8 and NDV I S2 for each pair of identified points. This value must be less or equal to the average value of difference. With the sampling points selected, a scatter plot for both agricultural covers was plotted, with the x-axis representing the NDV I S2 variable and y-axis the NDV I L8 variable. The fitting equation for the two variables was obtained from linear least squares regression, to which the 1:1 line represents perfect agreement between the NDV Is. In addition, for each graph, the Pearson correlation coefficient (r) was determined. The fitting equation that relates the harmonized vegetation index ( NDV I L8 ) is presented in Equation (3).
where i is the type of agricultural land-cover considered (orchards or annual crops), and a NDV I L8 i and b NDV I L8 i are the coefficients of linear fit between NDV I L8 i and NDV I S2 i for the analyzed season. Thus, the multimodal NDV I time series is composed of a total of three NDV I L8 and 14 NDV I L8 maps, and is called NDV I .
Following the workflow outlined in Figure 3, linear fitting between ET L8 and NDV I L8 for each L8 image ( Table 2) was carried out. The ET L8 maps at 30 m spatial resolution were generated from L8 images throughout the agricultural season using the SEB methodology proposed in Allen et al. [33] and implemented and validated by Fonseca-Luengo et al. [44]. The fit equation for five defined type of annual crops and six types of orchard covers, considering all L8 images is as follows: where k represents the different types of agricultural land use considered in Figure 2 (plums, olives, almonds, blueberries, table grapes and vines, industrial tomatoes, maize, wheat, cereals, and alfalfa), and j is the number L8 satellite images (DOY 4, DOY 52, and DOY 84, Table 2). Meanwhile, a ET L8 k,j and b ET L8 k,j are the coefficients of linear fit. Thus, for each map of ET L8 and NDV I L8 during the analyzed season, different fit equations for each orchard and annual crop type in the study area were obtained.
Finally, with Equation (5), it is possible to estimate evapotranspiration ( ET) for different agricultural land uses considered, by means of the NDV I time series.
where m represents the total amount of S2 images that were captured 30 days before and 30 day after the day j when the L8 image were captured.
To determine daily evapotranspiration during the season, it was necessary to obtain the actual crop factor (Fc) from the ratio between ET k,m and reference evapotranspiration (ET r ), estimated using the Penman-Monteith model [19] at the time the images were captured. In Equation (6), the Fc for each type of cover (k) and the whole series (m) is described: Thus, daily estimated evapotranspiration ET c was the product of ET r and the Fc obtained in each of the satellite images and extrapolated to the analyzed agricultural season.
Whereas, potential evapotranspiration (ET p ) was determined from the potential crop factor (Fcp), obtained from the literature [8,45], associated with the percentage of coverage during plant development. However, in our proposed method, this factor was adjusted to the NDV I values obtained from the satellite images during the season. This factor was called NDV I potential crop factor (Fcp NDV I ) and the way to obtain it was through the linear relationship between the Fcp and the NDV I of each crop. For this, the extreme values of Fcp were related to the NDV I values for each crop, while for the orchards, the difference was made between young trees and adult trees. The Fcp NDV I is described in Equation (7).
Thus, daily ET p was the product of ET r and adjusted Fc during the agricultural season for each crop in the study area.
To determine the water balance during the agricultural season, the estimated water demand and the potential demand of each of the annual crop and orchard types in the study area were determined. The two water demands are described in Equations (8) and (9), respectively.
where d is the start day and n the final day of the analyzed agricultural season. Thus, the seasonal water balance is composed of a comparison between estimated water demand and potential demand of each of the crops in the studied area.
The methodology evaluation was carried out by comparing the estimated water demand with the volume of water applied during the analyzed period. To this end, annual crops and orchards representative of the study site were selected. Meanwhile, information on the volumes applied was compiled from the records of the volume meters installed on the different farms. To compare the two variables, statistical indicators used root mean square error (RMSE), which measures the variation of the estimated values with respect to the observed values, and the bias indicator (BI AS), which provides information on the tendency of the model to over-or underestimate a variable.

Results
In Table 3, the image dates used to fit the harmonization between L8 and S2 are shown. Table 3. Image dates used to fit the harmonization between L8 and S2. The scatter graph between NDV I S2 and NDV I L8 for the three pairs of images for the orchards and annual crops are shown in Figure 4.
Using Equation (3), the harmonized relation for NDV I L8 and NDV I S2 to orchards is In addition, the harmonized relation for the annual crops is NDV I L8 = 0.8588 · NDV I S2 + 0.0997.
In Figure 4, it is possible to observe that the relationship between NDV I L8 and NDV I S2 is close to 1:1 for both cover types. In addition, the Pearson coefficient values (r) were 0.88 for the orchards and 0.94 for the annual crop cover. The difference observed between orchards and the annual crops could be due to the partial cover of the orchards, a topic that has been investigated by Souto et al. [46].
In addition, in Figure 4, it can be seen for both cover types that L8 tends to give large values of NDV I compared to S2, with values almost up to 0.6. Although Tello et al. [47], found a correlation coefficient (r) of 0.99 for the NDV I in agricultural crops, this relationship was determined without distinguishing between different cover types.  The available NDV I maps during the growing season are shown in Figure 5. The time series was composed by NDV I L8 and NDV I L8 , making a total of 17 images during the agricultural season. Using this multimodal NDV I time series, the ET was estimated for each orchard and annual crop in the study site.
The results show that there is a significant correlation between ET L8 and NDV I L8 in each of the orchards and annual crops on the three analyzed days. This is reflected by the Pearson correlation coefficient (r) values, which overall are above 0.8. Upon comparing this correlation to what appears in the literature, it can be seen that the r values are very similar in [48][49][50]; this being dependent on the type of cover under study. Thus, a total of 18 and 15 fitting equations for the orchards and annual crop agricultural covers, respectively, were obtained. Using each of these fitting equations, it was possible to estimate ET for the season which extends from December 2017 to March 2018, specifically for each type of orchard and annual crop type in the study area. Nonetheless, ET estimation could have been improved if more images from the season had been obtained for the ET L8 -NDV I L8 fitting, which in our study area was impeded by the presence of clouds.
By replacing the coefficients from Table 4 in Equation (3), the ET maps were obtained for each crop of the study site ( Figure 6). In Figure 6, it is possible to observe the temporal and spatial distribution of ET represented by monthly images, which was estimated using the different ET L8 -NDV I L8 fitting equations for each orchard and annual crop type in the study area.  The overall results show that ET presents the highest values on 20 December and 19 January, in the range of 5 and 6 mm·day −1 . Meanwhile, on 13 February and 10 March the values are between 4 and 5 mm·day −1 This lessening of ET starting on 13 February can be attributed to the senescence period, particularly for the annual crops, which present different phenological periods than the orchards.
ET of industrial tomatoes varies from 3 mm·day −1 in December to up to 6 mm·day −1 in January, a variation that coincides with the development of the crop as it reaches maturity in late February and early March. However, the difference in ET of plums is lower than that of industrial tomatoes, with values between 3 and 5 mm·day −1 .
In general, orchards presented a relatively homogeneous water demand over the irrigation season because the canopy cover does not change during the season. It is important to point out that ET is closely related to the cover percentage of orchards, as well as annual crops [8,45].
The selected crops are identified in Figure 7 to show the behavior of the daily ET c during the analyzed period, which were obtained from the ET maps that were generated during the agricultural season following the proposed methodology. The results shown in Figure 8 show the variability in water consumption through ET c of the different orchards and annual crops, showing that in general the analyzed orchards present a similar ET c behavior during the season. However, the behavior of table grapes presents higher values than those of the other orchards, reaching 7 mm·day −1 in December. Additionally, it is observed that during the first days of December orchards obtained the highest ET c values recorded during the season; attributed to better irrigation management during this period or possibly the contribution of soil evaporation rather than the transpiration of the trees themselves.
In the case of the annual crops, the behavior of the ET c is extremely variable among them, reflecting the natural differences in the water needs of each crop, with corn and the industrial tomato registering the highest values of ET c , in a range varying between 3 to 6 mm·day −1 . The global results of water demand for the crops estimated with our proposed methodology ( Figure 9) show that the applied volume of water to the crops was lower than the potential demand for all the season. This indicates that the water applied to different crops was insufficient or inadequate during the season. In the particular cases of olives, table grapes, and maize ( Figure 9, B, D and H), the estimated water demand was practically identical to the potential demand, showing that the amount of water applied during the season was the required one. This scenario is the least frequent in the studied area, even though the availability of the resource in the sector does not show restrictions. However, the rest of the crops showed an estimated water deficit, ranging from 6% to 48% for olive trees (B) and wheat (I), respectively.
The water deficit obtained by the different crops during the analyzed period can generate considerable damage to its yield. In the case of wheat, a water deficit causes a decrease in grain yield, especially when water scarcity occurs during the sensitive stages of growth (the flowering and grain filling stage) [51]. In blueberries, water deficit in any period of its phenological stage affects the vegetative growth, substantially reducing the potential yield, and in addition, it is a determining factor for the production in the following season [45,52]. Almond trees are considered drought tolerant, but irrigation scarcity is critical to the production of high yields with high quality nuts, mainly in the preharvest stage [53]. In plum trees, water stress in the final stages of fruit growth significantly decreases their size, but accelerates ripening, and increases the concentration of sugars [54]. For vines, water deficit in the stages after grapes veraison can improve the quality of the wine with almost no reduction in yield [55,56]. On the other hand, the imposition of severe drought stress in inappropriate phenological stages can cause a significant decrease in yield and even a decrease in quality in extreme cases [57]. As we have seen, the response of annual crops and orchards to a lack of water is different in each case and the decisions that must be taken in irrigation management must be appropriate for each one of them.
The results shown in Figure 10 show that the estimated volume of water applied for maize, plums, and olive trees (A, B, C, and D) satisfies the estimated water requirements through the potential demand. In the case of the adult plum tree (B), the water applied during the season is less than the amount of water required by the orchard. In addition, it is shown that the potential demand obtained for the adult plum tree was lower than that of the young plum, which can be attributed to the lack of water applied during the season or inadequate water management. Olive trees in general are resistant to a lack of water and show a high recovery capacity after prolonged periods of drought, but adequate irrigation management is required to produce economical yields. However, a certain degree of stress improves oil quality [58,59].
For industrial tomatoes (F), the applied volume exceeded the estimated volume by 1300 m 3 ·ha −1 , but this is almost identical to the potential volume. This shows that the plants only evapotranspire the amount of water that is able to be retained in the soil at the root zone and not the excess applied water. Excessive irrigation in tomato plants can cause excessive leaf growth, plants with high vegetative vigor tend to produce low quality fruit, and even in some varieties, large variations in moisture levels in the soil during fruit ripening can cause fruit cracking, spots, rot, and variation in size and shape [60].
Statistically, the total estimated water demand through the proposed methodology based on the volume applied had an RMSE of 0.6 mm·day −1 and a BIAS of −0.4 mm·day −1 .
Although the methodology implemented is able to effectively reproduce the amount of water applied, the considerable lack of water during the vital stages of plant development is demonstrated. Both water scarcity and inadequate applications can be combined with an uneven distribution of water in the field, often due to the low uniformity of irrigation systems. This especially occurs when irrigation systems are not properly designed, operated, and maintained. This can cause a substantial decrease in crop yield and economic income [61].

Conclusions
On the basis of the obtained results, we conclude that it is possible to construct an NDV I time series by integrating L8 and S2 data, allowing the estimation of ET during an agricultural season and suitably demonstrating the responses of crops to irrigation excess and shortage problems associated with water management. This is based on the Pearson correlation coefficient (r > 0.8) values that were obtained for the orchards and annual crop cover. In addition, a comparison between estimated water demand and information on potential requirements meant that the water balance was kept up to date, according to the water needs of different agricultural cover types. This can encourage better management of applied water levels and production in different stages of the irrigation season. Further research on this topic, focusing on potential water demands for orchards and annual crops, will allow for better and more precise application of the proposed method.