Comparison of machine learning techniques and spatial distribution of daily reference evapotranspiration in Türkiye

Reference evapotranspiration (ET0) estimates are commonly used in hydrologic planning for water resources and agricultural applications. Last 2 decades, machine learning (ML) techniques have enabled scientists to develop powerful tools to study ET0 patterns in the ecosystem. This study investigated the feasibility and effectiveness of three ML techniques, including the k-nearest neighbor algorithm, multigene genetic programming, and support vector regression (SVR), to estimate daily ET0 in Türkiye. In addition, different interpolation techniques, including ordinary kriging (OK), co-kriging, inverse distance weighted, and radial basis function, were compared to develop the most appropriate ET0 maps for Türkiye. All developed models were evaluated according to the performance indices such as coefficient of determination (R2), root mean square error (RMSE), and mean absolute error (MAE). Taylor, violin, and scatter plots were also generated. Among the applied ML models, the SVR model provided the best results in determining ET0 with the performance indices of R2 = 0.961, RMSE = 0.327 mm, and MAE = 0.232 mm. The SVR model’s input variables were selected as solar radiation, temperature, and relative humidity. Similarly, the maps of the spatial distribution of ET0 were produced with the OK interpolation method, which provided the best estimates.


Introduction
Accurate estimation of crop water consumption (evapotranspiration, ET) and reference evapotranspiration (ET 0 ) is critical for water resource management, especially in arid and semiarid regions (Banda et al. 2018;Tunca et al. 2022). ET 0 is a complicated hydrologic cycle component crucial to agricultural water management. Knowledge of ET 0 would make it possible to reduce water consumption, increase irrigation efficiency, and allow properly scheduling of the irrigation application (Sattari et al. 2021). ET 0 can be determined by conducting laboratory experiments, or it can be estimated using mathematical models. Lysimeters, which provide accurate measurements of ET 0 , are typically used to develop and test other indirect ET 0 measurement methods (López-Urrea et al. 2006). However, it is not always possible to measure ET 0 with lysimeters because this method is costly, labor-intensive, and requires careful measurements (Allen et al. 1998). Therefore, using mathematical models by collecting meteorological data from the weather stations is preferable to estimate ET 0 (Ferreira et al. 2019).
In indirect methods, ET 0 is calculated using equations calibrated with direct methods; ET 0 is then converted to standard plant ET using the crop coefficient for different growing seasons. Generally, the equations used to calculate ET 0 values are complex, nonlinear, contain random factors, and rely on several assumptions. According to the literature, there are about 20 methods for estimating ET 0 based on meteorological variables. Of the methods, Penman (Penman 1948), Thornthwaite (Thornthwaite 1948), Blaney and Criddle (Blaney and Criddle 1950), Priestly and Taylor (Priestley and Taylor 1972), Hargreaves and Samani (Hargreaves and Samani 1985), and FAO Penman-Monteith (FAO56PM) (Allen et al. 1998) have been widely and successfully used mathematical models to estimate ET 0 . For these mathematical methods to solve the equations, rigorous optimization procedures, accurate spatiotemporal data, and knowledge of initial conditions are required (Prasad et al. 2017).
On the other hand, with recent developments in computer technology, researchers worldwide have gained attention from machine learning (ML) techniques. ML represents input-output relations without understanding the physical process, making them effective tools for modeling nonlinear systems. They have provided an alternative to estimating ET 0 based on a minimum number of weather input parameters (Krishnashetty et al. 2021). Recently, researchers have investigated the potential of various ML techniques, such as artificial neural networks (ANN), support vector machine (SVM), and genetic programming (GP), to determine ET 0 (Abrishami et al. 2019;Mohammadrezapour et al. 2019;Chia et al. 2020;Ayaz et al. 2021;Jang et al. 2021;Achite et al. 2022;Dimitriadou and Nikolakopoulos 2022;Makwana et al. 2022;Tejada et al. 2022). It is worth noting that the inputs to these models lend themselves well to sensitivity analyses and obtaining optimal structures and that they do not rely on mathematical relationships even for complex phenomena. Although these models have been trained and tested in the literature for many meteorological stations and climatic conditions, it is impractical to train them for every station in a large region. Therefore, regional models need to be created (Citakoglu et al. 2014).
Accurate field-level assessments of ET 0 can be critical for understanding the spatiotemporal patterns of water demand and improving water resource management and agricultural practices, such as decisions about the timing and quantity of irrigation (Kustas et al. 2018). The spatial variation of ET 0 is important for climate analyses and assessment of regional climate change scenarios such as droughts and floods (Güler 2014). In recent years, the geographic information system (GIS) has been extensively used to evaluate the spatial and temporal variation of ET 0 (Mardikis et al. 2005;Vicente-Serrano et al. 2007;Sabziparvar and Tabari 2010). Using GIS tools, spatial interpolation techniques can be applied to create ET 0 maps by expanding ET 0 values from points to areal/regional estimates. Interpolation techniques generally fall into two main categories: deterministic and geostatistical (stochastic). Deterministic interpolation generates surfaces using mathematical functions from sample points based on their similarity as inverse distance weighting (IDW) or smoothing as radial basis functions (RBF). In contrast, stochastic methods (kriging) use mathematical and statistical functions to evaluate the uncertainty of estimates (Webster and Oliver 2007). These interpolation methods are easily computerized and can be combined with contour mapping of regional variables. Over the past 2 decades, a number of spatial interpolation methods for ET 0 mapping have been tested with varying degrees of complexity. For instance, da Silva Júnior et al. (2019) evaluated methods for spatial interpolation of ET 0 data in terms of precision and performance. Gentilucci et al. (2021) used geostatistical methods to calculate ET 0 and calibrate the Hargreaves equation over the last 10 years in Italy. Dalezios et al. (2002) studied the optimal interpolation of ET 0 over large regions of complex terrain in Greece. Hodam et al. (2017) investigated the potential of IDW and kriging to develop an appropriate method for spatial interpolation of point ET 0 .
The studies above have contributed significantly to the knowledge base on using ML and interpolation techniques in estimating ET 0 . However, the combined application of ML and different interpolation techniques for ET 0 estimation problems is currently limited, and the knowledge on this topic is still incomplete and fragmented. To fill this gap, this current study aims to investigate the performance of ML and interpolation techniques and to determine the best model for accurately predicting daily ET 0 values. The main contributions of this study are as follows: The potential estimation of daily ET 0 with k-nearest neighbor algorithm (KNN), support vector regression (SVR), and multigene genetic programming (MGGP) are investigated. The accuracy, applicability, and reliability of the models are examined in comparison to the ASCE Penman-Monteith (ASCE-PM) method. Various interpolation techniques, including ordinary kriging (OK), co-kriging (COK), inverse distance weighted (IDW), and radial basis function (RBF), are compared to find the most suitable spatial maps of ET 0 .

