Modeling Actual Evapotranspiration with MSI-Sentinel Images and Machine Learning Algorithms

: The modernization of computational resources and application of artiﬁcial intelligence algorithms have led to advancements in studies regarding the evapotranspiration of crops by remote sensing. Therefore, this research proposed the application of machine learning algorithms to estimate the ET r F (Evapotranspiration Fraction) of sugar can crop using the METRIC (Mapping Evapotranspiration at High Resolution with Internalized Calibration) model with data from the Sentinel-2 satellites constellation. In order to achieve this goal, images from the MSI sensor (MultiSpectral Instrument) from the Sentinel-2 and the OLI (Operational Land Imager) and TIRS (Thermal Infrared Sensor) sensors from the Landsat-8 were acquired nearly at the same time between the years 2018 and 2020 for sugar cane crops. Images from OLI and TIR sensors were intended to calculate ET r F through METRIC (target variable), while for the MSI sensor images, the explanatory variables were extracted in two approaches, using 10 m (approach 1) and 20 m (approach 2) spatial resolution. The results showed that the algorithms were able to identify patterns in the MSI sensor data to predict the ET r F of the METRIC model. For approach 1, the best predictions were XgbLinear (R 2 = 0.80; RMSE = 0.15) and XgbTree (R 2 = 0.80; RMSE = 0.15). For approach 2, the algorithm that demonstrated superiority was the XgbLinear (R 2 = 0.91; RMSE = 0.10), respectively. Thus, it became evident that machine learning algorithms, when applied to the MSI sensor, were able to estimate the ET r F in a simpler way than the one that involves energy balance with the thermal band used in the METRIC model.


Introduction
Evapotranspiration is the phenomenon of transferring water from liquid to gas to the atmosphere from the evaporation of water from the soil and water bodies, as well as the transpiration of plants through the leaves [1]. According to Allen et al. [2], quantifying water consumption in large areas, particularly in extended irrigated agricultural areas, is of great relevance for planning and managing water resources, mitigating impacts on the water bodies' streamflow rate, establishing use rights and consumption regulation, and avoiding conflicts of water use.
Evapotranspiration can be determined directly in the field by lysimeter techniques, which are quite reliable but costly, or estimated with (1) full-physical models based on the principles of conservation of mass and energy; (2) semi-physical models that use conservation of mass or energy; and (3) black-box models based on artificial neural networks, By aiming to develop models with higher prediction capacity and reliability, this work, through the application of machine learning (ML) algorithms, sought to model variables of interest. Such algorithms, based on artificial intelligence, have a robust structure that allows the identification of relationship patterns between the variables to be modeled and the so-called predictor variables (independent). According to Cervantes et al. [12], ML is an interdisciplinary area based on computer science, statistics, mathematics, and optimization, among several other areas. Many ML algorithms, each with different characteristics, are used for prediction and classification, being, in recent years, applied in research focused on agricultural sciences and remote sensing to develop models with greater ability to represent phenomena [13][14][15][16][17]. In evapotranspiration prediction, these models used remote sensing data and climate data to estimate actual evapotranspiration [18,19], as well as only climate data to estimate reference evapotranspiration [20][21][22][23]. The various available studies applied the algorithms in predicting reference evapotranspiration or actual evapotranspiration using remote sensing with thermal data. Thus, this work sought to train machine learning algorithms to estimate the evapotranspiration fraction (ETrF) used in the METRIC model from the data of the Sentinel-2 satellites constellation that does not use thermal data on the sugar cane crop used in this case study.

Study Area
The study was performed in the northern region of the state of Minas Gerais, Brazil, classified as AS climate-tropical with dry summer [24] (Figure 1). The sugar cane crop is under 150 irrigated fields using center pivot irrigation used in this study. This area was chosen because it is a research partnership area with little to no precipitation ( Figure 2) that favors the acquisition of a larger number of cloud-free images, thus being able to acquire a volume of spectral information of the area, mainly coinciding with Landsat-8 and Sentinel-2. This area was chosen because it is a research partnership area with little to no precipitation ( Figure 2) that favors the acquisition of a larger number of cloud-free images, thus being able to acquire a volume of spectral information of the area, mainly coinciding with Landsat-8 and Sentinel-2.

Landsat-8 and Sentinel-2 Data
The data from the sensors onboard the Landsat-8 and Sentinel-2 satellites were acquired to generate the observed variable through the METRIC model, and the explanatory variables used for training and testing the machine learning algorithms can be seen in a simplified way in the flowchart in Figure 3.
A total of 56 satellite images were acquired, of which 22 were from OLI (Operational Land Imager) sensor and TIRS (Thermal Infrared Sensor) from the Landsat-8 satellite, and 34 came from the MSI (Multispectral Instrument) sensor of the Sentinel-2A and 2B satellites. Images from OLI and TIRS sensors were necessary to estimate the evapotranspiration fraction (ETrF) using the METRIC model, this being the response variable, whereas data from the MSI sensor were used as predictor variables in the machine learning models. The spectral and spatial characteristics of the sensors used are available in Tables 1 and 2.

