Estimation of Reference Evapotranspiration Using Spatial and Temporal Machine Learning Approaches

Abstract

Evapotranspiration (ET) is widely employed to measure amounts of total water loss between land and atmosphere due to its major contribution to water balance on both regional and global scales. Considering challenges to quantifying nonlinear ET processes, machine learning (ML) techniques have been increasingly utilized to estimate ET due to their powerful advantage of capturing complex nonlinear structures and characteristics. However, limited studies have been conducted in subhumid climates to simulate local and spatial ET_o using common ML methods. The current study aims to present a methodology that exempts local data in ET_o simulation. The present study, therefore, seeks to estimate and compare reference ET (ET_o) using four common ML methods with local and spatial approaches based on continuous 17-year daily climate data from six weather stations across the Red River Valley with subhumid climate. The four ML models have included Gene Expression Programming (GEP), Support Vector Machine (SVM), Multiple Linear Regression (LR), and Random Forest (RF) with three input combinations of maximum and minimum air temperature-based (Tmax, Tmin), mass transfer-based (Tmax, Tmin, U: wind speed), and radiation-based (Rs: solar radiation, Tmax, Tmin) measurements. The estimates yielded by the four ML models were compared against each other by considering spatial and local approaches and four statistical indicators; namely, the root means square error (RMSE), the mean absolute error (MAE), correlation coefficient (r²), and scatter index (SI), which were used to assess the ML model’s performance. The comparison between combinations showed the lowest SI and RMSE values for the RF model with the radiation-based combination. Furthermore, the RF model showed the best performance for all combinations among the four defined models either spatially or locally. In general, the LR, GEP, and SVM models were improved when a local approach was used. The results showed the best performance for the radiation-based combination and the RF model with higher accuracy for all stations either locally or spatially, and the spatial SVM and GEP illustrated the lowest performance among the models and approaches.

1. Introduction

Evapotranspiration (ET) is a combination of two separate processes that transfer huge volumes of water and energy from the soil (evaporation) and vegetation (transpiration) to the atmosphere [1,2]. ET is the second greatest component of hydrological cycle and a major component of the Earth’s surface energy balance. ET closely relates to plant growth, droughts, gas efflux, and production. Since ET plays a crucial role in watershed and agricultural water management, accurate spatial and local estimation of crop water requirements (ET_a) at the scale of human influence is a critical need for a wide range of applications [3]. Quantifying ET_a from agricultural lands is vital to the management of water resources. Measurement methods of ET_a are available through water vapor transfer methods (e.g., Bowen ratio, eddy covariance) [4,5,6], water budget measurements (e.g., soil water balance, weighting lysimeters) [7,8,9], or remote sensing techniques [10,11,12]. However, the application of field scale methods is limited due to cost, complexity, and maintenance. Estimating ET_a using remote sensing models has developed in recent years, but cloud existence in many areas during the satellite passing dates limits the imagery’s usefulness. Therefore, due to difficulties in acquiring direct ET_a measurements, ET_a can be estimated using by multiplying the calculated reference ET (ET_o), using different calculation methods and meteorological data, with the crop coefficient (K_c). However, the required meteorological data for ET_o calculation methods are not readily available for many locations, and methods with fewer variables must be considered when basic meteorological data are available [13]. However, the simplified and basic models are suited to estimating ET_o on a weekly or monthly basis instead of a daily basis [14].

The ET_o calculation is a complex process due to a large number of associated meteorological variables, and it is hard to develop an accurate empirical model to overcome all the complexities of the process [15]. Over the last few decades, machine learning (ML) techniques have attracted the interest of streams of researchers around the globe to overcome the ET_a estimation complexity. These methods are evolutionary computation techniques that can achieve the best relationship in a system with a data driving tool. Due to the capability of the ML methods in tackling the nonlinear relationship between dependent and independent variables [16], numerous ML techniques have been applied and proposed to predict ET_o for agricultural purposes including genetic programming (GP) [17,18]; kernel-based algorithms, e.g., support vector machine (SVM) [19,20]; artificial neural network [21,22,23]; wavelet neural network [14,24]; random forest (RF) [24,25]; and multiple linear regression (LR) [26,27].