Research site and data sets
Türkiye has located in the northern hemisphere and is bounded by latitudes of 36-42° N and longitudes of 26-45° E. The country is closer to the equator than the poles and is in a temperate climate zone. It is approximately 1,600 km long and 800 km wide, with a total geographic area of 783,562 km 2 . The total water area is 9820 km 2 . The lowest altitude is 0 m at the Mediterranean Sea, and the highest altitude is 5137 m at Mount Ararat. Daily meteorological measurement data for a 52-year period (between 1962 and 2014) were collected from 213 meteorological stations (Fig. 1). The ASCE Penman-Monteith method was used to calculate daily ET 0 using meteorological data such as temperature (T), relative humidity (RH), wind speed (W S ), and solar radiation (R S ) from these meteorological stations. Information on the daily ET 0 values and the locations of the stations representing different geographic regions and climatic conditions are presented in Table 1.

Method of ET 0 Estimation
To compute ET 0 on a daily scale, the ASCE-PM equation (Allen et al. 2006) (Eq. 1) was employed.
where ET 0 is the reference evapotranspiration (mm d −1 ), R n is net solar radiation reaching the soil surface (MJ m −2 d −1 ), G is soil heat flux density (MJ m −2 d −1 ), T is mean daily air temperature at the height of 2 m (°C), u 2 is the wind speed at the height of 2 m (m s −1 ), e s is saturated vapor pressure (kPa), e a is actual vapor pressure (kPa), Δ is the slope of vapor pressure curve (kPa °C −1 ), and γ is psychometric constant (kPa °C −1 ).

