Investigating the roles of different extracted parameters from satellite images in improving the accuracy of daily reference evapotranspiration estimation

Agricultural water management, crop modeling, and irrigation scheduling are all dependent on the accurate estimation of reference evapotranspiration (ET0). A satellite image can also compensate for the lack of reliable weather information. So, in this study, stochastic gradient descent (SGD) has been implemented for optimizing multilayer perceptron (MLP) and developing SGD-MLP to estimate daily ET0 in Tabriz (semi-arid climate) and Babolsar (humid climate) stations, Iran, using extracted data from satellite images. The estimated ET0 values were compared to the determined ET0 based on the FAO-Penman–Monteith equation. Based on satellite image data collected from 2003 to 2021, the database was constructed. During the development of the abovementioned models, data from 2003 to 2016 (70%) were used for training purposes, and residual data (30%) were used for testing purposes. Additionally, the input variables, including land surface temperature (LST) day and night, normalized difference vegetation index (NDVI), leaf area index (LAI), and a fraction of photosynthetically active radiation (FPAR) from MODIS sensor, were utilized to estimate the daily ET0. Thus, there are three studied models; first is based on the LST, second on the vegetation indices, and third on the combination of the LST and the vegetation indices. Additionally, four performance indexes, including the coefficient of determination (R2), the root-mean-square error (RMSE), Willmott’s index of agreement (WI), and Nash–Sutcliffe efficiency, were utilized in order to measure the implemented model’s accuracy. According to the obtained results, the SGD-MLP-3 with input parameters of LSTday&night, LSTmean, LAI, NDVI, and FPAR gave the most accurate results with RMSE and WI values of as 0.417 mm/day, 0.973, for Tabriz and 0.754 mm/day, 0.922 for Babolsar stations, respectively. Conclusively, LST of daytime, nighttime, and average may be suggested as the most influential parameter for ET0 estimation.


