Hybrid Statistical and Machine Learning Methods for Daily Evapotranspiration Modeling

Abstract

Machine learning (ML) models, including artificial neural networks (ANN), generalized neural regression networks (GRNN), and adaptive neuro-fuzzy interface systems (ANFIS), have received considerable attention for their ability to provide accurate predictions in various problem domains. However, these models may produce inconsistent results when solving linear problems. To overcome this limitation, this paper proposes hybridizations of ML and autoregressive integrated moving average (ARIMA) models to provide a more accurate and general forecasting model for evapotranspiration (ET₀). The proposed models are developed and tested using daily ET₀ data collected over 11 years (2010–2020) in the Samsun province of Türkiye. The results show that the ARIMA–GRNN model reduces the root mean square error by 48.38%, the ARIMA–ANFIS model by 8.56%, and the ARIMA–ANN model by 6.74% compared to the traditional ARIMA model. Consequently, the integration of ML with ARIMA models can offer more accurate and dependable prediction of daily ET₀, which can be beneficial for many branches such as agriculture and water management that require dependable ET₀ estimations.

1. Introduction

Efficient irrigation management is a critical aspect of modern agricultural techniques, and accurately estimating evapotranspiration (ET) is essential for effective water resource management, irrigation scheduling, watershed management, and drainage system planning [1]. ET is a method used to measure the water requirements of a crop, which comprises the movement of water vapor from the soil into the air through evaporation from the soil and transpiration from the plants [2]. The first step in determining the ET of an agricultural system is to calculate the reference evapotranspiration (ET₀), which is a widely accepted method for quantifying the water requirements of a crop. However, estimating ET₀ is a complex task that can be conducted either by direct measurement with lysimeters or indirectly through mathematical models.

Typically, lysimeters that provide accurate measurements are used to develop and evaluate other indirect methods [3,4,5]. However, this method is often deemed impractical due to its time-consuming and precise measurement requirements. Moreover, the high cost and complexity associated with lysimeters typically limit their application to research settings. Therefore, for practical purposes, mathematical models based on weather station data have become the preferred alternative [6].

In indirect methods, the equations used to calculate ET₀ values are usually complex, nonlinear, contain random factors, and rely on several assumptions. The literature reports about 20 methods that can be used to estimate ET₀ based on meteorological variables. Of these methods, FAO Penman–Monteith (FAO56PM) [2], Thornthwaite [7], Blaney and Criddle [8], Priestly and Taylor [9], and Hargreaves and Samani [10] have been widely and successfully applied. However, accurate estimates obtained by these methods require the application of rigorous optimization procedures, accurate spatiotemporal data, and a thorough understanding of initial conditions [11].

In recent years, the development of artificial intelligence methods based on machine learning (ML) has gained worldwide attention as an alternative approach to estimate ET₀ using a minimal number of weather input parameters [12]. Numerous research papers have found that ML techniques such as artificial neural networks (ANN), generalized neural regression networks (GRNN), and adaptive neuro-fuzzy interface systems (ANFIS) perform better in ET₀ estimation than empirical and semi-empirical methods [13,14,15,16,17,18,19,20,21,22,23].For instance, Kumar et al. [13] compared the performance of ANN and the Penman–Monteith method in predicting ET₀ and concluded that ANN produced better estimates. Similarly, Ladlani et al. [21] found that the GRNN model outperformed the Priestley–Taylor and Hargreaves–Samani models in ET₀ estimation. Pour et al. [22] used ANFIS and ANN to estimate daily ET₀ with various combinations of independent variables consistent with the Hargreaves–Samani, Priestley–Taylor, and FAO56-PM equations. The results indicated that both ANFIS and ANN provided satisfactory performance in estimating ET₀ from available climate data.