k-nearest neighbor algorithm (KNN)
In many cases, the KNN is an effective nonparametric ML approach for classification and regression (Guo et al. 2003;Cheng et al. 2022). Since the KNN is a distance-based model, the basic logic of the KNN algorithm is to determine the correct class for the test data by calculating the distance between the test data and all training points. In the method, the letter "k" stands for the number of neighbors, and the success of the classification depends heavily on this value. A simple way to determine the "k" value is to run the algorithm repeatedly with different values, select the best performance on the training data, and then classify the unknowns (Guo et al. 2003). The majority of its neighbors classify data, and its class is determined by the frequency of its "k" nearest neighbors, measured by a distance function. Although there are many distance functions, such as Euclidean distance, Mahalanobis distance, and Minkowsky distance and their variants, the Euclidean distance function has been used most frequently in the literature, and therefore this function was preferred in the present study (Mittal et al. 2019;Prasad et al. 2019;Xu et al. 2020;Bayram and Çıtakoğlu 2023). The Euclidean distance function can be expressed as shown in Eq. (2).
where x i and x j are the ith and jth variables, and n is the data number. KNN has a major disadvantage that a comprehensive record search is required for each sample to identify the most relevant samples (Sahoo and Ghose 2022). On the other hand, the KNN performs classification in a simple, fast, easy-to-understand, and competitive way (Imandoust and Bolandraftar 2013).

Support vector regression (SVR)
SVR is widely used in linear and nonlinear regression for prediction and curve fitting. SVR is based on the elements of the SVM, where support vectors are closer points in the direction of the generated hyperplane in an n-dimensional feature space that separates the data points around the hyperplane. SVR uses a part of the given data set to create a function estimator (Eq. 3).
where w ∈ R n is a weighted feature vector, and x is the input vector in R n . ϕ represents the mapping, and b is the intercept. This study uses ε-SVR because the ET 0 prediction model does not limit the number of support vectors (Ceperic et al. 2013). The optimization problem for the ε-insensitive problem is formulated using Eq. (4).
in which ½║ω║ 2 indicates the parameter for regularization, C is the empirical error penalty factor, ξ i and ξ i * are the slack variables, and ε represents the loss function.
To solve this optimization problem, a Lagrangian dual function can be constructed. The performance of SVR depends heavily on the kernel functions since SVR is a kernel-based algorithm. Although there are several kernel functions, the RBF is one of the most recommended kernels in the literature (Hosseinzadeh et al. 2021). Its formulation is as defined in Eq. (5).
where ║x − y║ 2 is identified as the squared Euclidean distance between the two feature vectors, and σ is a free parameter.
Despite its high accuracy and excellent generalizability, the main shortcoming of the SVR is that the kernel function is used only once to map the sample data to the highdimensional features (Zhong et al. 2019).

Multigene genetic programming (MGGP)
GP is a supervised ML technique that searches a program space rather than a data space. Symbolic regression is performed over GP to develop trees. GP solves complex problems by constructing and modifying function trees (Gandomi et al. 2010). This technique has the advantage that one can create prediction equations without assuming existing relations. MGGP is a robust variant of the GP model designed to overcome shortcomings. MGGP produces mathematical "multigene" models of predictor response data, i.e., linear combinations of low-order nonlinear input variable transformations. Unlike MGGP, the traditional GP evaluates a single tree expression. In the MGGP model, the input-output relation is determined by a random selection of functions and variables (Niazkar and Niazkar 2020). This model uses one or more gene trees and calibrates the coefficients using least-squares regression analysis (LSRA), a statistical method. Based on the MGGP, evolutionary processes are referred to as highlevel transitions and mutations, as opposed to the processes described in the standard GP, which are referred to as lowlevel processes (Bayram and Çıtakoğlu 2023). In MGGP, linearly combining nonlinear terms without functional structure generates a set of equations. It also provides a mathematical representation of the problem in the form of a tree structure or equations that are formal and explicit . This study used the open-source code of MGGP, which was obtained from the literature (Searson et al. 2010). For the MGGP model, the basic mathematical functions, such as " + , − , × , ÷ , √, x 2 , exp, ln." are used to obtain the optimum models. In the MGGP method, two important control parameters (multigene and tree-building parameters) are set by the user. The accuracy of the model can be increased by increasing the values of these two parameters. Detailed explanations of the MGGP parameters and operators can be found in the literature (Searson et al. 2010;Gandomi and Alavi 2012).