Introduction
Water resources management, especially irrigation agriculture, relies heavily on reference evapotranspiration (ET 0 ) (Valipour 2015). In irrigation, ET 0 has the most direct and significant application. To calculate crop evapotranspiration in the field, ET 0 is the key parameter. There have been several studies showing that an irrigation system that is both detailed and rational can increase irrigation water efficiency (Cifre et al. 2005;Du et al. 2010). So, understanding the relationship between climatic conditions and evapotranspiration is fundamental for determining an agricultural crop's optimal water requirement (Maeda et al. 2011). Therefore, an effective irrigation system is based on the accuracy of evapotranspiration estimation. Several approaches to calculating ET 0 have been developed in recent decades using physical principles, such as mass transfer and energy balance (Maeda et al. 2011). However, the Penman-Monteith FAO-56 (FAO-PM) approach is the most commonly used method for calculating ET 0 (Djaman et al. 2015). Although FAO-PM is achievable, its implementation remains inconvenient due to its requirement for large amounts of meteorological data, which originate from standard meteorological observation stations (Hobbins 2016;Mokari et al. 2021). Therefore, it is necessary to collect information on temperature, wind speed, relative humidity, solar radiation, and other factors for implementing the FAO-PM method. Thus, a model with fewer input climatic data is highly desirable, especially in the absence of complete climatic parameters. Eddy covariance devices or lysimeters can be used to determine ET 0 (Vickers 2017). Nevertheless, field measurements of ET 0 via lysimeters and Eddy Covariance systems are inefficient, costly, and time-consuming. In addition, their effectiveness depends on the precise and accurate planning of experiments (López-Urrea et al. 2006). Consequently, remote sensing has been used and improved to estimate ET 0 at various spatial scales for decades (Borel 2003). There has been a significant amount of research concerning the estimation of ET 0 using limited data on the ground. These studies led to the development of several methods for estimating ET 0 , including Hargreaves and Samani, Priestley-Taylor, and Thornthwaite equations (Jabloun and Sahli 2008;Tomas-Burguera et al. 2017). These equations, based on limited meteorological data, can be used in the estimation of ET 0 as a result of reasonable correlation coefficients and root-mean-square errors. Due to their nonlinear, non-stationary, and stochastic nature, meteorological variables are dynamic. The complexity of evapotranspiration requires techniques beyond those available in empirical models. Therefore, to develop ET 0 models in a new pattern, artificial intelligence approaches were introduced. Several machine learning practices have been introduced in recent years to enable the accurate calculation of ET 0 by using a variety of approaches (Bai et al. 2021;Samadianfard and Panahi 2019;Shamshirband et al. 2020). The ET 0 model proposed by Falamarzi et al. (2014) was based solely on temperature and wind speed data and was successful in modeling ET 0 using an artificial and wavelet neural network in Australia. Using a limited set of meteorological data, Feng et al. (2016) demonstrated acceptable accuracy in ET 0 estimation using extreme learning machines (ELMs) and artificial neural networks optimized by genetic algorithms (GANNs). Additionally, a number of practical approaches have been developed that utilize remote sensing data rather than meteorological data to estimate ET 0 despite the lack of meteorological stations. Leaf area indices (LAI), land surface temperatures (LST), fraction of photosynthetically active radiation (FPAR), normalized difference vegetation index (NDVI), and vegetation types can be remotely sensed from satellites to provide continuous temporally and spatially continuous information on the biophysical properties of the land surface (Huete et al. 2002;Los et al. 2000;Myneni et al. 2002;Wan et al. 2002). Satellite remote sensing provides an effective method for measuring long-term surface temperature over wide areas compared to traditional ground-based measurements (Duan et al. 2017). The results of Maeda et al. (2011) showed a correlation coefficient of 0.67 using MODIS (Moderate Resolution Imaging Spectroradiometer) LST as a replacement for atmospheric temperature in temperaturebased ET 0 models. To the best of our knowledge, there have been few types of research evaluating the potential of machine learning models for estimating ET 0 based on remote sensing data in regions with a variety of climate zones. Most of the researches focuses on meteorological data. So, in the current study, (1) a comprehensive evaluation was conducted on two models, namely the optimization algorithm of stochastic gradient descent (SGD-MLP) and multilayer perceptron artificial neural network (MLP) implementing only remote sensing data for estimating daily ET 0 within two different climate zones.
(2) Various combinations of remote sensing data were used to evaluate proposed ET 0 estimation models. (3) Examining the relationship between input parameters obtained from satellite images and ET 0 was done in two different climates.

Study area
In the present study, two different agro-climatic zones, including Tabriz and Babolsar stations, were selected in Iran. The climatic characteristics and location specifications of the two mentioned stations are presented in Table 1. The location maps of the selected stations are also shown in Fig. 1

Overview
To estimate daily ET 0 in Tabriz and Babolsar stations, the following input variables were used: land surface temperature (LST) day and night, normalized difference vegetation index (NDVI), leaf area index (LAI), and a fraction of photosynthetically active radiation (FPAR) from MODIS sensor. So, data were collected in the time period of 2003-2021. In other words, ET 0 is computed using machine learning algorithms based on remote sensing data to identify the relationships between remote sensing data and ET 0 (calculated by meteorology station data). Figure 2 illustrates the structure of the model. A daily dataset from the MODIS sensor is required to estimate ET 0 in this study. This research utilized 4-day LAI, FPAR, and 16-day NDVI data, so they need to be reconstructed and converted into daily data. In this case, cubic spline interpolation has been applied. To reconstruct the MODIS LST day/night data and correct the missing data due to atmospheric disturbances and dust, Terra and Aqua satellites' 1-day data were combined first to reduce the number of missing data; then, the Kalman filter was used to reconstruct the 1-day data.

