Hybrid Neural Network Based Models for Evapotranspiration Prediction Over Limited Weather Parameters

Evapotranspiration can be used to estimate the amount of water required by agriculture projects and green spaces, playing a key role in water management policies that combat the hydrological drought, which assumes a structural character in many countries. In this context, this work presents a study on reference evapotranspiration (<inline-formula> <tex-math notation="LaTeX">$ET_{o}$ </tex-math></inline-formula>) estimation models, having as input a limited set of meteorological parameters, namely: temperature, humidity, and wind. Since solar radiation (SR) is an important parameter in the determination of <inline-formula> <tex-math notation="LaTeX">$ET_{o}$ </tex-math></inline-formula>, SR estimation models are also developed. These <inline-formula> <tex-math notation="LaTeX">$ET_{o}$ </tex-math></inline-formula> and SR estimation models compare the use of Artificial Neural Networks (ANN), Long Short Term Memory (LSTM), Gated Recurrent Unit (GRU), Recurrent Neural Network (RNN), and hybrid neural network models such as LSTM-ANN, RNN-ANN, and GRU-ANN. Two main approaches were taken for <inline-formula> <tex-math notation="LaTeX">$ET_{o}$ </tex-math></inline-formula> estimation: (i) directly use those algorithms to estimate <inline-formula> <tex-math notation="LaTeX">$ET_{o}$ </tex-math></inline-formula>, and (ii) estimate solar radiation first and then use that estimation together with other meteorological parameters in a method that predicts <inline-formula> <tex-math notation="LaTeX">$ET_{o}$ </tex-math></inline-formula>. For the latter case, two variants were implemented: the use of the estimated solar radiation as (ii.1) a feature of the neural network regressors, and (ii.2) the use of the Penman-Monteith method (a.k.a. FAO-56PM method, adopted by the United Nations Food and Agriculture Organization) to compute <inline-formula> <tex-math notation="LaTeX">$ET_{o}$ </tex-math></inline-formula>, which has solar radiation as one of the input parameters. Using experimental data collected from a weather station (WS) located in Vale do Lobo (Portugal), the later approach achieved the best result with a coefficient of determination <inline-formula> <tex-math notation="LaTeX">$(R^{2})$ </tex-math></inline-formula> of 0.977. The developed model was then applied to data from eleven stations located in Colorado (USA), with very distinct climatic conditions, showing similar results to the ones for which the models were initially designed (<inline-formula> <tex-math notation="LaTeX">$R^{2}>0.95$ </tex-math></inline-formula>), proving a good generalization. As a final notice, the reduced-set features were carefully selected so that they are compatible with free online weather forecast services.

theft, since normally there is no comprehensive security or surveillance in such locations. Further, the available water holding capacity changes significantly with soil type [2], requiring that for each specific soil, in order to determine the soil humidity values corresponding to the field capacity and wilting point (WP), samples would need to be sent to a laboratory for analyses.
Crop evapotranspiration (ET c ), also known as crop water use, is the water that is used by a crop [3]. The Food and Agriculture Organization of the United Nations (FAO) recommends using the FAO-56 Penman-Monteith (FAO-56PM) formula as a reference method for computing reference evapotranspiration (ET o ) [4], being ET c and ET o related by a crop coefficient (K c ). FAO56-PM formula uses four main meteorological parameters: temperature, humidity, wind, and solar radiation (SR). In this context, several computational studies show that SR is the main driver of ET o [5], however, its measurement requires sensors like pyranometers, which are typically associated with expensive WSs, that need to be properly maintained and calibrated [6]. Also, solar radiation forecast application programming interfaces (APIs) are not common (at least freely) and present a high system cost penalty.
As part of a framework for the computation of optimal crop water irrigation scheduling requirements (with special emphasis to green spaces), this paper presents the computational models being prepared to estimate the ET o using machine learning, deep learning, acquired intelligence, meteorological data from WSs on the field, as well as meteorological data and forecasts from APIs available on the internet. This will optimize water and energy expenditure, improve the well-being of the crop, reduce reaction time in solving problems, improve anomalies detection methods, and maintain the quality of green spaces. Concisely, the framework being developed will be an intelligent irrigation solution, technologically differentiated from other platforms on the market. The development of the full framework is being done under project GSSIC -Green Spaces SMART Irrigation Control.
In short, this study compares several ET o and SR estimation models that use Artificial Neural Networks (ANN), Long Short Term Memory (LSTM), Gated Recurrent Unit (GRU), Recurrent Neural Network (RNN), and hybrid neural network models such as LSTM-ANN, RNN-ANN and GRU-ANN. Two approaches were taken for ET o 's estimation: (i) directly use those machine learning models to estimate it and (ii) first estimate solar radiation, and then use the obtained value and other meteorological parameters in a method that predicts ET o . Furthermore, for the latter case, two variants were addressed, namely: the use of the estimated SR as (ii.1) a feature of a neural network regressor, and (ii.2) the use of FAO-56PM method to compute ET o , which has SR as one of the input parameters. Figure 1 schematizes what we have just detailed.
Using experimental data collected from a WS located in Vale do Lobo, south Portugal, approach (ii.2) achieved the best result with a coefficient of determination (R 2 ) of 0.977, improving results previously achieved in [7]. The quality and generalization of the proposed models was then tested using 11 weather stations of the Colorado Agricultural Meteorological Network (CoAgMET), which is composed of automatic weather stations distributed across the state of Colorado, USA. The selected weather stations have data available starting from August 1992, which comprises almost 30 years of data. Also, with temperatures ranging from −39.06 to 44.26 degrees Celsius ( • C), they represent different climatic conditions than the Mediterranean's Vale do Lobo WS. This allows the assessment of model performance under different climate conditions, and for a longer period of historical data. The achieved metrics values, when comparable with the ones from other authors, showed that our proposal has in general a better performance in the analysed metrics, besides presenting alternatives to most of the existing models, namely through the use of recurrent neural networks and hybrid methods. The latter model is even more unusual in this category of studies, being by itself another evolution to the state of the art.
The paper is structured as follows. The next section presents the problem's background and the methodologies used by others to tackle the problem in study. Section III starts by detailing the computational setup and exploring the dataset, followed by section III-B where the proposed neural network architectures, hyperparameters, and overall training approach is presented. Then, on section III-C the proposed models and associated performance analysis is weighed. Section IV presents the results for the 11 CoAgMET agricultural weather stations. Finally, the last section presents the conclusion and establishes some future work.