Several authors have applied ML methods to detect ET_o values with minimum variation from the observed values [16,21,28]. Among these models, Gene Expression Programming (GEP) and SVM illustrated better estimation than other models [29,30]. Gene Expression Programming (GEP) is comparable to GP and both involve different sizes and shapes of computer programs encoded in linear chromosomes of fixed length [16]. The SVM method is a regression procedure that has been used successfully in the hydrology context [30,31,32] and agro-hydrology for ET_o modeling [19,33]. GEP has several advantages compared to GP such as generating valid structures, its multigenic nature, and the ability to surpasses the old GP system. Shiri et al. [17] evaluated GEP to estimate daily evaporation through spatial and local data scanning Kişi and Çimen [34] studied the potential of the SVM model in ET_o prediction and observed that the SVM model could be useful for ET_o estimation and hydrological modeling studies. Shiri and Kişi [35] compared the GEP and Adaptive neuro-fuzzy inference system (ANFIS) techniques for predicting short and long-term river flows. They also used a similar comparison to predict groundwater table fluctuations.

The LR model is one of the commonly used ML methods in hydrology [28,36]. It has been used to cover the study of the relationship between two or more hydrological variables and investigate the dependence and relationship between the successive values of hydrologic data. Some researchers have used the LR method to estimate the ET_o rate for different regions due to the simplicity of the method compared to other numerical methods [37,38]. Tabari et al. [27] compared the LR and SVM models’ performance for the semi-arid region and found that the SVM model was superior to empirical and LR models.

Due to the high computational costs and complexity of the ML models, the tree-based ensemble models attract people by their simplicity and estimation power. The RF model as an ML tree-based model can produce a great result compared to the other ML models [25,39]. This model is known for its simplicity and the ability to perform both classification and regression tasks. The RF method is also widely used to predict the ET rate of different climate regions [40]. Feng et al. [26] applied the RF model for ET_o estimation daily for southwest China and indicated that the RF model performed slightly better than the general regression neural network model. Shiri [41] evaluated the coupled RF with wavelet algorithm to estimate ET_o for the southern part of Iran and obtained results stating that the coupled RF model showed a great improvement compared to the conventional RF and empirical models. To our knowledge, this model has not been applied in the Northern US for ET_o studies.

According to the literature, GEP and SVM models have been frequently applied across the world in various climate conditions for ET_o estimation, while the LR and RF models’ applications were minimal. Besides, these models have not been compared with commonly used SVM and GEP models for subhumid climate conditions. Since the limited studies have been conducted to evaluate ML models for the Northern part of the US (which experiences a high variability of weather conditions and a huge amount of agricultural production), the objectives of this study were to (1) investigate the effect of different input combinations of meteorological data on the accuracy of daily ET_o estimation in subhumid climate using the GEP, SVM, LR, and RF methods, (2) compare the spatial and local prediction capabilities of the different ML models in ET_o estimation, and (3) evaluate the performance of the models based on the various study years and meteorological stations.

2. Materials and Methods

2.1. Study Area Climate and Reference Evapotranspiration (ET_o)

The weather data for the current study were obtained from the North Dakota (ND) Agricultural Weather Stations [42] located at Prosper (ND), Galesburg (ND), Leonard (ND), Sabin (MN: Minnesota), Perley (MN), and Fargo (ND) for 17 study years (January 2003–December 2016). Geographical location of the study stations is shown in Figure 1.

To calculate the daily reference evapotranspiration (ET_o) for each study station, the ASCE standardized reference evapotranspiration equation (ASCE-EWRI) was used for the alfalfa reference crop [43]. This equation provides a standardized calculation of ET_o demand that can be used in developing K_c and comparing it with other methods. Equation (1) presents the form of the standardized ET_o equation for all hourly and daily calculation time steps.

$$\mathrm{E}{\mathrm{T}}_{\mathrm{o}}=\frac{0.408 \Delta \left(R_n-G\right)+\gamma \frac{{\mathrm{C}}_{\mathrm{n}}}{T+273}{u}_2\left(e_s-e_a\right)}{\Delta +\gamma \left(1+C_d{u}_2\right)}$$

(1)

where ET_o is reference evapotranspiration rate (mm d⁻¹), R_n is net solar radiation (MJ m⁻² d⁻¹), G is soil heat flux (MJ m⁻² d⁻¹), $\gamma$ is psychrometric constant (KPa °C⁻¹), T is the mean daily air temperature, U₂ is the average wind speed at 2 m height (m s⁻¹), es is saturation vapor pressure (KPa), e_a is actual vapor pressure (KPa),$\Delta$ is the slope of the saturation vapor pressure–temperature relationship (KPa °C⁻¹), C_n and C_d are coefficients which are related to the crop and time step. The value for the constants C_n and C_d are 1600 and 0.38 for the alfalfa reference crop.