FAO-Penman-Monteith (FAO-PM)
A standard method for estimating ET 0 is the FAO-PM equation. FAO-PM was calculated using Eq. (1).
where ET 0 is the reference evapotranspiration (mm d −1 ), Δ is the slope of the saturation vapor pressure curve (kPa ∘ C −1 ) at the daily mean air temperature ( ∘ C), Rn and G are the net solar radiation and soil heat flux density (MJ m −2 d −1 ), γ is the psychrometric constant (kPa ∘ C −1 ), T is the daily mean temperature ( ∘ C), U 2 is the mean wind , e s is the saturation vapor pressure (kPa), and e a is the actual vapor pressure (kPa) (Allen et al. 1998).

Remote sensing data
MODIS is a NASA satellite sensor used to monitor Earth's environment onboard Terra and Aqua (launched in 1999 and 2002, respectively). Terra passes the equator at approximately 10:30 AM local time from north to south. Aqua crosses the equator from north to south at approximately 1:30 PM local time (Pagano and Durham 1993). Multiple bands of the MODIS sensor are capable of resolving spatial resolutions of 250 m, 500 m, and 1 km. Land surface conditions have been monitored by satellite observations for many years. As well as providing access to valuable data for regions without land observations, these observations also have the advantage of offering wide coverage. It would be possible to access a wide range of commonly Landsat and MODIS were used as remote sensing products, easily through their respective Google Earth Engine (GEE) applications (Amani et al. 2020; Moore and Hansen 2011;  Tamiminia et al. 2020). Despite the lower temporal resolution of remote sensing data, ET 0 can be estimated with the daily MODIS datasets (Kim et al. 2020). Terra and Aqua MODIS products, which were used in the current study, include LST day and night as well as NDVI, LAI, and FPAR (Table 2).

LST
An important parameter in the study of land surface processes is land surface temperature (LST) (Wan and Li 1997). NASA's Terra and Aqua Earth Observation Systems generate multiple LST products daily with MODIS sensors (Wang et al. 2008). During the period 2003 to 2021, LST was downloaded using the thermal infrared bands of the Daily Aqua (MYD11A1) and Terra (MOD11A1) satellites. LST is available as a daytime and nighttime product with a spatial resolution of 1 km using the products MYD11A1 and MOD11A1 (Zhu et al. 2013). Three types of LST were tested to assess the remote sensing-based ET 0 estimation model. The first model used a combination of Aqua and Terra products daytime LST (LST day ), the second model was based on the combination, of Aqua and Terra products nighttime LST (LST night ), and the third model utilized the average LST of day and night (LST mean ). Also, because LST values are often affected by clouds and other atmospheric disturbances, the results have more missing data (Yu et al. 2015). Due to the fact that clouds are constantly moving in both position and extent, one pixel that is cloudy in the morning (10:30 AM, Terra overpass time) may appear clear in the afternoon (1:30 PM, Aqua overpass time). The combination of Terra and Aqua products, under most meteorological conditions, would therefore provide cloudless composited images every day (Xie et al. 2009). Thus, Aqua and Terra products were combined, implementing the imputeTS library in R programming (Moritz and Bartz-Beielstein 2017), and time series data were reconstructed using Kalman filter interpolation (Andronis et al. 2022). The mentioned  There are a number of vegetation indices available, but NDVI is one of the most widely used (D'Odorico et al. 2013). In the current research, the cubic spline method was used to convert the 16-day product into daily values since vegetation greenness generally changes over a year in the form of a cosine curve (Kim et al. 2020). An alternative method of estimating leaf biomass is to use the leaf area index (LAI), which is based on a crop-specific growth coefficient and maximum gross primary production (GPP). The fraction of photosynthetically active radiation (FPAR) is the portion of photosynthetically active radiation absorbed by green vegetation. In this study, to obtain LAI and FPAR values, the MCD15A3H V6 level 4, 4-day product with a resolution of 500 m was used. Dataset MCD15A3H V6 level 4 contains pixel sizes of 500 m for a composite 4-day period. Based on all the measurements taken by NASA's Terra and Aqua satellites over a period of 4 days, the algorithm in its algorithmic form picks the best pixel available from all the acquisitions from both MODIS sensors located on those satellites (Zhen et al. 2021). The 4-day LAI and FPAR data were converted to daily data using cubic spline interpolation (Chen et al. 2006).