II. REFERENCE EVAPOTRANSPIRATION PROBLEMATIC AND KNOWN SOLUTIONS
Reference evapotranspiration, ET o , is the evapotranspiration of a reference surface, defined as hypothetical grass with a uniform height of 0.12 m, a fixed surface resistance of 70 sm −1 , and an albedo (reflection coefficient) of 0.23 [4]. Besides, crop evapotranspiration, ET c , represents the crop's water requirements and is proportional to reference evapotranspiration by means of the crop coefficient, K c [4]. Therefore, ET o prediction plays an important role, making it one of the fundamental parameters for smart irrigation scheduling, since it is proportional to the amount of water that needs to be restored during the irrigation period [8].
Some of the main characteristics that distinguish ET c from ET o are (i) the crop cover density and total leaf area, (ii) the resistance of foliage epidermis and soil surface to the flow of water vapor, (iii) the aerodynamic roughness of the crop canopy, and (iv) the reflectance of the crop and soil surface to short wave radiation [9]. In this context, known the value of K c , the ET c value for a specific time period can be estimated by (1) The crop coefficient can be simple or have two components, one representing the basal crop coefficient (K cb ) and another representing the soil surface evaporation component (K e ), being computed by where K s ∈ [0, 1] is used to introduce a K c reduction in cases of environmental stresses, such as lack of soil water or soil salinity [9]. So, it becomes clear that in order to make crop water requirement predictions, accurate estimation of ET o is required. Historically, several deterministic methods have been developed to estimate reference evapotranspiration using single or limited weather parameters and being generally categorized as: temperature, radiation, or combination based. For example, temperature based methods include Thorntwait [10], Blaney-Criddle [11], and Hargreaves and Samani [12] formulas; radiation methods include Priestley-Taylor [13] and Makkink [14] formulas; and combination methods, requiring both temperature and radiation, include Penman [15], modified Penman [16], and FAO-56 Penman-Monteith (FAO56-PM) [4] formulas. Shahidian et al. [6] give an overview of several of these methods and compare their performance under different climate conditions. The authors concluded that, when applied to climates different from those on which the methods were developed and tested, most of them yield a poor performance and may require the adjustment of empirical coefficients to accommodate to local climate conditions, which is not ideal.
FAO recommends the use of the FAO-56PM formula as a reference method for estimating ET o [4]. To give a deeper idea of the involved parameters, measured in millimeters per day [mm/day], the formula to estimate ET o is given by . Being based on physical principles, the formula has become widely adopted as a standard for ET o computation since it performs well under different climate types [6]. However, to compute ET o using FAO56-PM the following main meteorological parameters are required: temperature, solar radiation, relative humidity, and wind speed. The remaining meteorological parameters are constants being derived from latitude and elevation above sea level.
Except for solar radiation, all parameters (real or estimated) required by the FAO56-PM formula for ET o 's computation, can be freely obtained from common weather forecast APIs. Solar radiation forecasting APIs are, currently, not common and the ones available present a high-cost penalty. Therefore, a major asset would be to (i) develop alternative methods for ET o estimation using limited meteorological parameters, that do not require solar radiation and are compatible with the weather parameters obtained by freely available weather forecast and historical weather data APIs, and/or (ii) to estimate the solar radiation itself and use it as an approximation on the solar radiation dependent methods. This is also important since in most situations a proper functioning, maintained, and calibrated WS, with solar radiation measurement capability, is not close to the area of interest.
Recently, as an alternative, several authors have used machine and deep learning to estimate ET o , being that, as stated by Chia et al. [17], machine learning has proved itself to be a promising solution for ET o estimation with the common meteorological data. For instance, using data collected in Central Florida, a humid subtropical climate, Granata [18] compared three different evapotranspiration models which differ in the input variables. In their work, four variants of machine learning algorithms were applied to each model, namely: M5P Regression Tree, Bagging, Random Forest, and Support Vector Machine (SVM). Their best results are achieved using M5P Regression Tree and Bagging with a coefficient of determination of 0.987 and a mean absolute error of 0.14 mm/day (see the metrics definition in Sec. III-A). However, among other features, all models included as input variable the net solar radiation. The same author also studied three recurrent neural networkbased models for the prediction of short term ahead evapotranspiration [19]. Two variants of each model were developed, selecting the employed algorithm between LSTM and nonlinear autoregressive network with exogenous inputs (NARX). The prediction models were trained and tested using data from two sites with different climates: Cypress Swamp, southern Florida, and Kobeh Valley, central Nevada. With reference to the subtropical climatic conditions of South Florida, LSTM models proved to be more accurate than NARX models, while some exogenous variables such as sensible heat flux and relative humidity did not affect the results significantly.
Wu and Fan [20] evaluated eight machine learning algorithms divided in four classes: neuron based (MLP -Multilayer Perceptron, GRNN -General Regression Neural Network, and ANFIS -Adaptive Network-based Fuzzy Inference System), kernel-based (SVM, KNEA -Kernelbased Non Linear Extension of Arps decline model), treebased (M5Tree -M5 model tree, Extreme Gradient Boosting -XGBoost), and curve based (MARS -Multivariate Adaptive Regression Spline). The methods were applied to data collected from 14 WSs in various climatic regions of China and used only temperature or temperature and precipitation as input to the models. Ferreira et al. [8] used six alternative empirical reduced-set equations, such as Hargreaves and Samani [12], and compared the estimated values with the ones from an ANN and a SVM model. Data was collected from 203 WSs and used for daily ET o estimation for the entirety of Brazil. Temperature or temperature and humidity were used as input features. They concluded that, in general, ANN was the best performing model when including, as input features, data from up to four previous days. With the best algorithms reaching an R 2 median value around 0.80 considering all stations, results were weighed good given that only temperature or temperature and humidity were used as input.
Keshtegar et al. [21] used high-order response surface method to compute ET o . They included daily weather information which include the maximum temperature, maximum humidity, wind speed, solar radiation, and vapor pressure deficit, obtained from three observation stations in Burkina Faso, West Africa. Ten models were evaluated with the determination coefficient (R 2 ) and root-mean-square error values ranging between 0.2068 and 0.9966 and between 0.7237 and 0.9948, respectively. However, one of the used features is solar radiation, which is not always easy to obtain for most of the locations.
Muhammad et al. [22] studied the implementation of evolutionary computing models, namely the gene expression programming (GEP), for the simulation of daily ET o in different locations of Peninsular Malaysia. The models used various input combinations of meteorological variables including air temperature (mean, maximum, and minimum), relative humidity, solar radiation, and mean wind speed. Compared to other black box artificial intelligence algorithms, the authors' alleged major advantage of GEP is that it provides a set of equations that can be used by practitioners for reliable estimation of ET o at the field, with fewer meteorological variables and, thus, can have wide applicability in water resources management. As in the previous case, the model also requires the use of solar radiation with the already mentioned drawback.
Khosravi et al. [23] studied nine models, including five data mining algorithms and four adaptive neuro-fuzzy inference systems, for their ability to predict ET o at meteorological stations in Baghdad and Mosul, Iraq. As parameters for the models, they considered wind speed, sunshine hours, rainfall, maximum and minimum temperature, and relative humidity. Investigations on the modeling accuracy with different input parameter combinations showed that no single input combination showed a consistent modeling outcome. Moreover, hybrid models showed a higher predictive power than the individual models. A large part of the result shown had R 2 value around 0.9, with the best one achieving 0.951 for one station and 0.97 for the other.
The geographical robustness of various inter-model ensembles in estimating daily ET o was assessed in [17]. The study aimed to develop inter-model ensembles that consist of ANN, support vector regressors, ANFIS via Bayesian modeling approach, and the non-linear neural ensemble, trained for different meteorological stations with differential geographical characteristics in the Peninsular Malaysia. This work allowed to infer some aspects of the effect of the geographical characteristics on the performance of inter-model ensembles and examine the effect of data management strategies applied to solve the data-hungry issue (both qualitative and quantitative) of inter-model ensembles. Sanikhani et al. [24] explore 6 artificial intelligence models for modeling ET o using minimum and maximum temperatures of the air and extraterrestrial radiation. The models include MLP, GRNN, radial basis neural networks (RBNN), integrated adaptive neuro-fuzzy inference systems with grid partitioning and subtractive clustering (ANFIS-GP and ANFIS-SC), and GEP. The Hargreaves-Samani equation and its calibrated version were used to perform a verification analysis of the established models. In the best cases, results show R 2 values around 0.95.
In [25], Chen et al. study the estimation of daily reference evapotranspiration using three deep learning models (namely, deep neural network -DNN, temporal convolution neural network -TCN, and LSTM neural network). The performance of the three models was compared with the results of an SVM algorithm, a random forest algorithm, and empirical equations like Hargreaves and modified Hargreaves(temperature-based), Ritchie, Priestley-Talor, and Makkink (radiation-based), and Romanenko and Schendel (humidity-based) empirical models. Several other works can be found in the literature [26], [27], [28], [29], [30], [31], [32], [33], which vary the number and type of features, estimation models etc.
Vaz et al. [7] used data from a Vale do Lobo WS, in south Portugal, and explored the use of machine learning for ET o estimation. They concluded that instead of directly estimating ET o , the best result was obtained by using machine learning for SR estimation, having as input a limited set of meteorological features, and then use that result together with temperature, humidity and wind speed as input to the FAO56-PM equation, achieving an R 2 of 0.975, a MAE of 0.18 mm/day, and a MAPE of 5.51 %. The work here presented explores and develops deep learning based ET o prediction models, supported on the same data, improving the previous results, as it will be detailed in the next sections. Furthermore, the proposed models are applied to 11 other WS stations maintaining, without further hyperparameter tuning, similar performance in terms of attained metrics values. So, based on a set of features generally available in public WS, which do not include solar radiation, our solution produces a stable set of results for different weather environments presenting a state-of-the-art methodology.
Data from Vale do Lobo WS, in south Portugal, was collected starting from February 2019 up to and including September 2021. The WS is composed of sensors from Davis Instruments, where the following weather parameters are measured periodically throughout the day and stored with a daily resolution: temperature (minimum, maximum, and average), dew point (minimum, maximum, and average), relative humidity (minimum, maximum, and average), solar radiation (maximum and average), wind speed (minimum, maximum, and average), wind direction, atmospheric pressure (minimum, maximum, and average), rain intensity, and precipitation. This is the same dataset previously used by Vaz et al. [7], where the use of machine learning algorithms was explored to create ET o and SR estimation models. A ratio of 75 % to 25 % of train and test data was used, respectively, resulting in train data starting from February 1st, 2019 up to February 3rd, 2021, and test data from February 4th, 2021 up to September 30th, 2021. Being a time series, no shuffling was made to the train and test data, and the train data was further divided into 10 folders, used to implement time series cross validation [41]. Furthermore, a hyperparameter grid search strategy was used to tune the proposed machine/deep learning methods, as will be presented in the corresponding sections.
For model statistical evaluation and performance comparison, the mean absolute error (MAE), the mean absolute percentage error (MAPE), the mean square error (MSE), the root-mean-square error (RMSE), and the coefficient of determination (R 2 ) were used [42]. Just to recall, considering y t the actual value andŷ t the estimated value at instants t = 1, 2, . . . , n, and y the mean value of the actual samples, the evaluation metrics are defined as and In this set of metrics, the MAE measures the average of the absolute residuals in the original unit and the MAPE measures the average absolute percent error. The MSE represents the average of the squared difference between the original and predicted values in the data set, measuring the variance of the residuals in the squared unit of the original data, which can be returned to the original unit by computing the RMSE (the standard deviation of residuals). Finally, R 2 represents the proportion of the variance in the dependent variable, which is explained by the linear regression model, being a scale-free score. While values should be near to zero for the MAE, MAPE, MSE, and RMSE, they should be near to 1 for the R 2 metric. Depending on the regression problem in study (ET o vs. SR), two targets were considered: (i) ET o was computed using the FAO56-PM formula, as per Eq. (3), and using as input the data measured by the already referred WS, namely temperature, humidity, wind speed and solar radiation (see Sec. III-C1, III-C3, and III-C4); and (ii) the average solar radiation measured by the WS was used as SR target (see Sec. III-C2). VOLUME 11, 2023