Table 1. Statistical summary of the weather data for the study locations for the study period.

Station	Parameter	Unit	Xmax	Xmin	Xmean	SX	CV	CSX
Prosper, ND	Tmax	°C	37.9	24.3	11.3	14.3	1.27	−0.37
	Tmin	°C	−29.8	−38.1	−0.8	13.0	−16.73	−0.28
	WS	m s−1	14.2	0.9	4.2	1.8	0.43	0.55
	Rh	%	100	13.8	68.6	15.6	0.23	−0.14
	RS	MJ m−2	31.1	0.3	13.2	7.9	0.60	0.51
	ETo	mm	11.4	0	2.4	2.03	0.84	0.92
Galesburg, ND	Tmax	°C	36.8	23.6	10.9	14.2	1.30	−0.33
	Tmin	°C	−28.9	−37.3	−1.0	12.7	−12.41	−0.28
	WS	m s−1	12.8	0.7	3.9	1.6	0.41	0.45
	Rh	%	100	18.8	68.1	15.2	0.22	−0.09
	RS	MJ m−2	30.7	0.2	12.8	7.9	0.61	0.51
	ETo	mm	10.6	0	2.3	1.97	0.85	1.03
Leonard, ND	Tmax	°C	38.3	23.6	11.5	14.2	1.23	−0.39
	Tmin	°C	−28.6	−37.7	−0.6	12.9	−21.05	−0.28
	WS	m s−1	13.2	0.9	4.2	1.7	0.42	0.50
	Rh	%	100	17.85	67.40	15.3	0.23	−0.02
	RS	MJ m−2	31.6	8.1	13.6	8.1	0.60	0.51
	ETo	mm	10.6	0	2.5	2.09	0.85	0.77
Sabin, MN	Tmax	°C	37.8	24.3	11.2	14.1	1.26	−0.33
	Tmin	°C	−30.2	−38.5	−0.2	13.0	−73.34	−0.24
	WS	m s−1	12.7	0.5	4.0	1.7	0.42	0.46
	Rh	%	100	18.70	68.80	14.9	0.22	−0.08
	RS	MJ m−2	31.6	0.4	13.0	7.9	0.61	0.51
	ETo	mm	10.1	0	2.4	2.02	0.86	0.85
Perley, MN	Tmax	°C	37.3	24.1	10.9	14.3	1.31	−0.36
	Tmin	°C	−30.5	−40.7	−0.7	13.1	−18.07	−0.30
	WS	m s−1	11.8	0.8	4.1	1.7	0.41	0.48
	Rh	%	100	17.22	69.10	14.9	0.22	−0.08
	RS	MJ m−2	31.3	0.4	12.8	7.9	0.61	0.51
	ETo	mm	10.9	0	2.3	2.02	0.84	1.12
Fargo, ND	Tmax	°C	39.6	25.6	11.4	14.2	1.24	−0.36
	Tmin	°C	−29.5	−36.8	0.6	13.0	21.89	−0.23
	WS	m s−1	11.3	0.8	3.8	1.5	0.39	0.40
	Rh	%	100	15.55	66.19	14.9	0.23	−0.05
	RS	MJ m−2	31.0	0.1	12.8	7.9	0.61	0.52
	ETo	mm	10.5	0	2.5	2.07	0.84	0.92

summarized weather parameters of the study locations for the study period.

2.2. Models Structure and Application

To process the GEP and SVM algorithm, we used the GeneXpro program in Matlab, and to process the LR and RF models, a scikit-learn module embedded in the Python 3.2 programming language was used. GEP is an extension of GP [44] developed by Ferreira [45] that creates a computer program to investigate a relationship between input and output variables. GEP was developed to find a better solution for solving a particular problem relating to the understudied phenomena [46].

The application of GEP requires several steps. GEP is, like GAs and GP, a genetic algorithm, as it uses populations of individuals, selects them according to fitness, and introduces genetic variation using one or more genetic operators. The procedure to model daily evapotranspiration (as a dependent variable) by using weather variables (as independent variables) is as follows: 1. Selection of fitness function; 2. Choosing the set of terminals T and the set of functions F to create the chromosomes; 3. Choosing the chromosomal architecture; 4., Choosing the linking function; 5. Choosing the genetic operators.