Multilayer perceptron (MLP)
In hydrologic modeling, artificial neural networks (ANNs) have been extensively applied to capture and demonstrate intricate input/output interactions. An ANN consists of a collection of nodes (neurons) organized into layers and linked via weights. Three layers are required for a simple neural network to function: the input layer (information), the hidden layers (secret information), and the output layer (yield). It is the input layer that presents the input information in the modeling process. The neural network then functions on this information to produce the output, which is achieved at the layer of output. Neurons with weighted networks provide many inputs to each neuron. As the weighted inputs are added up, a transmission function is constructed, for example, a linear, logistic, or hyperbolic tangent function (Talebizadeh et al. 2010). Because the ANN does not require comprehensive statistics regarding the physical process prevailing in the systems, it is an effective method for analyzing complex hydrological processes. Several studies have utilized a variety of ANNs and training algorithms (Shadkani et al. 2021). The multilayer perceptron (MLP) method is the most commonly used ANN method in hydrological research. There are two categories of ANNs, feed-forward and backward, depending upon the direction of data transmission and processing. MLPs are part of a feed-forward network. According to previous studies, MLP networks can properly estimate nonlinear functions with a high level of accuracy (Schalkoff 1997). Multilayer perceptron networks were used in the current research, and backpropagation was used for training the neural network.
The MLP method is expressed as follows: W kj represents the weights between hidden layers and output layers; W ji is used to represent weights between hidden and output layers; X i represents the input variables; m refers to the number of hidden layer neurons; the number of neurons in the input layer is n; in the hidden layer, B j represents the bias amount of the neuron; in the output layer, B k represents the bias amount of the neurons; Y represents the output function; F represents the transmission function. ANNs are often activated by using a sigmoid function. Sigmoid (Sig) functions were used in the current study. For different Z values, the Sig transmission function is as follows:

Multilayer perceptron-stochastic gradient descent (SGD-MLP)
A stochastic gradient descent algorithm (SGD) simplifies a problem by computing the gradient of the objective function based on a small subset that is chosen by chance (Bottou 2010). The number of training datasets used for a single iteration determines the batch size. In order to accelerate convergence, a small batch size can enable parameter updates more frequently than gradient descent. The maximum frequency of updates and simple perceptron-like algorithm are provided by batch sizes of 1. SGD approximates the gradient of the objective function given by Eq. 4 by randomly selecting a small subset of the training samples. The batch size is the number of training samples used. SGD training is slow to converge when the batch size is N. One can speed up convergence by updating the parameters more frequently with a small batch size than with gradient descent. When the batch size is 1, the maximum frequency of updates is achieved, and it is a straightforward perceptron-like operation.
where C is the meta-parameter that controls the degree of regularization, which is usually tuned by cross-validation or using the held-out data, N is the number of training samples, and i is the weight of the feature. As long as a single training sample is used to approximate the gradient, the optimization procedure is identical to simple gradient descent, so the weights of the features are updated as follows at training sample j: (Eq. 5) The SGD-MLP method is a neural network-based model which is optimized using SGD. Among the advanced features of SGD's algorithm are epochs, rho, L1 or L2 adjustment, momentum training, adaptive learning rate, and rate annealing that provide high accuracy in prediction when using MLP models. Aside from optimizing MLP results with SGD, the network also contains hidden layers containing neurons with tanh (hyperbolic tangent function), rectifier (where x is the input value, select the maximum of (0, x)), max out (select the maximum input vector coordinates), and ExpRectifier (exponential rectifier function). When the adaptive learning rate is inactivated, the user-determined learning rate describes the size of the weight updates based on the difference between the predicted and target values. The variance is referred to as the delta, which is only shown at the output layer. The output of each hidden layer is estimated accurately using backpropagation. Considering that redundant momentum may cause oscillations, momentum is ramped up gradually. The rate annealing, momentum training, dropout rate annealing, and momentum training parameters are activated if the adaptive learning rate is disabled.