One strategy for achieving performance improvements through ML techniques in estimating ET₀ is to use previous days’ meteorological data as input to the models in addition to the current day’s data. The autoregressive integrated moving average model (ARIMA), commonly known as the Box–Jenkins methodology, is a stochastic model with a strong and reliable performance in predicting a wide range of climatic, hydrologic, and meteorological variables [24,25,26]. However, a single model for ET₀ prediction using time series data is prone to significant errors. Therefore, hybrid approaches that combine multiple models to improve accuracy could be an alternative solution. Arca et al. [27] reported the superiority of the ARIMA model over the ANN model in estimating daily ET₀ in Italy. Kishore et al. [28] employed ANN and ARIMA models to predict ET₀ in Kanchipuram, India, and noticed that the ARIMA model was very effective and reliable for short-term forecasts. Landeras et al. [29] found that applying ARIMA and ANN models improved the performance of predicting weekly ET₀ one week ahead compared to the model based on annual averages.

These studies have contributed significantly to the knowledge base on using ML and ARIMA models in estimating ET₀. However, the applications of combined ML techniques to reduce the errors of the ARIMA model for ET₀ estimation problems are currently limited, and the knowledge on this topic is still incomplete and fragmented. To fill this gap, this study aims to present a new hybrid model that combines the advantages of ARIMA and ML techniques to achieve accurate daily ET₀ prediction. The main contributions of this study are summarized as follows:

-

A time series ARIMA model and various ML techniques, including ANN, GRNN, and ANFIS, are built to predict the daily ET₀ of Samsun, Türkiye, based on climate parameters.

The models’ accuracy, applicability, and reliability are compared to the ASCE Penman–Monteith method.

-

A novel approach combining ARIMA and ML models is developed for ET₀ predictions.

The remainder of this paper is structured as follows. Section 2 details the ARIMA and ML models used for ET₀ estimates. Section 3 provides details on the methodology employed, including the study area, dataset, model selection, and performance criteria for the models. Section 4 compares the predictive performance of the models and presents the numerical results. Section 5 discusses the study’s findings and provides information on the limitations of this study. Finally, Section 6 summarizes the conclusions and offers proposals for future work.

2. Background

2.1. Linear Process: Box–Jenkins ARIMA Model

The Box–Jenkins method is a statistical forecasting technique used to predict and control univariate time series. It was developed assuming that time-dependent events are random and the associated time series are stochastic processes. This method applies to discrete and stationary time series of uniformly distributed observed values [30].

The ARIMA model with linear stochastic structure is commonly used in many research fields for predicting desired data. An ARIMA (p, d, q) model is a combination of autoregression AR (p) (an additive linear function of p past observations), an integer I (d) (d makes a series stationary), and moving average MA (q) (q random error). The partial autocorrelation function (PACF) and autocorrelation function (ACF) of the data are used to determine the order of q and p of the ARIMA model. The model for each environmental factor can be described mathematically in Equation (1).

$$\begin{array}{cc}\hfill X_t^i& ={\alpha }_t^i+{\phi }_1^i{X}_{t-1}^i+{\phi }_2^i{X}_{t-2}^i+\cdots +{\phi }_p^i{X}_{t-p}^i\hfill \\ & -{\epsilon }_t^i-{\theta }_1^i{\epsilon }_{t-1}^i-{\theta }_2^i{\epsilon }_{t-2}^i-\cdots -{\theta }_q^i{\epsilon }_{t-q}^i\hfill \end{array}$$

(1)

where $X_t^i$ is the ith influencing factor resulting from the differentiation of d times; ${\phi }_1^i$, ${\phi }_2^i$, …, ${\phi }_p^i$ and ${\theta }_1^i$, ${\theta }_2^i$, …, ${\theta }_p^i$ are autoregressive; and moving average coefficients have to be calculated for the ith factor. It is assumed that ${\epsilon }_t^i$, ${\epsilon }_{t-1}^i$, ${\epsilon }_{t-2}^i$, …, ${\epsilon }_{t-p}^i$ have a mean of zero and a constant variance.

The fitting of an ARIMA model is conducted in 4 steps as follows:

-

Identification of the ARIMA (p, d, q) model structure;

-

Estimate parameters for ARIMA (p, d, q) model;

-

Check residuals to determine model adequacy;

-

Predict future data from existing data.

2.2. Artificial Neural Networks (ANN)