B. NEURAL NETWORK MODELS' ARCHITECTURE
The following neural network types were used for the conduction of this work: ANN, LSTM [43], GRU [44], and Artificial RNN [45], as well as hybrid models like LSTM-ANN, RNN-ANN, and GRU-ANN.
Artificial neural networks typically have an input layer, one or more hidden layers, and one output layer. The input layer dimension is defined by the number of inputs to the model, i.e., the number of features that are used on a particular model. Hidden layers are internal layers that can vary in quantity and number of neurons. The output layer can have one or more neurons, depending on the number of outputs a model has. Artificial neurons are activated using activation functions, which have an important role on the performance of neural networks, since they provide the non-linearity that is needed to learn a complex problem [46]. During the conduction of this work, Rectified Linear Unit (ReLU), Tangent Hyperbolic Function (Tanh) and Sigmoid activation functions were tested, however ReLU always yielded better results, being therefore fixed [47].
LSTM, GRU, and RNN algorithms are deep learning methods that, unlike feed-forward neural networks, implement feedback connections. These feedback connections provide the ability to add memory to the models, making its use justified through the fact that the dataset (and target) is a time series, where some patterns might be cyclic or information from the previous days might play an important role [43], [44], [45].
Furthermore, in the work presented by Vaz et al. [7], for the same problem and dataset here in study (using limited weather parameters as input features), several machine learning regression models were compared in their performance, namely: Ordinary Least Squares, Ridge, Lasso, k-Nearest Neighbors, Support Vector Machine, Decision Tree, and Random Forest (RF). Since Random Forest [48], [49] gave the best results for ET o and SR estimation models, RF will also be evaluated, so that the proposed neural network based models can be directly compared with other top performing methods.
Systematic testing (varying the number of neurons in powers of 2) allowed us to set the non-hybrid models to have two hidden layers: the first one consisting of 512 neurons, followed by a second layer that has 32 neurons, represented as having [512,32] neurons. This configuration consistently gave the best results, being observed that augmenting the number of layers and/or neurons of the neural network would not increase model performance, while reducing would impact on model's performance. Following the same arrangement, allowed us to decide the hybrid models' architecture (namely, LSTM-ANN, RNN-ANN, and GRU-ANN). In this case, to balance the training computational requirements and the models' performance, it was decided to use four hidden layers. E.g., as depicted in Fig. 2, the LSTM-ANN model has the first two hidden layers consisting of [32,64] neurons of LSTM type, followed by two more layers  [32, 64; 64, 32] neurons. Also, as can be seen in Fig. 2, between each layer there is the possibility of using dropout regularization.
Solver, loss function, and hyperparameters like learning rate, number of epochs, and batch size play an important role on model training and performance [50]. For all models, the Adam optimizer [51] was selected as a solver, mean squared error (MSE) was used as the loss function, kernel initialization was done using the Glorot normal  initialization [52], and the number of epochs was set to a high value of 1000, however, early stopping was used, with a patience value set to 150. Dropout regularization [53] was applied during training and its value was selected using grid search with values ranging from 0.0 up to 0.6, in 0.1 steps. During initial model development and training, the batch size was also selected using a grid search approach, with a range from 64 up to 256, in multiples of 32. However, a batch size of 128 always gave a good and consistent result, so later it was fixed to 128 for all models, reducing models training time. Input features were normalized to values between 0 and 1. Finally, max norm weight constraint was applied, and its value was grid searched between 1 and 4, in steps of one.
It should also be noticed that time lag was introduced on all models, i.e., data from the previous days was introduced as an input feature together with current day data. In this case, time lags from 1 up to 10 days were tested, however it brought no improvement on models performance. In this experimental context, attempts of introducing new features through the use of polynomial features and features inverse were also made, but did not present any further improvement in models' metrics and, as such, it is not presented here.
In appendix, Tables 5 and 6 summarize the tested and tuned hyperparameters, respectively, for each of the proposed models.