Evaluation criteria
In the current study, four performance indexes, including the coefficient of determination (R 2 ), the root-mean-square error (RMSE), Willmott's index of agreement (WI), and Nash-Sutcliffe efficiency (NS), are utilized to measure the model's accuracy. The mathematical description of (4) R 2 , RMSE, WI, and NS can be expressed, respectively, as follows: where N is the number of observed ET data, ET i(observed) and ET i(model) are observed and estimated ET, respectively, and ET i (observed) is the mean of observed ET.

Developing algorithms for two different agro-climatic zones
Input variables from satellite images were combined to estimate ET 0 at Tabriz (semi-arid climate) and Babolsar (humid climate) stations using an optimization algorithm of SGD and MLP. Table 4 shows the input parameters for each model. There are three models: first is based on the land surface temperature, second on the vegetation indices, and third is a combination of the land surface temperature and the vegetation indices. The results of the statistical parameters for the two mentioned stations are shown in Table 5 and Fig. 3. With all considered parameters, ET 0 was estimated with the best accuracy at the two studied stations. The first scenario performs better than the second based on a comparison of the input parameters at both stations. Accordingly, the effect of LST on estimating ET 0 is higher than the effect of vegetation indices. Also, model performance can be affected by cloudy weather. Cloudy weather reduces the accuracy of satellite image parameters. However, LST data are more affected by sky cloudiness than vegetation indexes. Table 6 shows the number of missing LST data for the two mentioned stations obtained from MODIS sensors. It also shows the average of cloudy hours derived from meteorological station data for 2003-2021 for these two stations. Based on Tables 5 and 6, the lower accuracy of the statistical parameters at Babolsar station is attributed to a large number of missing LST data and more cloudy hours than at Tabriz station. In addition, the SGD-MLP-3 provided the most accurate results during the testing period, i.e., R 2 , RMSE, NS, and WI, as 0.949, 0.417 mm/day, 0.894, and 0.973 for Tabriz and 0.864, 0.754 mm/day, 0.742, and 0.922 for Babolsar station, respectively. In addition, the MLP-3 algorithm showed weaker performance than the hybrid SGD-MLP models with the performance indices, R 2 , RMSE, NS, and WI, as 0.939, 0.449 mm/ day, 0.877, and 0.969 for Tabriz and 0.847, 0.804 mm/day, 0.706, and 0.911 for Babolsar station. In order to evaluate the algorithms and models, a radar chart was used to evaluate the best-calculated values of RMSE (mm/day) for all constructed models. The obtained results for Tabriz and Babolsar stations are shown in Fig. 2. A hexagon with a smaller RMSE value indicates the best performance of the station. Based on Fig. 4, all models (SGD-MLP and MLP) at Tabriz station predicted ET 0 values more accurately than those at Babolsar station. Further, Fig. 3 shows that SGD-MLP-3 performed better at both stations than the other models in estimating ET 0 values. Moreover, Fig. 4 displays scatter plots of FAO-PM ET 0 values and estimated ET 0 by SGD-MLP and MLP models. The regression line in the scatter plots provides the coefficient of determination (R 2 ) for both models. At both stations (Tabriz and Babolsar), SGD-MLP-3 (all input parameters) produced more accurate ET 0 values. A comparison of all input parameters for two stations with the MLP algorithm demonstrated a lower level of accuracy. Further, the models   overestimated ET 0 values at Babolsar station as the regression line is above the 1:1 line. In order to evaluate the models, a radar chart was used to evaluate the best-calculated values of RMSE (mm/day) for all constructed models (Fig. 5). A hexagon with a smaller RMSE value indicates the best performance of the station. So, it can be concluded that all models (SGD-MLP and MLP) at Tabriz station estimated ET 0 values more accurately than those at Babolsar station. Moreover, Fig. 4 reveals that SGD-MLP-3 estimated ET 0 values better than other established models at both stations. By comparing root-mean-square error, standard deviation, and correlation coefficient, Taylor's diagram provides a reliable comparison of different models' performance. At Tabriz and Babolsar stations, the superiority of SGD-MLP and MLP models with different input variables was investigated using the Taylor diagram (Fig. 6). In the Taylor diagram, the red square represents FAO-PM's daily ET 0 . Generally, models near the red square are considered to be more accurate. Based on Fig. 6, MLP-3 and MLP-SGD-3 offer the weakest and most accurate estimates of ET 0 in both stations. Additionally, the accuracy of ET 0 estimation at Tabriz station is better than that at Babolsar station across all models. Figure 7 shows the estimated daily ET 0 by SGD-MLP-3 at Tabriz and Babolsar stations in the time period of 2017-2021. The estimate of ET 0 illustrated in Fig. 7 reflects the same trend as the observed ET 0 calculated using the FAO-PM and ground meteorological observations at two locations.
In Table 7, the cross-correlation between the input parameters and the target is shown. There is a strong relationship between the LST mean in Tabriz station and the FAO-PM ET 0 based on all of the studied inputs (R = 0.792). Moreover, the NDVI at Babolsar station has a weak relationship with the   (R = 0.196). According to the results of current research, the land surface temperature (day, night, and average) in both stations have a stronger relationship with ET 0 than vegetation indices. Also, due to the high number of cloudy days at the Babolsar station and the low quality of the satellite images, the correlation of all parameters is lower than correspondent values at Tabriz station. The results of previous studies were also compared with the findings of the current research. Using artificial neural networks (ANNs) and remote sensing data, Chen et al. (2013) developed a model to estimate daily evapotranspiration (ET). Surface net radiation, normalized difference vegetation index, and land surface temperature were derived from satellite remotely sensed data. According to their results, there is a coefficient of determination of 0.770 between estimated and measured ET values. Moreover, Zhang et al (2018) compared the three algorithms of support vector machine (SVM), backpropagation neural network (BP), and adaptive neuro-fuzzy inference system (ANFIS) with different combinations of MODIS sensor input variables for estimation of daily ET 0 in Northwest China. Among the three models, the SVM model with a correlation coefficient of 0.917 had the best results. The obtained results in the current study are consistent with the mentioned findings.