The fitness function must be determined in the first step with a random generation of chromosomes of a certain program (initial population) and evaluated against a set of fitness cases [47]. Using weather station data as input variables (terminals) to model daily ET_o involves the next general step. The selection of fitness functions (i.e., absolute error, relative error, and correlation coefficient) depends on the experience and intuition of the user. The GEP model in the current study was developed based on the recommended functions by Shiri et al. [17]. In the third step, the chromosomal architecture can be defined by having the weather variables as terminal and function set as chromosomes. The fourth step was to select a linking function that relates genes to each other in addition to linking the parse trees [17]. Finally, genetic operators’ corresponding rates were chosen.

Table 2. Parameters used per run of gene expression programming (GEP).

Number of Chromosomes	30	One-Point Recombination Rate	0.3
Head of the size	8	Two-point recombination rate	0.3
Number of genes	3	Gene recombination rate	0.1
Linking function	Addition	Gene transposition rate	0.1
Fitness function error type	RMSE	Insertion sequence transposition rate	0.1
Mutation rate	0.044	Root insertion sequence transposition	0.1
Inversion rate	0.1	Penalizing tool	parsimony pressure

summarizes the commonly used parameters for each run.

The SVM was developed by Cortes and Vapnik [48] and is known as the classification and regression method [34] to solve problems by applying a flexible representation of the class boundaries and implementing an automatic complexity control to reduce overfitting. In SVM, the dependency of the dependent variable to a set of independent variables is evaluated. In regression estimation with Support Vector Regression (SVR), which is used to define SVMs in the literature, a functional dependency $f\left(x\right)$ between a set of sampled input points $X=\left\{x_1,x_2,x_3,\dots ,x_l\right\}$ (here, input sampled refer to meteorological variables) taken from $R_n$ (input vector of n dimension) and target values $Y=\left\{y_1,y_2,y_3,\dots ,y_l\right\}$ (ET_o as target values) with$y_i\in R_n$. More detail on SVM can be found in Vapnik [49].

The LR is a statistical method used to describe a quantitative relationship between a dependent variable and one or more independent variables [27,50]. In LR, the function is a linear equation and is expressed as:

Y_i = b_o + b_{1 × 1} + …. + b_kx_k

(2)

where b_o–b_k are the fitting constants, y_i is the dependent variable, and x₁–x_k are the independent variables for this system.

The RF method combines a group of decision trees for either classification or regression purposes. Although each decision tree may not be capable of learning well, the combination of the decision trees results in a strong learner. Each decision tree predicts the outcome individually, and RF votes among the outcomes for classification or averages the outcomes for regression. Each decision tree is trained on a different subset of samples by a bagging extension of the RF model to reduce the risk of overfitting. Moreover, a different subset of input variables can be used in each tree to make it more useful in prediction for datasets with higher dimensions [51]. For this study, a small subset of data was used to find a good combination of parameters for the RF model. As a result, there were several trees in the forest and the minimum number of samples required for the leaf nodes were 50 and 35, respectively. The mean square error criteria are used as a procedure for estimation.

The calculated daily ET_o was used to feed the GEP, SVM, LR, and RF models. Three treatments including temperature, radiation, and mass transfer-based combinations were used as input to feed the models, and each model of the combinations was assessed for spatial and local approaches. Different statistical analysis was performed to evaluate the accuracy and performance of the different combinations and approaches for each studied station. The combinations were as follows:

(i)

T_min, T_max: temperature based (GEP1, SVM1, LR1, RF1)

(ii)

T_min, T_max, R_s: radiation-based (GEP2, SVM2, LR2, RF2)

(iii)

T_min, T_max, W: mass transfer based (GEP3, SVM3, LR3, RF3).

2.3. K-Fold Cross-Validation

Splitting the data into the sets of data for testing and training is a usual procedure for assessing the ML techniques. Using 10–30% of the complete dataset as a single test set is a common method for GEP evaluation. Therefore, the K-fold cross-validation technique was used to increase the evaluation performance and set of data for either training or testing purposes. Using K-fold cross-validation, the dataset was divided into K subsets, and the training process was repeated K times leaving each time a distinct set of patterns for testing until a complete testing scan for the dataset was fulfilled. Computation cost defines the minimum assembly size of the test set. Here, the minimum test size was fixed as one year for local modeling and one station for spatial modeling. Consequently, at a local scale, one year was held out each time for testing while the models were trained using the remaining 16 years; hence, a total of 612 models (17 years × 6 stations × 3 input combinations × 2 models) were established for the local k-fold testing. At the spatial scale, one station was considered as a test block each time and the models were trained using the patterns from the remaining stations; hence, a total number of 36 models (6 stations × 3 input combinations × 2 models) were constructed. The local and spatial approaches were noted with T and S in the figures.