C. ET o ESTIMATION USING ANN METHODS WITH A LIMITED SET OF FEATURES
Solar radiation is not normally available either as a measurement or as a forecast, therefore the need to develop models to estimate ET o that do not require it as an input feature. In this context, ET o estimation models that have as input a limited set of parameters are explored in Sec. III-C1. Furthermore, since solar radiation is the main factor for the determination of ET o [4], in Sec. III-C2 solar radiation estimation models that use a limited set of features are explored. Finally, two approaches are taken: (i) use the previously estimated solar radiation as an input to another neural network model (Sec. III-C3) or (ii) use the FAO56-PM formula to compute ET o , using as an input feature the estimated solar radiation (Sec. III-C4). Please refer to Fig. 1 for scheme of the proposed flow.

1) ET o ESTIMATION USING ANN METHODS (EXCLUDING SOLAR RADIATION)
In this section neural network based models are used to directly estimate ET o , using as inputs limited weather parameters. The set of features used are Month ∈ {1, 2, . . . , 12}, maximum and minimum temperature (TempMax and TempMin), average humidity (HumididtyAvg), and average wind speed (WindAvg). Table 1 shows the results obtained by the non-hybrid models and RF (included for comparison purposes), being    Fig. 3 shows that the estimator follows relatively well the ET o target value. In the same figure, the shadowed region corresponds to the test data, being visible the expectable slight increase of the absolute error, when comparing to the train region (as stated earlier, all presented metrics are calculated using only the test data). Table 3 also presents the results obtained using the hybrid models, namely, LSTM-ANN, RNN-ANN, GRU-ANN. The results are slightly worse than the ones obtained using the GRU, with the LSTM-ANN model attaining an R 2 of 0.959 against the 0.962 of GRU (MAE, MAPE, MSE, and RMSE are also similar).