The Multilayer Perceptron (MLP) is a popular type of ANN that is increasingly being used to handle applications with imprecise data or complex attribute relations. The MLP comprises an input layer, one or more hidden layers, and an output layer, with each neuron having continuous inputs and outputs, a sum input function, and a nonlinear activation function. The number of hidden nodes needs to be determined to optimize the MLP and obtain the best performance. This has been studied in numerous research papers, and it is generally accepted that one hidden layer is sufficient for approximating complex nonlinear functions [31,32].

The MSE function (Equation (2)), a popular measure of error applied to ANN and ML algorithms, was utilized to evaluate the quantitative error of the network. The equation calculates the average difference between the actual outputs (x) and the expected outputs (y) over n inputs.

$$MSE=\frac{1}{n}{\displaystyle \sum _{j=1}^n{\left(x_j-y_j\right)}^2}$$

(2)

The stop criteria were set at MSE = 1 × 10⁻⁶ and 5000 epochs. Comprehensive information on ANN can be found in Refs. [33,34,35].

2.3. Generalized Regression Neural Networks (GRNN)

The GRNN proposed by Specht [36] does not require an iterative training process such as the backpropagation technique. Instead, it computes the function between the input and output vectors from the data itself to obtain a more accurate approximation. As the amount of training data increases, the estimation error decreases with minimal constraints on the function. Similar to other regression techniques, GRNN can be used to determine continuous variables and is based on a basic statistical method called kernel regression.

A GRNN is a powerful tool that consists of four layers: input, pattern, summation, and output. In the input layer, data are taken in by a number of observed parameters equal to the number of input units. The pattern layer contains training patterns, while the summation layer comprises two types of neurons: single-division and summation neurons. The single-division neurons are connected to the pattern layer, and the summation neurons are connected to the output layer. The hidden and output layers use Radial basis and linear activation functions, respectively. Lastly, the output layer normalizes the output set by dividing the output of each S-summation neuron by the output of each D-summation neuron, producing the predicted value Y_i for the unknown input vector x [37] (Equations (3) and (4)):

$$Y_i=\frac{{\displaystyle {\sum }_{i=1}^n{y}_i.\exp \left[-D\left(x,x_i\right)\right]}}{{\displaystyle {\sum }_{i=1}^n\exp \left[-D\left(x,x_i\right)\right]}}$$

(3)

$$D(x,x_i)={\displaystyle \sum _{k=1}^m{\left(\frac{x_i-x_{ik}}{\sigma }\right)}^2}$$

(4)

where $y_i$ represents the weighted connection between the ith pattern layer neuron and the S-summation neuron; n is the training pattern numbers; D is the Gaussian function; m is the number of input vector elements; and x_k and x_ik are the jth elements of x and x_i, respectively. The optimal value of the spread parameter, denoted by σ, is determined experimentally.

2.4. Adaptive Neuro-Fuzzy Inference System (ANFIS)

ANFIS is a powerful tool that combines the learning capabilities of neural networks with the decision-making process of fuzzy logic [38,39]. It constructs an optimized neural network based on a training dataset suitable for the task. The accuracy and applicability of the model to the task can then be evaluated by testing on unseen data. ANFIS has an advantage over other ANN, as it is more interpretable thanks to the inclusion of fuzzy logic. The fuzzy logic can justify the weighting values set for the model, which is a major shortcoming of ANN. The ANFIS architecture consists of five layers: fuzzification, rule, normalization, defuzzification, and summation. The fuzzification layer assigns membership functions to the input values to form fuzzy clusters, which are determined by a set of parameters called premise parameters. The degrees of membership are then calculated based on these parameters, as shown in Equations (5) and (6).

$$O_i^1={\mu }_{A_i}(x)=\frac{1}{1+{\left|\frac{x-c_i}{a_i}\right|}^{2b_i}}$$

(5)

$$O_i^1={\mu }_{B_i}(y)=\frac{1}{1+{\left|\frac{y-c_i}{a_i}\right|}^{2b_i}}$$

(6)

where x and y are the inputs to node i; A is a linguistic label associated with this node function; O_i is the membership function of A_i; and a_i, b_i, and c_i are parameters of the membership function.

The fuzzification layer takes the input data and assigns a membership value to each fuzzy set. They are converted into firing strengths (w_i) for the rules. The membership values can be calculated using Equation (7):