Landsat-8 and Sentinel-2 Data
The data from the sensors onboard the Landsat-8 and Sentinel-2 satellites were acquired to generate the observed variable through the METRIC model, and the explanatory variables used for training and testing the machine learning algorithms can be seen in a simplified way in the flowchart in Figure 3.

Response Variable
The response variable, evapotranspiration fraction (ETrF), was obtained for each pixel, with a 30 × 30 m spatial resolution, from the METRIC model. ETrF is reckoned by dividing instantaneous evapotranspiration (ETinst) in each pixel by the hourly reference evapotranspiration (ETr) given by the meteorological station. Allen [7] standardized ETr  A total of 56 satellite images were acquired, of which 22 were from OLI (Operational Land Imager) sensor and TIRS (Thermal Infrared Sensor) from the Landsat-8 satellite, and 34 came from the MSI (Multispectral Instrument) sensor of the Sentinel-2A and 2B satellites. Images from OLI and TIRS sensors were necessary to estimate the evapotranspiration fraction (ET r F) using the METRIC model, this being the response variable, whereas data from the MSI sensor were used as predictor variables in the machine learning models. The spectral and spatial characteristics of the sensors used are available in Tables 1 and 2.  Images were obtained between the years 2018 and 2020. During this period, eight images from the Landsat-8 satellite and eight from the MSI sensor covered the study area on the same day. Images of those dates were used to train and test the models and the application of residual analysis. The choice for images with coinciding dates between the two satellites was made to avoid different reflectance of the same study area since vegetation, especially crops, change quickly. Thus, obtaining images on different days to perform the modeling could have introduced bias and not a faithful representation of the modeled product. In Table 3, images with date, time (GMT-3), and path/row used to train the models are available.  Table 3) and 30 from MSI sensors (including 4 images in Table 3). Those images were used to determine crop coefficient (K c ) and actual evapotranspiration (ET a ) during the sugar cane cycle. The estimation was performed by both the METRIC model and the suggested model.

Response Variable
The response variable, evapotranspiration fraction (ET r F), was obtained for each pixel, with a 30 m × 30 m spatial resolution, from the METRIC model. ET r F is reckoned by dividing instantaneous evapotranspiration (ET inst ) in each pixel by the hourly reference evapotranspiration (ET r ) given by the meteorological station. Allen [7] standardized ET r for alfalfa at 0.5 m high. According to the authors, when using this condition, the ET r F can be considered equivalent to the crop coefficient (K c ). It also enables the extrapolation of the crop's actual evapotranspiration when the satellite switches to the daily 24 h level. Thus, the ET r F is determined by Equation (1).
where ET inst is the instantaneous evapotranspiration (mm·h −1 ) and ET r is the reference evapotranspiration (mm·h −1 ) standardized to alfalfa at 0.5 m height at the moment the satellite is passing. ET inst is calculated from the latent energy consumed in the evapotranspiration (ET) process and the latent heat of vaporization Equation (2).
where 3600 is the conversion from seconds to hours of the duration of the satellite pass; LE is the latent energy consumed in ET (W m −2 ); ρ w is the density of water (~1000 kg m −3 ); λ is the latent heat of vaporization (J kg −1 ). λ represents the heat absorbed when one kg of water is evaporated, and it is determined by Equation where T s is the surface temperature ( • K) determined by band 10 of the TIRS sensor (Table 1). LE is calculated from the Earth's surface energy balance, Equation (4), which involves net radiation (R n ), a sensible flux of heat transferred to the ground (G), and a sensible flux of heat convected to air (H). These three components, responsible for determining LE, are usually expressed in W m −2 .
Atmosphere 2022, 13, 1518 7 of 23 R n is the radiant energy of the surface that is partitioned into H, G, and LE and is determined by Equation (5).
where R s↓ is the input of shortwave radiation (W m −2 ); α is the surface albedo (adim.) determined by bands 2, 3, 4, 5, 6, and 7 of the Landsat-8 (Table 1); R L↓ and R L↑ are the input and output of long waves (W m −2 ), respectively; and ε 0 is the surface thermal emissivity. G is the heat storage flux in the ground due to heat conduction. When vegetation is present, G tends to have lower values, which means that the rate of heat storage in the soil is lower. The METRIC model provides two methodologies to quantify G, one developed by Bastiaanssen [4] and another by Tasumi [25] (the details can be seen in Allen [7]). For this work, the methodology of Tasumi [25] was chosen, as it was developed in irrigated areas. Thus, G was quantified by Equation (6a,b).
where LAI is the leaf area index (m 2 m −2 ) estimated by the methodology applied by Allen [2]. Finally, H represents the rate of heat loss to air by convection and conduction influenced by a temperature difference. H was calculated by Equation (7).
where ρ air is the air density (kg m −3 ); C p is the specific heat of air at constant pressure (J kg −1 K −1 ); dT is the temperature difference between two heights (z1 and z2) in an area close to the surface; r_ah is the aerodynamic drag (s m −1 ) between these two heights. All calculations to reach ET r F described in this topic were performed using the Water package developed by Olmedo [26] for the R language environment [27].