2) SOLAR RADIATION ESTIMATION USING ANN METHODS WITH A LIMITED SET OF FEATURES
As an alternative to the presented in the previous section, and since solar radiation is the main driver of evapotranspiration [4], [5], this section studies solar radiation estimation models that have as input a limited feature set, compatible with the limited number of features returned by the common weather forecast APIs. The estimated result will be later injected as a feature in other machine learning model or used in the FAO56-PM formula (Eq. 3).
So, in this section, the solar radiation measured by the WS was used as target. Initially, the features used in the solar radiation models were the same as the ones proposed by Vaz et al. [7], namely: Month, Day, maximum and minimum temperature (TempMax and TempMin), average humidity (HumididtyAvg), average wind speed (WindAvg), dew point, and the polynomial feature Month 2 × Day. However, it was found that the model metrics were improved by dropping dew 970 VOLUME 11, 2023 point and the polynomial feature Month 2 × Day and, instead, adding the sunset hour angle (ω s ) and daylight hours (N ). As a note, the sunset hour angle (ω s ) was defined as [4] where ϕ is the latitude of the WS (in radians) and is the solar declination (also in radians), and J is the number of the day in the year. Furthermore, the daylight hours (N) is given by The results obtained are summarized in Table 2, where it can be seen that RNN and ANN are the best performing methods, with RNN being slightly better with an R 2 of 0.833, a MAE of 18.46 W /m 2 /day, a MAPE of 10.30 %, and a RMSE of 28.817. This result is better than what was obtained by Vaz et al. [7], using the same dataset, on which Random Forest gave the best results, with an R 2 of 0.814, a MAE of 21.31 W /m 2 /day, and a MAPE of 11.29 %. The inclusion of daylight hours and sunset angle improved performance for all neural network based models but, as can be seen on Table 2, it worsened Random Forest performance. Fig. 4 depicts the target solar radiation that was measured by the solar station (blue), the approximated solar radiation obtained with ANN method (orange), and the absolute error curve (green). Shadowed is the test set. This solar radiation estimation will be used next to predict the ET o values.
As in Sec. III-C1, the recursive hybrid LSTM-ANN, RNN-ANN and GRU-ANN models' results, summarized also in Table 2, give similar performance to the recursive non-hybrid methods, with the advantage of requiring less VOLUME 11, 2023 TABLE 4. Comparison of several regression methods for ET o estimation using a limited set of features. Data with a span of almost 30 years was collected from 11 WS spread across the Colorado state, and model performance metrics from the 11 WS was averaged for each of the developed models. computation power, since the number of trainable parameters is highly reduced due to the smaller network that is used. E.g., the LSTM model has 1,128,609 trainable parameters, while the LSTM-ANN model only has 35,841, requiring approximately 21 % of the training time, and 58 % of the inference time.