Spatial interpolation methods
Interpolation techniques can be divided into two main categories: deterministic and geostatistical (stochastic). Deterministic interpolation techniques use mathematical functions to generate surfaces from sample points based on the degree of similarity (inverse distance weighting, IDW) or smoothing (radial basis functions, RBF). On the other hand, stochastic techniques (kriging) use mathematical and statistical features to estimate uncertainty. Geostatistical methods are known as kriging. Semi-variograms are used in interpolation by kriging. OK is a common method used in geostatistical studies to predict the value of a variable at a single point or block. It is used when there are no trends in the data. The most common functions used to model semi-variograms are exponential, spherical, and Gaussian functions (Hodam et al. 2017;Küçüktopcu and Cemek 2022). COK is advantageous compared to kriging when working with variables that are expensive to sample (Orejuela et al. 2021). The deterministic methods, on the other hand, use mathematical functions for interpolation. IDW is used to estimate the value at unsampled points from the values of sampled points using linear combinations of values weighted by the inverse function of the distances between the points of interest and the surrounding known values. Some basic parameters are considered when applying the interpolation methods. The accuracy of the method is affected by exponential parameters (power 1, power 2, etc.) (Burrough and McDonnell 1998). The RBF methods are a series of exact interpolation methods in which the surface must pass through each measured sample value. The RBF method has five basic functions: thin plate spline (TSP), spline with tension (ST), completely regularized spline (CRS), multiquadric function (MQ), and inverse multiquadric function (IMQ) (Xie et al. 2011). The flowchart for selected ML methods and interpolation techniques of this study is presented in Fig. 2.

Statistical evaluation
Various statistical performance indicators, such as the mean absolute error (MAE), root mean square error (RMSE), and coefficient of determination (R 2 ), were considered (Cemek et al. 2022) and formulated in Eqs. (6-8).
in which Z i and Z p are the values calculated and predicted at the jth observation, Z i,avg is the mean value of the calculated variable, and n is the data number.

Graphical evaluation
For this study, in addition to the statistical analyses described above, we employed three different visualization methods, including scatter, violin, and Taylor plots to evaluate the effectiveness of the models examined.