Data Extraction for Training
Images from MSI sensors were divided into two groups according to their spatial resolution. Two approaches were applied: approach 1 using a 10 m spatial resolution and approach 2 with a 20 m spatial resolution. This division aimed to understand the effect of the number of predictor variables and spatial resolution on the performance of predictive models.
The primary variables selected for approach 1 were reflectance of blue (ρB), green (ρG), red (ρR), and near-infrared (ρNIR) bands; the variables for approach 2 were reflectance of blue (ρB), green (ρG), red (ρR), red edge 1 (ρRe1), red edge 2 (ρRe2), red edge 3 (ρRe3), near-infrared (ρNIR), shortwave infrared 1 (ρSWIR1), and shortwave infrared 2 (ρSWIR2). As approaches 1 and 2 have, respectively, a 10 and 20 m spatial resolution, it was necessary to match their dimensions with the dependent variable, that is, reduce their resolution to 30 m × 30 m. Therefore, the resample function was used with the bilinear method of the raster package [28] for the R language based on the images of the dependent variable. With images in the same spatial resolution, the central irrigation pivots were cut within each scene using a shapefile, and then the pivots were separated into two groups. The first group encompassed 70% of the total pivots and was used to train the model. The remaining 30% of center pivots were used to test the model ( Figure 1).
Due to computational limitations, it was not possible to use all pixels of each pivot of the seven images. Therefore, 25,000 pieces of information, totaling 175,000 pieces of data, were randomly extracted to train the model (Table 3). To test the model, 10,714 pieces of data were also randomly extracted from each image, totaling 74,998 pixels, maintaining the 30% data proportion for model training. Afterward, ET r F values lower than 0 that eventually contained some pixels generated by noise in some images were filtered out. In the end, data frame files with geographic coordinates, spectral bands, and ET r F were obtained.

Production and Selection of Predictor Variables
By aiming to increase the number of predictor variables, the NRPB (normalized ratio procedure between bands) technique was applied to the primary predictor variables in approaches 1 and 2. The NRPB normalized all primary variables, therefore increasing the explanatory variable number for the models (Equation (8)). Hence, the number of variables produced on each approach was higher than the primaries, especially in approach 2, which has the biggest variable number.
where ρ x and ρ y correspond to the surface reflectance of the wavelengths of the MSI sensor. The NRPB process was performed using the band ratio function of the labgeo package [29] for the R language, successfully applied by Filgueiras [14,30].
A two-step selection was applied to the NPRB variables to select the most important model variables. The first step removed correlated variables. Explanatory variables with a correlation above 95% were removed to decrease information redundancy. The second step removed the least important variables. Therefore, the recursive feature elimination algorithm (RFE) was applied through the caret package [31] in the R language. This algorithm removes the least important predictor variables from the model by building a base model with all predictors. Then, the importance of each predictor for the base model was calculated, and subsequently, the algorithm classifies these predictors and removes the least important ones, leaving a minor number of variables for the training steps. This is a key process, as it reduces unnecessary variables, aiding both training steps and future applications.

Training and Test
The machine learning model training consisted of the input of the remaining predictor variables after the variables selection process into four regression algorithms: Linear regression (LM); Cubist; eXtreme gradient boosting-linear method (Xgblinear); and eXtreme gradient boosting-tree method (Xgbtree). These models were chosen for having an elevated predictive potential and for performing fast on the training, which was carried out with the aid of the caret package [21]. After training, 8 possible predictions were obtained, 4 for approach 1 and 4 for approach 2, which were submitted to testing for evaluating the assertiveness and statistical error.
Thirty percent of data extracted from the pivots were used to test the trained algorithms. Statistical analyses for the test were: the coefficient of determination (R 2 ), Equation (9); root mean squared error (RMSE), Equation (10); mean absolute error (MAE), Equation (11); mean bias error (MBE), Equation (12).
Atmosphere 2022, 13, 1518 where Pi is the value predicted by the model; Oi is the observed value; O is the average observed value; n is the number of data pairs.

Residual Analysis
In addition to statistical metrics, a residual analysis was performed to evaluate the error behavior of predicted values for all models trained for both approaches. Therefore, an image was selected exclusively for this purpose (Table 3). In this image and inside the study area, 11 pivots were selected (highlighted in yellow in Figure 1) for residual analysis between the observed and the predicted values. The 11 pivots were chosen for being close to each other, making it easier to elaborate thematic maps, and because there was no cloud coverage near them. It is noteworthy that the image destined for residual analysis was not in the training process, so the data from this image are new information for the models already trained. Thus, models were applied. Then the predicted ET r F (approaches 1 and 2) and the observed ET r F (METRIC) were extracted from the 11 pivots ( Figure 1) to reckon the residues (Equation (13)).