2.4. Evaluation Criteria

To investigate the performance of the models for each combination and approaches, four statistical indicators were used, namely, the root mean square error (RMSE), the mean absolute error (MAE), correlation coefficient (r²), and scatter index (SI), defined as follows:

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(\mathrm{E}{\mathrm{T}}_{\mathrm{e}}-\mathrm{E}{\mathrm{T}}_{\mathrm{o}}\right)}^2}$$

(3)

$$\mathrm{MAE}=\frac{{{\displaystyle \sum }}_{i=1}^n\left|\mathrm{E}{\mathrm{T}}_{\mathrm{e}}-\mathrm{E}{\mathrm{T}}_{\mathrm{o}}\right|}{n}$$

(4)

$$r^2={(\frac{{{\displaystyle \sum }}_{i=1}^n\left(\mathrm{E}{\mathrm{T}}_{\mathrm{o}}-\overline{\mathrm{E}{\mathrm{T}}_{\mathrm{o}}}\right)\left(\mathrm{E}{\mathrm{T}}_{\mathrm{e}}-\overline{\mathrm{E}{\mathrm{T}}_{\mathrm{e}}}\right)}{\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left(\mathrm{E}{\mathrm{T}}_{\mathrm{o}}-\overline{\mathrm{E}{\mathrm{T}}_{\mathrm{o}}}\right)}^2{{\displaystyle \sum }}_{i=1}^n{(\mathrm{E}{\mathrm{T}}_{\mathrm{e}}-\overline{\mathrm{E}{\mathrm{T}}_{\mathrm{e}}})}^2}})}^2$$

(5)

$$\mathrm{SI}=\frac{\mathrm{RMSE}}{\overline{{\mathrm{E}}_{\mathrm{o}}}}$$

(6)

where ${\mathrm{ET}}_{\mathrm{e}}$ and ${\mathrm{ET}}_{\mathrm{o}}$ are simulated and calculated reference evapotranspiration at the i-th time step, respectively, n is a number of time steps, $\overline{{\mathrm{ET}}_{\mathrm{e}}}$ and $\overline{{\mathrm{ET}}_{\mathrm{o}}}$ are mean values of simulated and calculated ET_o, respectively.

The RMSE describes the average magnitude of errors and can take on values from 0 to ∞ indicate perfect and worst fit, respectively, and the MAE scores the error magnitudes without any specific weight to larger/smaller errors. Therefore, the lower value of the RMSE and MAE is desirable. The r² values around 1 indicate a perfect linear relationship between estimated and calculated values, where the closer a value is to zero, the more it demonstrates the poor relationship between simulation and calculation. Finally, SI is a dimensionless index of RMSE that gives a good insight to compare the performance of different models.

3. Results and Discussions

The local and spatial analysis of four models for six studied stations is shown in

Table 3. Global average performance indicators of the gene expression programming (GEP), support vector machine (SVM), multiple linear regression (LR), and random forest (RF) methods for three input combinations of local (T) and spatial (S) approaches.

Evaluation Criteria		Input Combination
		1 (Temperature-Based)				2 (Radiation-Based)				3 (Mass-transfer-based)
	Approach	GEP	SVM	LR	RF	GEP	SVM	LR	RF	GEP	SVM	LR	RF
R2	T S	0.75 0.78	0.80 0.75	0.77 0.77	0.85 0.84	0.85 0.87	0.91 0.85	0.88 0.88	0.92 0.93	0.77 0.77	0.86 0.77	0.78 0.78	0.86 0.88
RMSE (mm/day)	T S	0.90 1.07	0.97 1.13	0.97 0.98	0.82 0.80	0.71 0.76	0.72 0.77	0.68 0.69	0.57 0.55	0.72 0.91	0.73 0.93	0.94 0.95	0.73 0.69
MAE (mm/day)	T S	0.64 0.84	0.71 0.82	0.76 0.77	0.58 0.57	0.50 0.57	0.54 0.61	0.51 0.52	0.38 0.36	0.64 0.69	0.62 0.67	0.75 0.76	0.53 0.49
SI	T S	0.38 0.44	0.41 0.46	0.40 0.40	0.34 0.33	0.29 0.32	0.30 0.33	0.28 0.29	0.24 0.23	0.35 0.38	0.33 0.36	0.39 0.39	0.30 0.28