Machine learning results
In this study, a Pearson correlation analysis was performed between meteorological input variables including T, RH, R S , and W S , and the output variable (ET 0 ). According to the correlation analysis, R S has the strongest effect on ET 0 , while W S has the weakest one (Fig. 3). Previous studies also support this argument (Gong et al. 2021;Ge et al. 2022). The highest linear relationships were found between R S , T, RH, and ET 0 ; therefore, R S , T, and RH were used as the input parameters. Analyses were performed with the following seven combinations: (i) R S , (ii) T, (iii) RH, (iv) R S and T, (v) R S and RH, (vi) T and RH, and (vii) R S , T, and RH.
To eliminate unit problems, the data were normalized between 0 and 1. Then, the data were randomly divided into two groups: training (70%) and testing (30%). After the hyperparameters of the models were determined using the training data, the estimation performances of the models were obtained using the test data. For KNN, "number of neighbors (k)," "leaf size (ls)," and "power parameter (p)" are the main hyperparameters for model generation. A grid search procedure was used for each model to find the best values for these parameters. For the parameters k, ls, and p, the ranges 1-50, 1-20, and 1-10 were considered, respectively. In all models examined, the parameters ls and p did not significantly improve model performance, so default values of 30 and 2, respectively, were chosen. On the other hand, the parameter k significantly improved the accuracy of the prediction results. For example, in the model with R S , T, and RH as inputs, k equal to 9 gave the best estimation result (Appendix 1a). All KNN models' prediction performances are also presented in Table 2. In the testing stage, the lowest RMSE and MAE values (0.424-0.309 mm) were obtained for the KNN7 model where the inputs were R S , T, and RH; however, the highest values (1.284-1.015 mm) were obtained for the KNN3 model in which the input was RH.
The three parameters of the SVR, including the error term penalty parameter (C), radius (ε), and kernel coefficient (γ), are called the primary hyperparameters. According to Chang and Lin (2011), the parameter γ is equal to 1/K, where K is the number of inputs. The tuning parameters in the SVR model were determined using a grid search approach with tenfold cross-validation in this paper (Küçüktopcu 2023). For example, the best accuracy was achieved when C and ε in the SVR7 model were 500 and 0.0001, respectively (Appendix 1b). The predictive performances of all SVR models are also shown in Table 3. The SVR7 model achieved the best results (RMSE = 0.327 mm, MAE = 0.232 mm), while the worst estimate was obtained by the SVR3 model (RMSE = 1.237 mm, MAE = 1.029 mm) in the testing phase.  With the MGGP, input and output data can be expressed nonlinearly. MGGP equations and performance criteria for different input models are presented in Table 4. The best results (RMSE = 0.391 mm, MAE = 0.298 mm) were obtained with the MGGP7 model, whereas the worst estimates were obtained with the MGGP3 model (RMSE = 1.797 mm, MAE = 1.535 mm) in the testing phase.
The comparisons of the three models of the best combination of input data (R S , T, RH) for the testing data were presented on the scatter, Taylor, and violin plots (Figs. 4 and 5). In scatter plots (Fig. 4a, b, and c), most data were close to the 1:1 line. As shown in Fig. 4, all three models proved successful in estimating ET 0 values. Among the models examined, the SVR7 model performed best with respect to the R 2 , RMSE, and MAE criteria. In contrast, the KNN7 model lagged behind the other models (SVR7 and MGGP7) in estimation performance. Figure 5a shows the violin plots for each model. The SVR7, KNN7, and MGGP7 models produced similar plots to the calculated model. As shown in Fig. 5a, the SVR7 model provided similar estimates of ET 0 as the calculated values. The Taylor diagram, which includes statistical analysis, was also used in comparing the three ML techniques. The Taylor diagram evaluates the agreement of the estimated data with the reference data. By using the Taylor diagram, further comparisons of the models were obtained. As shown in Fig. 5b, the SVR7 model gave a closer estimate of the calculated value than the MGGP7 and KNN7 models.