Application
Among the selected pivots (Figure 1), one pivot (highlighted in black) was used for the application of the trained models. This pivot was chosen due to having the largest number of Sentinel-2 images during the crop cycle. It had images from before sugarcane budding to harvest. The application consisted in quantifying the crop coefficient (K c ) and the actual evapotranspiration (ET a ) during the sugarcane crop cycle for the best METRIC model. In addition to K c , water consumption during the sugar cane crop cycle also was estimated using 29 images of the MSI sensor from 11/23/2019 (DAE-days after emergency 05) to 09/03/2020 (DAE 290) and 17 images of the OLI and TRIS sensor between 12/07/2019 (DAE 19) and 09/04/2020 (DAE 291).
From the 29 images of the MSI sensor, the model with the best statistical results on approaches 1 and 2 was applied to determine the ET r F. According to Allen [7], the ET r F is equivalent to the K c when using the ET r of the alfalfa, which is the condition in which the ET r F was modeled. The 17 images of the OLI and TIRS sensor were used to calculate K c through the METRIC model. After quantifying the K c , the crop's actual evapotranspiration was established. Therefore, the accumulated reference evapotranspiration for a 24 h (ET r-24 ) period was calculated on the same day the image was acquired.
The ET r-24 in each day was calculated using the Penman-Monteith equation standardized by the ASCE (American Society of Civil Engineers) [32], and meteorological data were acquired in the A539 station, Mocambinho from the INMET (National Institute of Meteorology). With all the ET r-24 and K c , the ET a in each pixel was established by Equation (14).
Finally, the sugar cane crop total evapotranspiration was determined during the cycle through METRIC and approaches 1 and 2 applying integral function through time, where x was the dates referring to images, and y was the ET a of each date. Thus, the area under the curve, which corresponds to the total evapotranspiration during the cycle, was calculated. This process was performed using function auc from the MESS package [33] for the R language.

Models
In approach 1, 11 predictors were generated by normalizing the spectral bands of spatial resolution of 10 m. Only seven predictors remained after removing predictors with a correlation above 95% and five when using the recursive feature elimination (RFE). Figure 4 shows the results after the application of the RFE algorithm. Notably, the potential number of predictor variables (ideas) is five. Thus, in the case of more than five variables, the ones with low importance can be removed without significantly impacting the trained models. On the other hand, less than five variables negatively affect the model prediction.
where x was the dates referring to images, and y was the ETa of each date. Thus, the area under the curve, which corresponds to the total evapotranspiration during the cycle, was calculated. This process was performed using function auc from the MESS package [33] for the R language.