. According to the three combinations’ performance, the radiation-based method illustrated the highest accuracy for either local or spatial approaches compared to the other combinations. The mass-transfer-based combination was the next most accurate combination. The results showed that the local trained models surpassed the spatially trained models because of using the same patterns for both training and testing models. However, the spatial models gave comparable results compared to the local model in some cases, especially for radiation-based combinations. Differences in temperature among the stations have dramatically affected the performance of both the temperature-based and the mass transfer-based models. In all cases, the minimum differences between the performance accuracy of the local and spatial models belonged to the LR model. This can be inferred to the mathematical structure of this technique, where only linear relationships can be supposed between the input and target parameters with a lower degree of flexibility compared to heuristic data driven models.

Among four models with three input combinations, the models relying on radiation, mass-transfer, and temperature-based combinations showed the lowest RMSE and MAE, respectively (Table 3). Comparing the GEP, SVM, LR, and RF models, the RF model illustrated the lowest rate of RMSE and MAE with the best performance for radiation-based approaches. However, the RF model improved 4.37, 5.76, and 1.49 percent from local to spatial approaches for temperature, radiation, and mass-transfer-based combinations, respectively, which was in contrast with the improvement’s direction in the other models. Considering the models based on radiation combination, the spatial RF model exhibited the highest linear relationship (r² = 0.927) between calculated and estimated ET_o in comparison with the other models. The local RF method was the next accurate approach to estimate ET_o based on radiation-based data. This observation illustrated the ability of the RF algorithms to estimate ET_o using data from local stations for training. Furthermore, the LR model had significant improvement for RMSE and MAE from spatial to local approaches for all three types of input combinations. For the LR model, the r² value was not changed for radiation-based and temperature-based combination from spatial to local approaches and the change was 0.13 percent for the mass-transfer-based model. Therefore, the LR model illustrated almost similar results for both spatial and local approaches among all models.

The GEP and SVM models illustrated the great improvement rate for all three input combinations from spatial to local condition with the highest improvement of 21 percent for mass-transfer-based combination. Specifically, the GEP model showed an improvement from spatial to local approach, however, the percentage of improvement was 2.3, 3.9, and 0.13 for radiation, temperature, and mass-transfer-based combinations, respectively. In term of obtained improvement for RMSE and MAE from spatial to local approaches, both of the GEP and SVM models gained similar results. The correlation coefficient of the SVM model decreased from spatial to local approaches for about 6.3, 6.5, and 10.9 percent for radiation, temperature, and mass-transfer-based combinations, respectively. By using local radiation data for training the models, the SI indicator for the GEP and SVM models showed an improvement of 8.2 and 10 percent from spatial to local approach, respectively. This improvement was about 6.6 and 8.7 percent for mass-transfer-based and 15 and 10.9 percent for temperature-based combinations, respectively.

Statistical analysis revealed the similar performance of the local GEP and SVM models. For RMSE and MAE statistical variables, GEP and SVM models showed a greater improvement in performance for mass-transfer and temperature-based combinations, respectively. By considering correlation coefficient values, it can be concluded that the improvement in accuracy of either GEP or SVM approaches was not significant and all illustrated the ability to estimate with acceptable accuracy. Therefore, if temperature data are not available at a particular station, but they are for other stations, the GEP and SVM approaches can be useful to estimate ET_o. However, due to the higher mapping ability of the GEP models, using either local or spatial GEP are preferable.

The models relying on the mass-transfer combination had slightly higher accuracy than the temperature-based approach, but lower accuracy compared to the radiation-based combination. All of the local and spatial GEP and SVM methods illustrated lower improvement compared to that for the temperature-based approach. This showed that wind speed can have a significant effect on the accurate estimation of spatial and local ET_o. Due to the flat topography of the study area and being faced with lots of high-speed winds during the growing season and almost all other seasons, including the wind as a parameter to build the model architecture and estimating the ET_o can increase the accuracy of the approach.

Overall, the RF and the LR models illustrated the best performance among the four models, and comparing the GEP and SVM models, the GEP model showed better performance than the SVM model for all three input combinations.

A breakdown of the models’ performance accuracy at each station is shown in Figure 2, Figure 3 and Figure 4 for all of the three input combinations, respectively. In the case of the temperature-based combination (Figure 2), the local GEP and SVM models (shown as TGEP and TSVM) gave more accurate results than the spatial (shown as SGEP and SSVM) models. For the LR and RF models, the difference in accuracy between local (TLR and TRF) and spatial approaches (SLR and SRF) was not significant and both showed better performance than GEP and SVM models since they relied on the patterns of the same location used for training and testing the models. According to Table 2, station 6 (Fargo) had the highest, and station 2 (Galesburg) had the lowest range of recorded temperature among the study stations. This range may be caused to have the lowest performance for station 2; however, it was difficult to evaluate the model’s performance in the climate context of each station due to the few number of study stations. The RF model showed the best performance with higher accuracy for all stations either locally or spatially, and the spatial SVM and GEP illustrated the lowest performance among other models and approaches.