3) ET o ESTIMATION USING THE APPROXIMATED SOLAR RADIATION
This section describes the implementation of the ET o estimation using as features solar radiation (approximated using the ANN presented in Sec. III-C2) together with maximum temperature, average humidity, and average wind speed. To be more clear, the difference between the present models and the ones in Sec. III-C1 is the use of the (estimated) SR, which was not previously used.
The obtained results are summarized in Table 3. The ANN model was the best one with an R 2 of 0.968, a MAE of 0.22 mm/day, a MAPE of 6.66 %, and a RMSE of 0.303 (against the R 2 of 0.962, the MAE of 0.24 mm/day, the MAPE of 7.25%, and the RMSE of 0.331 previously obtained with GRU model). Similarly, with the advantages previously stated, the LSTM-ANN model achieved an R 2 of 0.967, just one thousandth worse than the ANN model, but improving the MAPE to 6.5%. Also, the obtained result is better than the ones presented in Vaz et al. [7], using a similar technique but based on ML algorithms, where the best performing model was Random Forest with an R 2 of 0.951, a MAE of 0.26 mm/day, and a MAPE of 7.44 %. Feature engineering as well as the use of other weather limited features was attempted, but no further improvements could be made. Fig. 5 (top) depicts the target ET o (blue), estimated ET o (orange), and absolute error (green) curves for the train and test (shadowed) dataset, allowing to observe the model closely follows the reference evapotranspiration target.