Models
In approach 1, 11 predictors were generated by normalizing the spectral bands of spatial resolution of 10 m. Only seven predictors remained after removing predictors with a correlation above 95% and five when using the recursive feature elimination (RFE). Figure 4 shows the results after the application of the RFE algorithm. Notably, the potential number of predictor variables (ideas) is five. Thus, in the case of more than five variables, the ones with low importance can be removed without significantly impacting the trained models. On the other hand, less than five variables negatively affect the model prediction. In approach 2, the number of predictors generated by normalizing the bands is 45, higher than in approach 1 because of the bigger number of spectral bands for the special resolution of 20 m. When establishing the baseline of 95% of correlation, 22 variables remained and were later reduced to 12 after the application of the RFE. Figure 5 shows the stability of the statistical metrics using 12 predictable variables to train the model. In approach 2, the number of predictors generated by normalizing the bands is 45, higher than in approach 1 because of the bigger number of spectral bands for the special resolution of 20 m. When establishing the baseline of 95% of correlation, 22 variables remained and were later reduced to 12 after the application of the RFE. Figure 5 shows the stability of the statistical metrics using 12 predictable variables to train the model. In Table 4, the 5 predictor variables selected for approach 1 and the 12 selected for approach 2 can be found. It is worth mentioning that the predictors in this table are not rated by importance, as each algorithm assigns different importance after the training. In In Table 4, the 5 predictor variables selected for approach 1 and the 12 selected for approach 2 can be found. It is worth mentioning that the predictors in this table are not rated by importance, as each algorithm assigns different importance after the training. In supplementary materials are the figures with the importance of the predictors for each model in approaches 1 and 2. Table 4. Selected predictor variables for training models in approaches 1 and 2.
For approach 1, the most common wavelengths among the predictors were B2 (blue) and B4 (red), both appearing twice, counted individually and in the NRPB, followed by B3 (green) and finally B8 (near infrared). B2 and B4, which stood out in this approach, are absorbed in higher intensity by plants from photosynthesizing pigments [11]. B3 is also absorbed by the same pigments, however, at a much lower intensity. While B8 tends to be reflected or transmitted because when there is absorption, the physiological structures of the plant heat up, and the plant tends to dissipate this energy through its structures. Therefore, the more nourished the dissipation structures are, the more efficient the reflectance and/or transmittance are. Thus, for approach 1, the wavelengths B2, B3, and B4 correspond to the pigments of the leaves, while B8 (Near-infrared) is related to the structural part. As a result, the selected predictors are important for both the photosynthesis process and the structure of the cellular components of the plants studied. Hence, three wavelengths have characteristics related to the photosynthesizing pigments and one related to the structure. In approach 2, the bands that stood out were B12 (short infrared), appearing five times, and B2, appearing four times. B12 is a short wavelength related to water content, which is an object of study in this work, through the physical-biological process of evapotranspiration. B5 and B6 bands, both red edges, also stood out, representing a transition among the visible region and near-infrared, and B8 and B11, shortwave infrared. Thus, it can be noted in this approach that the selected predictors are related not only to the photosynthetic and structural parts but also to water content. Hence, it can be inferred that three wavelengths with photosynthetic characteristics were assigned, one structural, two in the transition from photosynthetic to structural, and two with water content.
Test performances for approach 1 for each machine learning model can be seen in Figure 6. red. Thus, it can be noted in this approach that the selected predictors are related not only to the photosynthetic and structural parts but also to water content. Hence, it can be inferred that three wavelengths with photosynthetic characteristics were assigned, one structural, two in the transition from photosynthetic to structural, and two with water content.
Test performances for approach 1 for each machine learning model can be seen in Figure 6. Among the four models chosen for this approach, XgbLinear and XgbTree obtained the best statistical results, with R 2 = 0.80 and RMSE = 0.15; followed by the Cubist with R 2 = 0.77 and RMSE = 0.15; lastly, and with the worst result, the multiple linear regression (Lm) with the minimum R 2 = 0.73 and the major RMSE = 0.17. Regardless of the model, the high dispersion shown on the test draws attention, especially in values lower than 0.6 ( Figure 6). Ke et al. [34], when estimating evapotranspiration using machine learning in Landsat-8 data and MODIS for a heterogeneous environment, noted that in areas with crops, there was a higher dispersion between predicted vs. observed values than in areas with forest, grazing, and bushes. Thus, results from these authors agree with the ones found here, as, in crop areas, there is a higher surface movement, being more dynamic both in terms of vegetation cover and in terms of soil moisture. Figure 7 shows the statistical results of the test for approach 2, which used only a 20 m spatial resolution as a source of predictor variables. Landsat-8 data and MODIS for a heterogeneous environment, noted that in areas with crops, there was a higher dispersion between predicted vs. observed values than in areas with forest, grazing, and bushes. Thus, results from these authors agree with the ones found here, as, in crop areas, there is a higher surface movement, being more dynamic both in terms of vegetation cover and in terms of soil moisture. Figure 7 shows the statistical results of the test for approach 2, which used only a 20 m spatial resolution as a source of predictor variables. It is observed in Figure 7 that Lm, as well as in approach 1, obtained the worst result, while Cubist, XgbLinear, and XgbTree were the ones that best explained the response variable, with lower RMSE (0.10) and higher R 2 . R 2 , as in approach 1, had different values among these last three models, being 0.91 for XgbLinear and 0.90 for Cubist and XgbTree. When comparing Figures 6 and 7, approach 2 obtained more satisfying results than approach 1, with models with R 2 higher than 0.88 and RMSE lower than 0.11, whereas It is observed in Figure 7 that Lm, as well as in approach 1, obtained the worst result, while Cubist, XgbLinear, and XgbTree were the ones that best explained the response variable, with lower RMSE (0.10) and higher R 2 . R 2 , as in approach 1, had different values among these last three models, being 0.91 for XgbLinear and 0.90 for Cubist and XgbTree. When comparing Figures 6 and 7, approach 2 obtained more satisfying results than approach 1, with models with R 2 higher than 0.88 and RMSE lower than 0.11, whereas approach 1 had R 2 lower than 0.80 and RMSE higher than 0.15. It seems that approach 2 even reduced the dispersion seen in Figure 6 significantly.

Residual Analysis
The elaborated models were applied to the image of 02/12/2020 (Table 3), which was destined for the residual analysis. The results are shown in Figure 8. approach 1 had R 2 lower than 0.80 and RMSE higher than 0.15. It seems that approach 2 even reduced the dispersion seen in Figure 6 significantly.