The RF and LR methods showed the lowest range of SI compared to the spatial and local GEP and SVM methods. For temperature-based combinations, the spatial and local LR approaches had minimum SI ranges of 0.018 and 0.020, respectively, and the spatial SVM and GEP methods illustrated the highest SI ranges of 0.113 and 0.119, respectively. The spatial RF approaches with an average of 0.333, and spatial SVM, with an average of 0.457, showed the lowest and highest SI rate, respectively. Therefore, spatial RF approaches may be the most practical way to estimate the missing meteorological data of the study stations.

Figure 3 shows the evaluation result of the radiation-based combination for the four models with spatial and local approaches. The amount of received radiation for all study stations was similar. According to the global performance of the defined models, the radiation-based combination gained the best performance among the three input combinations. Besides, the radiation-based combination had the lowest rate of RMSE, MAE, and SI, and the highest rate of r² for each of the study stations. Among the spatial and local scenarios, the local approach had a better performance than the spatial approach. For the radiation-based combination, the spatial RF and local RF models had an accurate estimation of ET_o, respectively. For stations 3 and 4 (Leonard and Sabin) either spatial or local approaches of GEP and SVM models gained lower performance than the other stations. This could be due to the slightly higher magnitude and variations of solar radiation (Table 2) among the other stations during the study period.

Among the study stations comparison, the SI range of the spatial RF was 0.018, which showed the best performance compared to the other applied methods. As obtained from the temperature-based combination, the LR method performed well in the radiation-based combination too, with an SI indicator range of 0.021 and 0.024 for local and spatial LR approaches, respectively. The worst performance was observed for spatial GEP and SVM approaches, with SI indicator of 0.128 and 0.140, respectively. According to the GEP and SVM models, the local GEP performed well compared to other approaches of the SVM and GEP models. The statistical indicators were in agreement with the spatial RF performance in which they showed the lowest rate of RMSE and MAE and the highest value of r². However, comparing the MAE might not be a valid indicator due to taking into account the local order of magnitude of the target variable. The ranking of the SI indicators showed that spatial RF and LR could overcome the lack of meteorological data for the station. On the other hand, the averages of the SI values for all six study stations showed that the spatial RF and local RF had the lowest and the spatial GEP and spatial SVM had the highest rate of SI indicators. Therefore, either spatial or local RF methods could be useful to estimate the missing values for any of the stations.

Figure 4 shows the statistical indices of the mass-transfer-based combination. Similar to the previous combinations, the spatial and local RF gave a more accurate estimation than other methods. On the other hand, the local SVM approach showed better estimation than the spatial SVM and GEP methods for all stations except station 2, which had the lowest range of temperature variation. The fluctuations of the indices among the stations were higher than the radiation-based combination and lower than that for temperature-based combination, which showed mediocre accuracy compared to the other combinations.

By having wind speed and temperature data as the input variable for the mass transfer-based combination, the spatial RF approach gained the lowest SI and highest r² values for ET_o estimation compared to all other methods. The minimum and maximum SI values for mass transfer-based combination were obtained for the spatial RF and spatial SVM approaches, which were 0.011 and 0.120, respectively. According to the performance ranking of the models based on the SI indicator, spatial LR, local LR, and local RF showed better performance after spatial RF with SI values of 0.015, 0.018, and 0.018, respectively. The local SVM, local GEP, and spatial GEP had the SI values of 0.087, 0.10, and 0.119, respectively. The average of SI for all study stations showed that the spatial and local LR had the highest and spatial and local RF had the lowest SI values, respectively. Therefore, by having the lowest range of SI and lowest value of SI for the spatial RF approach, it might be more practical to apply the spatial RF for other stations without training a specific model for each station. Accordingly, no local dataset would be needed to train the local models. This could be helpful to estimate the ET_o for stations with partial or missing meteorological data.

To understand the yearly performance of the applied models based at each of the study stations, the models were assessed per test year. Figure 5 illustrates the SI values obtained from the three input combinations for each study year of the study stations. The SI values of the models fluctuated considerably for almost all stations during the test years.