$$O_i^2=w_i={\mu }_{A_i}(x).{\mu }_{B_i}(y)$$

(7)

The normalization layer derives a normalized value of firing strengths for each rule. The normalized value results from the ratio of the firing strength of the ith rule to the sum of all firing strengths (Equation (8)):

$$O_i^3=w_{i,avg}=\frac{w_i}{w_1+w_2+w_3+w_4}\begin{array}{cc},& i\in \left\{1,2,3,4\right\}\end{array}$$

(8)

The weights of the rules determined by using a first-order polynomial in each node of the defuzzification layer are calculated as shown in Equation (9).

$$O_i^4=w_{i,avg}{f}_i=w_{i,avg}(p_{i}x+q_{i}y+r_i)$$

(9)

where w_i,avg is the output of the normalization layer and p, q, and r are the parameter sets. These parameters are referred to as consequence parameters. Each rule has one more consequence parameter than the number input.

The summation of the outputs of the defuzzification layer gives the actual output of the ANFIS (Equation (10)).

$$O_i^5={{\displaystyle {\sum }_i{w}_{i,avg}f}}_i=\frac{{\displaystyle {\sum }_i{w}_i{f}_i}}{{\displaystyle {\sum }_i{w}_i}}$$

(10)

3. Methodology

The hybrid method for predicting daily ET₀ values is a four-step process (Figure 1).

Step 1: The ASCE Penman–Monteith method was used to calculate daily ET₀ values by utilizing meteorological data from the Atakum weather station;

Step 2: The ARIMA, ANN, GRNN, and ANFIS models used ET_0(t−1), ET_0(t−2), and ET_0(t−3) as inputs to predict daily ET₀;

Step 3: The hybrid models, including ARIMA–ANN, ARIMA–GRNN, and ARIMA–ANFIS, used the error values of the periods t − 1, t − 2, and t − 3 as input parameters to predict daily ET₀;

Step 4: Statistical and graphical methods were employed to evaluate the predictive accuracy of all models for ET₀ values.

3.1. Study Area and Dataset

Daily meteorological measurement data for a 10-year period (from January 2010 to January 2020) were obtained from the Atakum Weather Station (41°34′ north latitude, 36°25′ east longitude, and 4 m altitude), operated by the General Directorate of Meteorology of Samsun province in Türkiye. The city of Samsun is situated between the deltas of the Yeşilırmak and Kızılırmak rivers, covering an area of 9083 km². The climate in Samsun is generally temperate, with distinct characteristics on the coast and in the interior regions. The effects of the Black Sea climate are felt on the coast, leading to hot summers and rainy winters. Inland, the city is influenced by the Akdağ and Canik Mountains, with cold, rainy, and snowy winters and cool summers [40]. The average annual temperature in Samsun is 14.6 °C, which is higher than the Turkish average of 13.6 °C. The average minimum temperature is 3.9 °C in February, while the average maximum temperature is 27.1 °C in August. The ASCE Penman–Monteith method was used to calculate daily ET₀ using meteorological data from the Atakum weather station.

3.2. Data Pre-Processing

The process of data preprocessing is critical to achieve higher performance and accuracy for ML models. It includes noisy, missing, and inconsistent data. In this study, daily weather data (temperature, relative humidity, wind speed, and sunshine duration) were merged as inputs with calculated ET₀ data. Further preprocessing of the data included data cleaning, transformation, and splitting. Standardization of the data to a range of 0 and 1 was implemented to overcome the challenges associated with processing data with different units of measurement and to avoid complications associated with extreme values. To evaluate the predictive performance of the different models, each dataset was divided into two subsets: a training set and a testing set. The training set, comprising 70% of the data, was used exclusively for model development. After model development, the testing set, comprising 30% of the data, was used to evaluate the proposed models.

3.3. Selection of Component Models

A number of linear statistical models have been proposed in the literature for time series models, which are generally divided into two categories: linear and nonlinear models. ARIMA models are one of the most popular linear time series models widely used in hydrology and meteorology [41,42,43,44,45]. ARIMA is popular for its statistical properties and the well-known Box–Jenkins method [30]. By using a differencing process, the ARIMA model can effectively convert non-stationary data into stationary data, so it can be used to fit non-stationary time series. Therefore, in this paper, the ARIMA model is chosen for developing the proposed hybrid approach. ANN, GRNN, and ANFIS are employed to identify the nonlinear component of the time series, based on their successful applications in modeling time series data [46,47,48,49].