Interpolation method results
Descriptive statistics for calculated ET 0 values with the ASCE Penman-Monteith method are shown in Table 5. The highest ET 0 value (5.37 mm day −1 ) was calculated for July and the lowest (0.85 mm day −1 ) for January. The coefficient of variation (CV) measures the variation of an attribute. A CV ≤ 15% indicates a low variation, 16-35% a moderate variation, and a CV ≥ 36% high variation ).
Our CV values revealed that reference evapotranspiration values exhibited a low variation in May (16.10 mm) and moderate variation in January (33.85 mm). The highest coefficient of skewness for ET 0 occurred in November (1.28), January (1.01), and February (1.12), while the distribution was normal in the remaining months. Geostatistical methods are more reliable when the data exhibit a normal distribution. Therefore, data normality was checked before the initiation of geostatistical methods. In most of the cases, data showed a lognormal distribution. Semi-variogram models and model parameters of nugget (Co), Sill (Co + C), and range for daily ET 0 values are shown in Appendix 2. Besides, R 2 and RSS (residual sum of squares) indicated the model's availability, and nugget ratios indicating spatial dependence levels were also determined. Semi-variograms for ET 0 revealed that the exponential model was suitable for January, February, April, November, and December, and the spherical model was suitable for the other months. The highest r (0.98) and the lowest RSS (1.87 × 10 -5 ) values of the spherical model were observed in March. The highest range value (31.41 km) occurred in December and the lowest (2.87 km) in April. Accurate estimation of variables depends on the density of observed data points, their spatial variation, and proper modeling of semivariograms (Fotheringham et al. 2000).
The ratio of nugget semi-variance to total semi-variance is used to classify spatial dependence. The spatial dependence is classified as strongly spatially dependent when the ratio is ≤ 25%, moderately spatially dependent when 25-75%, and weakly spatially dependent when ≥ 75% (Cambardella et al. 1994). Experimental and theoretical semi-variograms for daily ET 0 values are presented in Appendix 3. Spatial Fig. 4 Scatter charts of a SVR7, b KNN7, and c MGGP7 models variation in short distances is generally related to the sampling scheme not capturing the short-distance spatial variation and/or experimental errors. Rehman and Ghori (2000) modeled spatial dependency of global solar radiation by geostatistical techniques where they fitted a spherical model to experimental monthly semi-variograms. Irmak et al. (2010) used predicted minimum and maximum temperatures using their monthly semi-variograms of maximum and minimum temperature in OK.
In the OK method, the Gaussian model outperformed others based on RMSE and MAE. Elevation was used as an auxiliary variable in COK, and the lowest RMSE and MAE values were with the exponential model. In IDW, the lowest RMSE and MAE values were found with linear interpolation (power = 1), and in the RBF method, where four different functions were evaluated, the lowest RMSE and MAE values were observed in the ST model (Appendix 4). Xu et al. (2006) interpolated seasonal and annual ET 0 calculated with the Penman-Monteith method and pan evaporation and concluded that kriging and IDW were the best methods.
Nugget variance is generally high, suggesting that there were fine-scale discontinuities in the ET 0 , which would be resulted from several factors, including the combined effect of differences in climate and topography, small-scale differences that the current sampling scheme could not capture, errors resulted from locations of climate stations, and errors resulted from calculations of ET 0 from climate data.
A greater nugget effect occurred for winter months, indicating a greater short-range variation of ET 0 in the winter. The range values for winter months were far greater than the remaining nine months, which showed that ET 0 was more spatially dependent in the winter than in other seasons. Sill  values representing overall variation in ET 0 were greater for winter months than the rest; this revealed that ET 0 is more variable in winter than in the other seasons. All these showed that the spatial variation of ET 0 in winter was considerably different from that of the rest of the year. Spatial distribution maps for ET 0 predicted by OK, COK, IDW, and RBF were presented in Fig. 6 for January (the coldest month in winter) and July (the hottest month in summer). In each case, the range of ET 0 was divided into ten classes. The spatial pattern of ET 0 predicted four techniques differed across the nation. Although the overall spatial patterns were similar, some significant differences occurred among predictions at some locations. Geostatistical techniques yielded similar spatial patterns, and similarly, IDW and RBF predicted similar spatial patterns of ET 0 across Türkiye in January and July. Figure 6 further shows that the acreage of ET 0 classes differs considerably among prediction methods in both cases (January and July). The validation results showed that OK outperformed the other three prediction methods, justifying its use for predicting ET 0 in Türkiye. Kriging techniques have advantages, including error maps and avoiding data clustering, and the basis provides for stochastic simulation for the possible realization of estimate variables (Karamouz et al. 2012;Geleta et al. 2019). It is not always likely to obtain data from every meteorological station. Therefore, unknown values can be estimated from data taken at missing data locations that can be sampled. Interpolation methods are important in bridging such gaps through spatial mapping of regionalized variables for ET 0 .
The spatial pattern of ET 0 was similar in December and January; March, April, and May; July and August; and September and October (Fig. 7) by the OK-Gaussian method. ET 0 values predicted by the OK agreed well with those calculated by the ASCE Penman-Monteith method. The OK under-predicted ET 0 consistently in Izmir (Aegean cost), Antalya (Southern Anatoli), and Diyarbakır (southeastern Anatolia) and overpredicted in Sivas (Central Anatolia). In the remaining provinces, there was no consistent over-or under-prediction. In general, over-/ under-prediction was more drastic in the summer months.
Due to the lack of knowledge of internal variables, researchers have developed various ML algorithms to predict evaporation because they provide straightforward solutions to nonlinear multivariate functions. For instance,  Keskin et al. (2004) successfully applied fuzzy logic models as an alternative to classical evaporation estimation formulas for Lake Egirdir in the western region of Türkiye. Lu et al. (2018) proposed the gradient-boosting decision tree model because it was the most accurate and stable model for estimating evaporation in the Poyang Lake watershed. Kumar et al. (2021) recommended the radial function-based SVM model for estimating evaporation under the same climatic conditions and with the same meteorological parameters. Ahi et al. (2022) found that the ANN model had high performance in evaporation estimates with few input parameters. Furthermore, numerous recent studies have been published in the literature to estimate ET 0 using independent variables from specific periods for one or more meteorological stations (Mehdizadeh et al. 2021;Mosre and Suárez 2021;Rashid Niaghi et al. 2021;Kadkhodazadeh et al. 2022;Kim et al. 2022;Zouzou and Citakoglu 2022). In this study, the daily data, including T, RH, W S , and R S , from 213 stations in Türkiye were used to estimate ET 0 . In the Pearson correlation analysis between the input parameters and the output parameter ET 0 , R S had the most significant influence on ET 0 , while W S had the least effect. The argument presented here is in line with previous research (Ge et al. 2022;Bayram and Çıtakoğlu 2023). This study used KNN, SVR, and MGGP models to estimate ET 0 . Based on the results, the SVR model with inputs R S , T, and RH had the best accuracy, while the KNN model had the least accurate results. Bayram and Çıtakoğlu (2023) found that the MGGP model is more robust than the M5Tree and KNN models in modeling monthly ET 0 . However, some authors have found that the performance of the KNN model is better (Yamaç and Todorovic 2020;Al-Mukhtar 2021).
In the second stage of this study, we also investigated the different interpolation techniques, including OK, COK, IDW, and RBF, to predict the daily ET 0 values of Türkiye. The findings of our study showed that OK is the best interpolation method to predict ET 0 in Türkiye. Different studies examining ET 0 have been conducted for spatial distribution using different methods. Dalezios et al. (2002) used geostatistical techniques to study ET 0 and identified spatial variations by kriging estimates over very complex and large Greek fields. Jadhav et al. (2017) used simple kriging (SK), OK, and IDW to determine the spatial distribution of monthly ET 0 values for ten districts in western Maharashtra. The SK interpolation technique was found to be suitable to map ET 0 for all months except February. Geleta et al. (2019) used the Penman-Monteith method to calculate ET 0 from 1991 to 2015 in the Horro Guduru Wollega zone from meteorological data and utilized the OK interpolation technique to estimate temporal and spatial variability. Okechukwu (2020) employed IDW and kriging interpolation techniques for monthly, annual, and seasonal ET 0 and obtained a positive correlation of prediction results by IDW (R: 0.83) and kriging (R: 0.50) for ET 0 .