As shown in Figure 5, the SI values fluctuated considerably within test years for all input combinations and approaches. Among the study stations, Prosper and Sabin stations showed the average maximum range for the SI values. The minimum average of the SI value 0.223 was observed for the RF radiation-based combination for the Fargo station, and the maximum average of the SI was obtained for the SVM temperature-based combination (SVM1) for the Galesburg station. The Galesburg station had the lowest temperature range among the study stations. According to Figure 4, the RF2 (radiation-based combination) model showed the lowest fluctuation compared to the other approaches with a similar trend among the study stations. For all of the study stations, test years, and input combinations, the RF models gave the best performance with the lowest SI values compared to the other models and combinations. The SVM and GEP models showed the worst SI averages for the temperature-based combinations. In this case, the order of the accuracy of the models was similar to that obtained from the station-based analysis. The radiation-based combination gave the most accurate results and estimation in comparison with the temperature or mass-transfer-based combination models.

Comparing the performance of the models relying on each of the combination methods, it can be observed that the performance of the approaches is similar. However, the range of the SI indicator for different approaches was different depending on the test years. For example, 2011, 2012, and 2013 were the dry, normal, and wet years among the study years, respectively. The result of SI per test years showed that 2012 had lower SI than 2011 and 2013 for all of the three input combinations, and the SI values of the various methods and approaches were close together. On the other hand, for the 2013 test year, some of the approaches illustrated a huge jump for the obtained SI from 2012 to 2013 due to the weakness of the model performance.

Due to the variability in the meteorological data and the climate pattern, the variability in performance accuracy at each station was expected. Similar variability in performance of the ML approaches was observed by Shiri et al. [17]. Selection of the training set and testing set plays an important role in model performance. The existence of any abnormal year in the test years in comparison with training datasets causes it to have an inaccurate estimation [19]. By lowering the number of the input values, the validity of the training set for estimation of test years decreases. Because of the lower input values, the variable rate would be low, and this type of input would only be valid for periods with very specific trends. This explanation may clarify the performance of spatial approaches, where the relationship encountered might be site-specific. Other researchers illustrated the site-specific performance of spatial approaches for the different study regions and climate conditions [34,35].

The comparison of the three input combinations showed that the performance of some approaches was similar for some years while the performance of methods for some years were far from each other. For example, the SVM model showed the most improvement from temperature to mass-transfer-based combination, which became like the RF method. However, depending on the station and test year, this similarity becomes even closer. All the test years and stations showed an improvement from temperature to radiation or mass-transfer-based combination except for the Prosper station, which is in agreement with the findings of Adnan et al. [15]. On the other hand, the Prosper station showed the best improvement for the SVM model for radiation-based approaches. Considering that all of the input combinations rely on the temperature and another variable (solar radiation or wind speed), it might be thought that the performance differences could be due to the effect of the second variable in the estimation of the output. Besides, when the performance of the models is similar, the impact of the secondary variable might be less than the primary variable (temperature). However, when the gap between the performance of the SI indicator increases, it proves the crucial impact of the second variable on the model performance for the specific test year or station. A similar conclusion could be drawn for the comparison between the input combinations. If each of the combinations shows a better performance than the other combination method for a specific year and station, the second variable effect should be important and might have a significant influence in the explanation of the output for that test year.

4. Conclusions

This paper aimed to evaluate the different ML methods including GEP, SVM, LR, and RF to estimate ET_o in the Red River Valley. The external approach exempts using local data like spatial approaches in the current paper for simulating ET_o values in a decisive way when local data is not available, reliable, or sufficient. Global comparison of the performance accuracy of the applied models revealed that the RF model was the best for all combinations among the four defined models. Furthermore, the RF model illustrated the best performance for spatial and local conditions for all input combinations. In general, the LR, GEP, and SVM models were improved when a local approach was used, except for the RF model, which was less accurate with a local approach. The radiation-based combination was the most accurate predictor among all models tested. As a result, this combination showed the lowest rate of improvement due to better performance in the first step.

The results showed that due to the flat topography of the study area with high wind speeds during the growing season, including the wind as a parameter to build the model architecture and estimate the ET_o can increase the accuracy of the prediction. Besides, it might be more practical to apply the spatial RF model for stations with missing meteorological data without the need for local training. The recommended application of spatial RF using radiation combination allows for a more reliable estimate of ET_o to fill the missing values for more precision water management purposes. Further research should confirm the current results in other geographical locations and for the various input combination methods.