3.4. Performance Criteria of Model

To accurately assess the performance of the models in this study, three commonly used metrics were employed: the coefficient of determination (R²), root mean square error (RMSE), and mean absolute error (MAE). These equations, outlined in Equations (11)–(13) by Waller [50], provide a comprehensive and quantifiable assessment of the models’ performances.

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^n{\left(Z_{i,m}-Z_{i,p}\right)}^2}}{{\displaystyle \sum _{i=1}^n{\left(Z_{i,m}-Z_{i,avg}\right)}^2}}$$

(11)

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(Z_{i,m}-Z_{i,p}\right)}^2}}{n}}$$

(12)

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|Z_{i,m}-Z_{i,p}\right|}$$

(13)

where Z_i,m is the measured value; Z_i,p is the predicted value; Z_i,avg is the average value measured; and n is the number of data.

Additionally, a Taylor diagram was used to analyze the standard deviation (SD) and correlation coefficients (R) between the predicted and observed values. This analysis provided further insight into how closely the model predictions matched the observed data.

4. Results

In this paper, ET₀ estimation was performed in two stages. In the first stage, ET_0(t−1), ET_0(t−2), and ET_0(t−3) were taken as input parameters for predicting daily ET₀ using the ARIMA, ANN, GRNN, and ANFIS models. In the second stage of the study, the error values of the periods t − 1, t − 2, and t − 3 were used as input parameters. The difference between the results of the ARIMA model and the observed ET₀ data served as the output value. Following that, the ARIMA combination value was calculated with re-estimated errors, considering the discrepancies between the errors estimated by ANN, GRNN, and ANFIS, the initial ARIMA, and the actual ET₀ values.

The implementation of the ARIMA model is based on a time series and should be stationary. It is influenced by climate, soil, and crop conditions that affect the daily ET₀ value. The stationarity of the daily ET₀ was studied considering the factors affecting crop water use. The distributed lag model is a time series model in which the effect of a single predictor, such as the lag of the predictor itself, on the response variable can vary over time. The lag length, or the time period that best explains the effects of the predictor and its lagged values on the response variable, can be determined by examining the model’s fit to the data and measuring the model’s predictive accuracy. The lag length can be determined by examining the model’s shape and computing the Akaike Information Criterion (AIC). In this study, the quadratic polynomial equation as a function of lag length was used.

The Box–Jenkins technique was employed to analyze the autocorrelation function (ACF) and partial autocorrelation function (PACF) of ET₀ values for the ARIMA model. The ACF and PACF plots indicated that the time series of ET₀ was stationary (Figure 2a,b). Therefore, the ARIMA (1, 0, 1) model was implemented for the ET₀ estimation. The coefficients of AR1 and MA1 parameters were statistically significant (p < 0.01) based on the results of Maximum Likelihood Estimation (MLE), as shown in

Table 1. AR1 and MA1 coefficients of final ARIMA equation.

Parameters	Coefficient	Standard Error	t-Statistic	p-Value
Constant	0.12072	0.00820	14.72	0.000
AR1	0.95476	0.00671	142.37	0.000
MA1	0.46610	0.01990	23.41	0.000

.

In the ANN technique, a feedforward backpropagation MLP with a single hidden layer was utilized, and the Levenberg–Marquardt training algorithm was used to train the network. The hidden and output layers were constructed using the tangent sigmoid (tansig) and linear transfer (purelin) functions, respectively. The number of hidden nodes was optimized by incrementally increasing it from three to seven, and the best-performing model was selected based on the lowest MSE value in the testing phase. The ANN (3, 5, 1) model provided the best accuracy during the test period compared to the other ANN models based on the R², RMSE, and MAE criteria. In this model, three represents the input variables, namely, ET_0(t−1), ET_0(t−2), and ET_0(t−3)), five is the number of hidden layer neurons, and one is the output variable (ET₀).