4) ET o ESTIMATION USING FAO56-PM EQUATION AND THE APPROXIMATED SOLAR RADIATION
In this section, another type of hybrid approach was tested, which uses the previously estimated solar radiation as an input, together with temperature, humidity and wind speed to compute ET o using the FAO56-PM formula (SR-ANN → FAO56-PM). The result obtained is presented in Table 3, with an R 2 of 0.977, MAE of 0.16 mm/day and MAPE of 5.05 %. This result is better than the previously obtained ones, and once again is also better than what was presented by Vaz et al. [7] which consisted of an R 2 of 0.975, an MAE of 0.18 mm/day, and an MAPE of 5.51 %. Fig. 5 (bottom) plots ET o target, estimated ET o , and absolute error, where in the bottom plot an improvement of the absolute error can be seen when compared with the previous plots.
Next section will generalize the current results to 11 WS ran by the Colorado Agricultural Meteorological Network (CoAgMET).

IV. CoAgMET AGRICULTURAL WEATHER STATIONS
The Colorado Agricultural Meteorological Network (CoAg-MET) is a network of automatic weather stations distributed across the state of Colorado. The collected weather data, such as temperature, humidity, wind speed and solar radiation, can be accessed via an API at www.coagmet.com. Data from 11 weather stations was collected, namely the stations having the following code names: alt01, avn01, ctz01, hyk02, idl01, ksy01, lcn01, oth01, pkh01, rfd01, yjk01. These weather stations were chosen because they are spread across the Colorado state and have data available starting from August 1992, up to and including July 2021 (when the test were run), which comprises almost 30 years of data. Also, with temperatures ranging from −39.06 to 44.26 degrees Celsius ( • C), and having desert zones, the Colorado state poses different climatic conditions than the Mediterranean based Vale do Lobo. This allows the assessment of model performance under different climate conditions than Vale do Lobo, and for a longer period of historical data. For each WS, the ANN, the LSTM-ANN and the RF models were trained using the same hyperparameters found for the Vale do Lobo Dataset (see Table 6). As previously, a ratio of 75 % to 25 % of train and test data was used, meaning that the developed models are making predictions for more than 7 years. The results obtained for the 11 WS were averaged and are presented in Table 4. In short, the best results were attained by the hybrid method where SR is estimated by an ANN and then injected in the FAO56-PM formula (SR-ANN → FAO56-PM), achieving an R 2 of 0.984, a RMSE of 0.261 mm/day, a MSE of 0.069 mm 2 /day, a MAE of 0.16 mm/day, and a MAPE of 4.84 %.
Furthermore, R 2 and RMSE variation across all stations is small, as can be observed in the boxplots of Fig. 6. Also, for all models, R 2 is always above 0.958 and RMSE is below 0.45 mm/day, which can be considered a good result. MAE metric variation across all WS is also small, as it can be observed in Table 4 and inferred from the boxplots on Fig. 7 which shows that, for the 7 years of test data, the maximum absolute error is below 1 mm/day, apart outliers. As a note, values are regarded as outliers if they are more than 1.5 times the interquartile range (IQR) below the first quartile (Q 1 ) or more than 1.5 × IQR above the third quartile (Q 3 ), i.e., values outside the interval [Q 1 − 1.5 × IQR, Q 3 + 1.5 × IQR].