Residual Analysis
The elaborated models were applied to the image of 02/12/2020 (Table 3), which was destined for the residual analysis. The results are shown in Figure 8.  Notably, Figure 8 maintains the same pattern found in Figures 6 and 7, in which approach 1 presented a result inferior to approach 2 due to the minor concentration of residual points close to the zero value. For Approach 2, a good performance was observed, as less than 10.3% of the values distanced from −0.2 to 0.2, which indicates a low prediction error, mainly XbgLinear. In both approaches, the explanatory variable tends to be less precise when ET r F is low, especially for values under 0.6, as previously discussed; this is better shown in Figure 9.
Notably, Figure 8 maintains the same pattern found in Figures 6 and 7, in which approach 1 presented a result inferior to approach 2 due to the minor concentration of residual points close to the zero value. For Approach 2, a good performance was observed, as less than 10.3% of the values distanced from −0.2 to 0.2, which indicates a low prediction error, mainly XbgLinear. In both approaches, the explanatory variable tends to be less precise when ETrF is low, especially for values under 0.6, as previously discussed; this is better shown in Figure 9.  Figure 9 shows the lower precision during the prediction of values lower than 0.6 in greater detail, with particular reference to pivot 9. On this pivot, when METRIC showed values close to 0.6, approaches 1 and 2 obtained predicted values close to zero, which can be seen in Figure 8 with the residual analysis. Nevertheless, Figure 9 also reveals the great assertiveness of the models, which can be seen especially in pivots 7 and 10, where prediction models were able to capture a wealth of details regarding METRIC estimates. On pivot 10, the southeast region shows low ETrF values through METRIC and machine learning models. However, for that pivot, this result was more prominent when using Sentinel-2 images, as they showed greater spatial resolution. Additionally, the thermal Figure 9. Spatial variability of ET r F for approaches 1, 2, and Metric for 02/12/2020. Figure 9 shows the lower precision during the prediction of values lower than 0.6 in greater detail, with particular reference to pivot 9. On this pivot, when METRIC showed values close to 0.6, approaches 1 and 2 obtained predicted values close to zero, which can be seen in Figure 8 with the residual analysis. Nevertheless, Figure 9 also reveals the great assertiveness of the models, which can be seen especially in pivots 7 and 10, where prediction models were able to capture a wealth of details regarding METRIC estimates. On pivot 10, the southeast region shows low ET r F values through METRIC and machine learning models. However, for that pivot, this result was more prominent when using Sentinel-2 images, as they showed greater spatial resolution. Additionally, the thermal band on Landsat-8, which has a coarser resolution (100 m), can mask nuances of the surface of the monitored area. Figure 10 shows the distribution of the crop coefficient (K c ) during the sugar cane cycle, determined through the XgbLinear, which presented better metrics for approaches 1 and 2 through the METRIC model and also through the FAO-56 report [1].