In the GRNN technique, various spread factors (from 0.1 to 1) were tested to find the optimal model using a trial-and-error approach. The GRNN (3, 0.1) model with inputs of ET_0(t−1), ET_0(t−2), and ET_0(t−3) and a spread value of 0.1 demonstrated higher accuracy than other GRNN models during the testing phase.

To evaluate the performance of the ANFIS technique, different types of membership functions (MFs), such as Gaussian, trapezoidal, and triangular, were tested with varying numbers of MFs. The ANFIS model with three triangular MFs for the inputs ET_0(t−1), ET_0(t−2), and ET_0(t−3) was found to be the most effective among all the ANFIS models during the testing phase.

Table 2. Comparison of ARIMA, ANN, GRNN, and ANFIS for daily ET₀ prediction using training and testing datasets.

Model	Model Structure	Training			Testing
		MAE	RMSE	R2	MAE	RMSE	R2
ARIMA	(1, 0, 1)	0.562	0.778	0.719	0.580	0.771	0.748
ANN	(3, 5, 1)	0.548	0.759	0.732	0.563	0.755	0.760
GRNN	(3, 0.1)	0.511	0.702	0.938	0.527	0.697	0.946
ANFIS	(3, trimf)	0.531	0.736	0.749	0.558	0.762	0.774

summarizes the training and testing performance of the best ARIMA, ANN, GRNN, and ANFIS models. The results indicate that the GRNN model achieved the highest accuracy during the test phase (R² = 0.946, RMSE = 0.697 mm, and MAE = 0.527 mm) compared to the other models, followed by ANFIS, ANN, and ARIMA. Figure 3 shows the scatter plots of these models for the observed and predicted daily ET₀ values in both training and testing periods.

The second stage of the study aimed to construct a combined model, which utilizes a residual from the ARIMA model as input to a new ML network. This is due to the nonlinear aspect of the ARIMA model, which allows for more complex features to be incorporated into the model. The training and testing results of the models obtained with the hybrid approach are shown in

Table 3. Comparison of combined approaches (ARIMA–ANN, ARIMA–GRNN, and ARIMA–ANFIS) for daily ET₀ prediction using training and testing datasets.

Model	Model Structure	Training			Testing
		MAE	RMSE	R2	MAE	RMSE	R2
ARIMA–ANN	(3, 5, 1)	0.521	0.715	0.763	0.541	0.719	0.782
ARIMA–GRNN	(3, 0.1)	0.254	0.381	0.934	0.269	0.398	0.935
ARIMA–ANFIS	(3, trimf)	0.508	0.700	0.773	0.534	0.705	0.791

. As seen from the table, the ARIMA–GRNN hybrid model achieved the lowest RMSE (0.400 mm) and MAE (0.270 mm) values in the testing stage, demonstrating its efficiency in predicting ET₀ values. Moreover, hybrid approaches were the most accurate overall, as evidenced by the statistical indices. The scatter plots of these hybrid models for the observed and predicted daily ET₀ values in the training and testing periods are displayed in Figure 4.

The hybrid models (ARIMA–GRNN, ARIMA–ANFIS, and ARIMA–ANN) were able to outperform both the traditional ML models (GRNN, ANFIS, and ANN) and the traditional ARIMA model by significantly reducing RMSE values. Specifically, the ARIMA–GRNN model reduced the RMSE by 48.38%, the ARIMA–ANFIS model by 8.56%, and the ARIMA–ANN model by 6.74%, compared with the traditional ARIMA model. Similarly, the ARIMA–GRNN model reduced RMSE by 42.90%, the ARIMA–ANFIS model by 7.48%, and the ARIMA–ANN model by 4.77% compared with the non-ARIMA models.

The Taylor diagram, presented in Figure 5, provides an evaluation of the different model’s performance in terms of bias, consistency, and scatter. The black point on the diagram represents the observed ET₀ value, while the other points correspond to the models’ predictions. The ARIMA (R = 0.87 and SD = 1.32 mm) and ANN (R = 0.87 and SD = 1.33 mm) models are significantly far apart, indicating their poor performance in estimating ET₀ values. On the other hand, the ANFIS (R = 0.88 and SD = 1.35 mm), ARIMA–ANFIS (R = 0.89 and SD = 1.40 mm), and ARIMA–ANN (R = 0.88 and SD = 1.39 mm) models are extremely close, with no noticeable difference between them. Lastly, the GRNN (R = 0.97 and SD = 1.46 mm) and ARIMA–GRNN (R = 0.97 and SD = 1.97 mm) models are the closest to the calculated ET₀ results.