Conclusions
This study consists of two stages. In the first stage, a daily ET 0 forecast was conducted with different meteorological inputs. In this scenario, the capabilities and accuracy of three ML techniques, namely KNN, SVR, and MGGP, were compared. In the second stage of the study, different interpolation techniques, such as OK, COK, IDW, and RBF, were used to predict the daily ET 0 values of Türkiye. Then, spatial maps of ET 0 were generated using the interpolation technique that provided the best estimation results. The conclusions reported below were drawn from the results. R S , T, and RH proved to be the most influential input combination for ET 0 estimation. Based on statistical and graphical criteria, SVR was the most successful model, followed by MGGP and KNN. The OK interpolation technique outperformed the other methods (COK, IDW, and RBF) in predicting ET 0 in Türkiye.
The study also has some limitations: (1) Data from 213 meteorological stations covering a 52-year period (between 1962 and 2014) in Türkiye were used. (2) Four input parameters were considered for the prediction of ET 0 . (3) Three ML methods (KNN, SVR, and MGGP) and four interpolation techniques (OK, COK, IDW, and RBF) were investigated. In future studies, Deep Learning and hybrid techniques need to be explored to improve the accuracy of ET 0 prediction further.
Author contributions DY, EK, BC, and HS were responsible for conceptualization and methodology; EK and BC contributed to software; EK and HS carried out the formal analysis; DY and BC curated the data; EK and DY took part in writing-original draft preparation and visualization; and BC and HS participated in writing-review and editing and supervision. All authors have read and agreed to the published version of the manuscript.
Funding This research received no external funding.

Conflict of interest
The authors declare that they have no competing interests.
Ethical approval Not applicable.

Informed consent Not applicable.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.