Conclusion
For agricultural activities, estimating ET 0 is crucial. Through satellite data with a high temporal resolution, new insights can now be gained into agriculture. Consequently, satellite images provide an opportunity to determine the Fig. 6 Taylor diagrams for ET0 estimates for both models based on three scenarios at two stations Fig. 7 Evaluation of best algorithms SGD-MLP-3 (for Model-3) at Tabriz (semi-arid climate) and Babolsar (humid climate) station most critical ET 0 -related parameters, including LST during the day, NDVI, LAI, and FPAR from the MODIS sensor. The present study used satellite images to estimate ET 0 in croplands of Tabriz (semi-arid climate) and Babolsar (humid climate) stations, Iran, using MLP and SGD-MLP models. A training dataset from 2003 to 2016 was used for training, and a testing dataset from 2017 to 2021 was used for testing. According to the obtained results, the SGD-MLP-3 gave the most accurate results with values of performance R 2 , RMSE, NS, and WI as 0.949, 0.417 mm/day, 0.894, and 0.973 for Tabriz and 0.864, 0.754 mm/day, 0.742, and 0.922 for Babolsar station, respectively. So, there was a constant trend that revealed greater accuracy of both SGD-MLP and MLP models at Tabriz station than at Babolsar station. Moreover, at both studied stations, ET 0 values were estimated more accurately by model SGD-MLP-3. Also, at the Babolsar station, the correlation between satellite parameters with ET 0 values was lower than Tabriz station due to the high number of cloudy days and the low quality of the satellite images. In conclusion, the implementation of hybrid SGD-MLP models utilizing satellite image parameters was recommended for ET 0 estimation, especially in semi-arid regions.
Funding The author's received no specific funding for this work.

Data availability
The datasets analyzed during the current study are available online at https:// figsh are. com/ artic les/ datas et/ babol sar_ xlsx/ 21674 033

Conflict of interest The authors declare no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.