5. Discussion

Both time series and ML models have been independently evaluated for their ability to predict ET₀ values. Time series models are useful because they account for trends in the data and can be used to predict future values. ML models have also been shown to be effective in predicting ET₀ values because they can better identify more complex patterns in the data [15].

To achieve accurate prediction of ET₀ values, a combination of time series and ML models can be used. This approach needs to be customized to the specific dataset to yield the best outcomes. Firstly, the time series model should be employed to identify any existing trends in the data. These trends are then used as the basis for the ML model. The ML model should then refine the predictions further by taking into consideration any other factors that could affect the outcome.

This paper proposed hybridizations of ML and ARIMA models to create a more general and accurate ET₀ forecasting model. Among the models used in this study, the GRNN and ARIMA–GRNN models had lower RMSE and MAE values, making them slightly better than the others. Previous studies have independently evaluated the performance of time series [51,52,53], ML techniques [54,55], and hybrid models [56,57] for predicting ET₀ rates. For example, Ferreira et al. [58] utilized an ANN model to estimate daily ET₀, yielding RMSE of 0.70 mm and R² of 0.82. This study produced similar results, with RMSE of 0.755 mm and R² of 0.760. Antonopoulos and Antonopoulos [59] found that ANN models could estimate daily ET₀ with an accuracy ranging from an RMSE of 0.574 to 1.33 mm. In the study by Pour-Ali Baba et al. [22], the ANFIS models showed good performance, with RMSE values ranging from 0.474 to 0.851 mm. Dou and Yang [60] suggested that the ANFIS method is an valuable complement to traditional methods for estimating ET due to its robustness and flexibility. Landeras et al. [29] found that ARIMA and ANN models were more accurate than a mean year model based on historical averages, with a significant improvement in accuracy of 6–8% in reducing the mean squared differences of forecasts. These findings are consistent with the results of the present study. Arca et al. [27] reported the superiority of the ARIMA model in predicting daily ET₀ rates over the ANN model, while e Lucas et al. [60] found no significant difference between ARIMA and ML models in predicting ET₀ rates.

Although this study has yielded significant findings, there are some limitations that should be taken into consideration. For instance, the study only used data from a single location (Samsun province of Türkiye) and for a limited period (2010–2020). Additionally, the study only considered three ML methods (ANN, GRNN, and ANFIS) and ARIMA models, and other variables may also play an important role. In future studies, Deep Learning models with hybrid techniques need to be explored to improve the accuracy of ET₀ prediction further.

6. Conclusions

A new combined ARIMA soft computing modeling approach is proposed for estimating the daily water consumption of plants. Time series data were used to estimate ET₀ in Samsun province, northern Türkiye. Before applying ARIMA and ML models, the data stationarity was checked, and then the best ARIMA and ML models were developed from these stationary data. The behavior of ACF and PACF criteria was analyzed to determine the appropriate ARIMA model.

The hybrid methods, which included ARIMA–ANN, ARIMA–GRNN, and ARIMA–ANFIS, showed the best accuracy compared to individual ARIMA model, as indicated by statistical performance criteria (RMSE, MAE, R²). The ARIMA–GRNN model showed the highest improvement compared to ARIMA with a 48.38% decrease in RMSE values. The results of the hybrid models were satisfactory, as demonstrated by the comparison between observed and predicted values using different methods.

The hybrid approach has several advantages by combining the strengths of ARIMA and ML models, resulting in more accurate predictions. It also simplifies the process and allows for faster decision making. The findings of this study are noteworthy; however, it has limitations in terms of the data (2010–2020) and models used (ARIMA, ANN, ANFIS, and GRNN). Therefore, future studies could further improve the accuracy of ET₀ predictions by incorporating different hybrid techniques and deep learning models.