Models Application
band on Landsat-8, which has a coarser resolution (100 m), can mask nuances of the surface of the monitored area. Figure 10 shows the distribution of the crop coefficient (Kc) during the sugar cane cycle, determined through the XgbLinear, which presented better metrics for approaches 1 and 2 through the METRIC model and also through the FAO-56 report [1].  Figure 10 evidences what was previously discussed; for the initial phases of the crop, with low-density plant coverage, machine learning prediction model values showed uncertainties when compared to METRIC, but, as the culture develops, this error tends to reduce. For phase I (initial), when there is a crop emergency, the Kc for approaches 1 and 2 were 0.3 and 0.32, respectively, while METRIC estimated a Kc of 0.20, and the Kc determined by the 56 FAO`s report [1] is 0.40. In phase II (development), the Kc curves in approaches 1 and 2 become closer to the ones of the METRIC as the phenological state progresses, as well as the Kc-FAO curve becomes closer to both approaches and distances a little from the METRIC curve. For phase III (mid-season) crop growth, the average Kc was 0.98 for approach 1 and 1.00 for METRIC, while for approach 2, it was 1.02, and FAO's 56 report recommends using a 1.25 value. Moreover, in phase III, the METRIC curve behaves similarly to approach 2, which corroborates the results found in the models' tests. In phase IV (late-season), maturation, approach 1 was the most distant from both FAO's and MET-RIC's Kc, with a value equal to 0.69, while approach 2 had a Kc of 0.76 and METRIC of 0.77. Dingre and Gorantiwar [35], when quantifying Kc for sugar cane through the water balance method, found medium values for phases I, III, and IV of 0.36, 1.20, and 0.78, respectively. Silva et al. [36], for the Brazilian semiarid region, obtained a Kc for ratoon cane in phases I, III, and IV equal to 0.65, 1.10, and 0.85, different from FAO's recommendation. This evinces a divergence in Kc values shown in the literature from the ones recommended, which can be explained by the specificity of the local weather in which each study was performed, as well as the physiological conditions of the crop at the time of the satellite imaging. Phase III, for this work, reached the maximum value of 1.09; however, it showed variations throughout the whole phase.

Models Application
Regarding Figure 10, it is noticeable that with Sentinel-2, due to its better temporal resolution, it was possible to obtain a larger number of Kc information throughout the crop cycle than Landsat-8, which is important to obtain a major temporal variability of  Figure 10 evidences what was previously discussed; for the initial phases of the crop, with low-density plant coverage, machine learning prediction model values showed uncertainties when compared to METRIC, but, as the culture develops, this error tends to reduce. For phase I (initial), when there is a crop emergency, the K c for approaches 1 and 2 were 0.3 and 0.32, respectively, while METRIC estimated a K c of 0.20, and the K c determined by the 56 FAO's report [1] is 0.40. In phase II (development), the K c curves in approaches 1 and 2 become closer to the ones of the METRIC as the phenological state progresses, as well as the K c -FAO curve becomes closer to both approaches and distances a little from the METRIC curve. For phase III (mid-season) crop growth, the average Kc was 0.98 for approach 1 and 1.00 for METRIC, while for approach 2, it was 1.02, and FAO's 56 report recommends using a 1.25 value. Moreover, in phase III, the METRIC curve behaves similarly to approach 2, which corroborates the results found in the models' tests. In phase IV (late-season), maturation, approach 1 was the most distant from both FAO's and METRIC's K c , with a value equal to 0.69, while approach 2 had a K c of 0.76 and METRIC of 0.77. Dingre and Gorantiwar [35], when quantifying K c for sugar cane through the water balance method, found medium values for phases I, III, and IV of 0.36, 1.20, and 0.78, respectively. Silva et al. [36], for the Brazilian semiarid region, obtained a K c for ratoon cane in phases I, III, and IV equal to 0.65, 1.10, and 0.85, different from FAO's recommendation. This evinces a divergence in K c values shown in the literature from the ones recommended, which can be explained by the specificity of the local weather in which each study was performed, as well as the physiological conditions of the crop at the time of the satellite imaging. Phase III, for this work, reached the maximum value of 1.09; however, it showed variations throughout the whole phase.
Regarding Figure 10, it is noticeable that with Sentinel-2, due to its better temporal resolution, it was possible to obtain a larger number of K c information throughout the crop cycle than Landsat-8, which is important to obtain a major temporal variability of this coefficient. This can lead to more assertive crop management than if only Landsat-8 data were used. Saleem and Awange [37] mention that Sentinel-2 represents a new age for obtaining more precise information about the Earth's surface, as it has a greater spatial and temporal resolution among satellites that provide images for free. Nevertheless, information referring to K c can be expanded when joining the two orbital platforms, as more information will be obtained with a larger frequency of images, as highlighted in Filgueiras et al. [9].
Crop actual evapotranspiration in a spatial resolution of 10, 20, and 30 m for approach 1, approach 2, and METRIC, are shown in Figure 11. this coefficient. This can lead to more assertive crop management than if only Landsat-8 data were used. Saleem and Awange [37] mention that Sentinel-2 represents a new age for obtaining more precise information about the Earth's surface, as it has a greater spatial and temporal resolution among satellites that provide images for free. Nevertheless, information referring to Kc can be expanded when joining the two orbital platforms, as more information will be obtained with a larger frequency of images, as highlighted in Filgueiras et al. [9].
Crop actual evapotranspiration in a spatial resolution of 10, 20, and 30 m for approach 1, approach 2, and METRIC, are shown in Figure 11. Figure 11. Sugar cane temporal-spatial actual evapotranspiration in three different spatial resolutions evincing, through rectangles, coincident dates between Sentinel-2 and Landsat-8.
In Figure 11, attention is drawn to DAE 030 with a ray with an ET a larger than for the rest of the pivot. On that day, crops were in the initial emergency phase, and soil was exposed in most of the area, receiving water from irrigation. This ray corresponds to the moisture in the exposed soil, and it is displaced clockwise from the irrigation equipment. Such information is only visible with 10 and 20 m resolution, and thus, the METRIC model, due to the use of thermal images, can not show. Spatially detailed information, as seen in the Sentinel-2 images of Figure 11, has great value for field professionals. Coinciding dates between Sentinel-2 and Landsat-8 show high similarities in the spatial ET a between approaches 1 and 2 with METRIC. Approach 2 stands out for having larger spatial proximity. Such proximity is evidenced in Table 5, in which the averages ET as in approach 2 have the smallest differences in the estimated averages by the METRIC model. In addition, the standard deviation for approach 2 was also smaller than approach 1 but larger than the METRIC model. This fact might be attached to the more detailed spatial resolutions of approaches 1 and 2 when compared to the method using the METRIC model. After calculating the integral of ET a , it was found that the sugar cane total water demand along the crop cycle was 1417.77 mm for approach 1, 1474.26 mm for approach 2, and 1544.11 mm for METRIC. It is perceptible that total evapotranspiration estimated for approaches 1 and 2 were close to the total evapotranspiration by METRIC, where percentual differences were 8.18 for approach 1 and 4.52 for approach 2 when compared to METRIC. Approach 2 had a closer value to the one estimated by METRIC, going against the values found for dates indicated in Figure 11. Sugar cane total evapotranspiration during its cycle, found by Dingre and Gorantiwar [35], was 1388 mm. The one found by Silva et al. [38] for the Brazilian Northeast conditions was 1600 mm. Thus, it can be noticed that the total evapotranspiration found in this work is close to the values found in Brazil.

Discussions
The spectral bands of blue and red stand out in approach 1 because these wavelengths are absorbed mostly by alpha and beta-chlorophylls, resulting in the release of oxygen from photosynthesis [11,39,40]. It means that the magnitude of the release of this chalcogen occurs in them. Thus, the greater the absorption of blue and red wavelengths, the greater the oxygen release and, consequently, the greater the water vapor release since water vapor is released along with oxygen during photosynthesis, explaining why these spectral bands stood out. The green wavelength of the MSI sensor ranges from 0.53 to 0.59 µm, and the wavelengths in this range are less absorbed by pigments [11]; thus, it impacts little on photosynthesis and water release by plants. The near-infrared is reflected and/or transmitted by the structure of the mesophyll as a form of protection of its physiological and molecular structure. Thus, this length is not related to photosynthesis but to the physiological structure that influences photosynthesis indirectly. As a result, this wavelength has less prominence than the others. The wavelengths in the short infrared are the ones that stood out in the second approach because they are spectral bands strongly absorbed by the water content present in plants [41,42]. Therefore, when the leaves are in hydric comfort, these wavelengths tend to be more absorbed. On the other hand, when occurring hydric stress, these wavelengths tend to be less absorbed and more reflected. Several studies demonstrated the importance of shortwaves for monitoring and quantifying water stress on plants [43][44][45][46]. The red edge wavelengths are also important bands in approach 2 for the same principles as the short infrared. They are absorbed by water content in plants as a result of the sensitivity of red edge spectral responses to water content in plant leaves [47]. These authors mention that vegetation indices that have red edges in their equation have a higher sensitivity to water content in maize crops. Santos et al. [42], studying water stress and spectral response in sugarcane crops, identified that the near-infrared wavelength was sensitive to the water content in leaves. Chandel et al. [43] also obtained similar results with wheat crops. Thus, it strengthens the evidence for the importance of these wavelengths in quantifying water content, even more so when they are normalized with other bands through the NRPB. This is one of the arguments stating the superior predictive capacity of approach 2 compared to approach 1. The presence of the short infrared and red edge spectral bands in approach 2 significantly improves the statistical metrics of the models.
The dispersion found in ET r F values lower than 0.6 in both approaches may be related to greater heterogeneity of the cultivated area. When ET r F values are lower than 0.6, the crop does not completely cover the soil, and there may be patches of soil with different colors, moisture, etc., a proliferation of weeds in parts of the area, and even a greater degree of the vigor of the crop in parts of the plot. The fact that approach 1 is more dispersed than approach 2 is linked, quantitatively and qualitatively, to their predictors. Approach 2 had a total of 12 predictors: the same as approach 1 and some others, mainly the short infrared wavelength and red edge.
The main results and discussion found in the literature when applying METRIC to predict the evapotranspiration fraction is to consider ET r F to be equal to K c when using the reference evapotranspiration of alfalfa [7,[48][49][50]. However, this interpretation is not correct because sensors on board satellites or other platforms capture the information that is occurring in the field under natural conditions, and one cannot be sure that the plant is in maximum water comfort and with nutritional, pest, disease, and weed management adequate for maximum water uptake. Thus both METRIC and the models trained in this study estimate the product between K c and K s (stress coefficient), the latter being responsible for reducing ET r F to values lower than K c when the crop is under stress or making it equal to K c when the plant is in favorable conditions for maximum water uptake.

Conclusions
Approach 2 had statistical results superior to approach 1, mainly for the XgbLinear model, which obtained R 2 of 0.91 and RMSE of 0.10, whereas the metrics for the same model considering approach 1 were 0.80 and 0.15 for the R 2 and RMSE, respectively. This result was mainly influenced by the greater number of spectral bands that are strongly related to water content as the short infrared and the red edge, thus being the model to be applied. However, all models developed showed limitations when the dependent variable presents values lower than 0.6, a condition in which the crop canopy has not completely covered the soil, and there is greater variability.
Despite the limitation when the soil is not fully covered, the application proved efficient in predicting ET r F since the analysis had values close to zero and the maximum distance only in the linear regression algorithms for both approach 1 and approach 2. Furthermore, it was possible to note that ET r F cannot be considered the same as K c because the onboard sensors capture the actual condition of the crop in the field, hydric comfort or not. Hence, ET r F can be understood as the product between K c and K s that best represents field conditions. The use of the MSI sensor and machine learning techniques proved to be a new and simple alternative to estimate ET r F through spectral information, complementing the METRIC model estimation using OLI and TIRS sensors and increasing the frequency that information is generated for areas of interest. The combination of these sensors is useful to obtain the highest temporal resolution of the crop, especially for irrigated agriculture, which requires K c and K s coefficients to be determined daily for adequate replenishment of the irrigation blade.
Finally, this work highlights the possibility of using this methodology for other remote sensors to calculate spatial and temporal evapotranspiration, enriching the scientific debate.
Further, it can be applied in other locations and different cultures since the data needed are from the orbital sensors and meteorological data from the area of interest. Thus, the methodology, despite being extensive, is easy to apply with a low financial cost. However, in hopes of advancements, this methodology can be performed with field instrumentation to collect the actual evapotranspiration, as this evapotranspiration would be replaced by the evapotranspiration of the METRIC model, but more funding would be required.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/atmos13091518/s1, Figure S1: Importance of the variables for the approach 1 algorithm; Figure S2: Importance of the variables for the approach 1 Cubist algorithm; Figure S3: Importance of the variables for the approach 1 XgbLinear algorithm; Figure S4: Importance of the variables for the approach 1 XgbTree algorithm; Figure S5: Importance of the variables for the approach 2 LM algorithm; Figure S6: Importance of the variables for the approach 2 Cubist algorithm; Figure S7: Importance of the variables for the approach 2 XgbLinear algorithm; Figure S8: Importance of the variables for the approach 2 XgbTree algorithm.

Conflicts of Interest:
The authors declare no conflict of interest.