V. CONCLUSION AND FUTURE WORK
Evapotranspiration can be used to estimate the amount of water required by agriculture projects and green spaces, playing a key role in water management policies that combat the hydrological drought. Some equipment can be used to measure ET o , but they are expensive to buy and maintain. As an alternative, many use meteorological data to estimate the ET o values, but again, in the most well accepted formula (FAO56-PM), solar radiation is required, which is not common to have in general meteorological APIs.
In this work, several neural network based regression models (namely, ANN, LST, GRU, and RNN) and hybrid approaches (namely, LSTM-ANN, RNN-ANN, GRU-ANN) for ET o estimation were developed with different degrees of success. Since solar radiation is the main ET o driver, as stated by several authors, models were also developed for estimating solar radiation using input features that are readily available in common weather forecast APIs. This allowed both the use of the previously estimated solar radiation in neural network regressors, to estimate ET o , but also the possibility to use the hybrid approach where solar radiation is previously estimated, and then the FAO56-PM method is used to finally compute ET o . For the Vale do Lobo WS, the latter yielded the best results, with an R 2 of 0.977, a RMSE of 0.256 mm/day, a MSE of 0.066 mm 2 /day, a MAE of 0.16 mm/day, and a MAPE of 5.05 %. The attained results were also compared against the ones in a previous work from the authors, were several other machine learning methods were used (namely, Ordinary Least Squares, Ridge, Lasso, k-Nearest Neighbors, Support Vector Machine, Decision Tree, and Random Forest), demonstrating the overall good result, considering the limited weather parameter features that were used. Being the  Also, our work shows that the recurrent hybrid network models (e.g., LSTM-ANN) give similar results to the nonhybrid recurrent network models (e.g., LSTM). However, its computational cost is lower due to the decrease of trainable parameters, requiring only 21 % of training, and 58 % of inference time when compared to the non-hybrid recurrent networks. This is an important result when using edge computing, where computational power is limited.
Future work will include a dataset collected from the existing WS infrastructure that is installed in the Algarve region, in south Portugal. The objective will be to develop local and pooled models of ET o predictors for the Algarve region. Also, since all limited feature models here presented are compatible with the freely available weather forecast APIs, the impact of using such APIs as input data to the ML models here developed needs to be further assessed.

APPENDIX DEEP LEARNING ALGORITHM PARAMETERS
In this section, the hyperparameters and value ranges that were used in the grid search procedure are presented in Table 5. Tuned hyperparameters, for models presented in Table 1, 2, and 3, are presented in Table 6. Please note that, fixed hyperparameter values (e.g., number of neurons, and neural network architecture) are described in Sec. III-B. From 2018 to 2021, he was a High School Teacher, lecturing classes in computer science, electricity, and electronics fields. Since 2021, he has been working as a Researcher under Project GSSIC for which he is a scholarship holder, and also as a Guest Assistant Lecturer, lecturing computer programming classes with the University of Algarve. His research interests include machine learning, deep learning, data science, electronics, audio processing, and augmented creativity.

VI. ACKNOWLEDGMENT
GABRIELA SCHÜTZ received the degree in mathematics and the M.Sc. degree in statistics and operations research from the Faculty of Science, University of Lisbon, and the Ph.D. degree in mathematics-operations research from the University of Algarve. She is currently a Full (Coordenador) Professor of mathematics with the Institute of Engineering, University of Algarve. Her research interests include combinatorial optimization problems, network design problems, machine learning, and evolutionary algorithms. She is a member of the Center of Electronics Optoelectronics and Telecommunications (CEOT) and the Portuguese Mathematical Society (SPM).
CARLOS GUERRERO received the M.Sc. degree in soil fertility and plant nutrition from the Agronomic Superior Institute, University of Lisbon, and the Ph.D. degree in agriculture sciencesagro-environment from the University of Algarve. He is currently an Assistant Professor of soil sciences and remote sensing with the Faculty of Sciences and Technology, University of Algarve. He is an Agronomist with the Agronomic Superior Institute, University of Lisbon. His research interests include turfgrass maintenance, turfgrass diseases and irrigation, remote sensing, and precision agriculture. He is an Integrated Member of the MED-Mediterranean Institute for Agriculture, Environment and Development-Research Group of Soil, Water, and Climate (SWC) and the European Turfgrass Society (ETS)-actually he is one of the Board Members.
PEDRO J. S. CARDOSO received the B.Sc. degree in mathematics/computer science from the University of Coimbra, Portugal, the M.Sc. degree in computational mathematics from the University of Minho, Portugal, and the Ph.D. degree in mathematics/operations research from the University of Seville, Spain. He is currently a Full (Coordenador) Professor and a Researcher with the Universidade do Algarve for more than two decades and a member of the Associate Laboratory of Robotics and Engineering Systems (LARSyS). He has broad knowledge in the fields of algorithms and data structures in general, with particular emphasis in the field of machine learning and multiple objective meta-heuristics. Over the past years, he has been involved in over a dozen scientific and development projects, authored more than eighty publications, and an editor of over a dozen of books.