Optimal Alternative for Quantifying Reference Evapotranspiration in Northern Xinjiang

: Accurate estimation of reference evapotranspiration is a key step in irrigation and water resources planning. The Penman Monteith (FAO56-PM) formula recommended by FAO56-PM is the standard for calculating the reference evapotranspiration. However, the FAO56-PM model is limited in the observation of meteorological variables, so it is necessary to choose an alternative ET 0 model which requires less meteorological data. Based on the daily climate data of eight meteorological stations in northern Xinjiang from 2000 to 2020, seven empirical models (Hargreaves, Berti, Dorji, Dalton, Meyer, WMO, Albrecht) and four optimization algorithms (RF model, LS-SVR model, Bi-LSTM model and GA-BP model) combined with seven different parameters were evaluated comprehensively. The results show that the accurate of the empirical model based on temperature is obviously better than the empirical model based on air mass transport. The annual and multi-year alternative ET 0 models of different input parameter combinations are: LS-SVR1, RF2, LS-SVR3, LS-SVR4, GA-BP5, LS-SVR6, GA-BP7. It can be used as a substitute for the reference evapotranspiration model without relevant meteorological data. Only the LS-SVR6 model and GA-BP7 model are recommended as the best alternative models for northern Xinjiang reference evapotranspiration at daily, monthly and seasonal scales.


Introduction
The quantification and accurate estimation of evapotranspiration are of great significance to the formulation of farmland irrigation systems, the study of hydrology and water balance, and the planning of water resources. The Penman Monteith (FAO56-PM) method recommended by the FAO [1] is the standard method for estimating reference evapotranspiration (ET 0 ) thus far [2,3]. This method does not need to be initially calibrated locally and has global applicability [4][5][6]. However, this method has many input parameters (air temperature, humidity, solar radiation and wind speed), so the acquisition of many meteorological parameters becomes the only limitation of its application [7]. In this case, relatively simple empirical models are usually used to estimate ET 0 , and the selection of the optimal empirical model is of great significance to water resource planning and management [8,9]. Existing empirical evapotranspiration models mainly include temperature-based, radiation-based, air-mass-transport-based and combination models. The results show that the combined models are better, followed by radiation-based, temperature-based and aerodynamic models [10]. However, the combined model is highly dependent on the input of meteorological variables, which requires the sufficient input of meteorological variables, such as net radiation, soil heat flux, air temperature, wind speed and relative humidity. Moreover, some scholars also pointed out that, compared with the climate in humid and semi-humid regions, ET 0 in arid and semi-arid regions is mainly affected by aerodynamics and water vapor pressure deficit rather than radiation availability [11,12]. Therefore, the ET 0 model based on temperature and air mass transport was selected in this study to verify its applicability in northern Xinjiang.
In recent years, more scholars have applied artificial intelligence (AI) to the prediction and estimation of ET 0 and published a number of papers. For example, Vahid Nourani et al. [13] used a feed-forward neural network (FFNN), adaptive neuro-fuzzy reasoning system (ANFIS), support vector machine regression (SVR) and other algorithms to simulate the ET 0 several climate zones of Turkey, Cyprus, Iraq, Iran and Libya. Ahmed Elbeltagi et al. [14] estimated long-term ET 0 in Egypt through a deep neural network (DNN). Babak Mohammadi and Saeid Mehdizadeh [15] used support vector regression (SVR) and random forest (RF) to simulate daily ET 0 from Isfahan, Urmiya and Yazd, Iran. Tongren Xu et al. [16] evaluated the adaptability of machine learning, remote sensing and land surface ET 0 products in the United States. Saman Maroufpoor et al. [17] simulated ET 0 in Iran by using artificial neural network optimization (ANN-GWO). Francesco Granata [18] used the M5P regression tree, bagging, random forest (RF) and support vector machine regression (SVR) to simulate ET 0 in central Florida. Juan Yin et al. [19] simulated ET 0 in the central Ningxia region of China by using the mixed bidirectional long-and short-term memory model (BI-LSTM). Jingran Liu et al. [20] predicted the actual evapotranspiration of green pepper by using an extended neural network (MEA-ENN, GA-ENN) optimized by the mind evolutionary algorithm and genetic algorithm. The results show that artificial intelligence can accurately predict daily ET 0 and is a powerful tool for modelling ET 0 using incomplete meteorological parameters.
In summary, this study selected eight widely used empirical models for performance evaluation, including four temperature-based models: Hargreaves model [21], Berti model [22] and Dorji model [23] and four mass transfer-based models: Dalton Model [24], Meyer Model [25], WMO Model [26] and Albrecht Model [27]. Four optimization algorithms (RF model, LS-SVR model, BI-LSTM model and GA-BP model) were combined with 28 optimization algorithm models to make a comprehensive ranking. On this basis, this study proposes the following hypotheses. First, the empirical model and the ET 0 model based on the optimization algorithm are significantly affected by time and space in northern Xinjiang. Second, simple linear regression and the global performance index GPI can effectively verify the performance effect of the ET 0 model in northern Xinjiang. The research objectives of this study are: (1) to use linear regression to verify seven empirical models to determine the best substitute for the FAO56-PM model in northern Xinjiang; (2) to use the global performance index GPI to rank the 28 models based on the optimization algorithm and to determine the best algorithm model under seven different parameter input combinations; and (3) to discuss the influence of time scale on the model and to recommend the most appropriate reference evapotranspiration estimation model for northern Xinjiang at different time scales. The innovation of this study lies in the first application of eight empirical models and four models based on an optimization algorithm in northern Xinjiang and the recommendation of corresponding models in the absence of meteorological data to fill the gap in knowledge regarding models suitable for northern Xinjiang under different meteorological parameters.

Study Area
Xinjiang is located in the hinterland of Eurasia and the northwest border of China. The Tianshan Mountains traverse the whole territory, dividing the whole territory into northern and southern Xinjiang. The part north of the Tianshan Mountains is called northern Xinjiang. In northern Xinjiang lies the Altai Mountains, and between the Altai Mountains and Tianshan Mountains lies the semi-closed Junggar Basin (Figure 1). Northern Xinjiang has a temperate continental arid climate, and the annual average temperature ranges from 1.96 • C to 9.06 • C. The annual precipitation ranges from 137.75 mm to 253.96 mm, and the annual evaporation ranges from 909.92 mm to 1328.76 mm (Table 1). In the eastern and western regions, the wind speed is high, and the number of windy days is high. The mountainous terrain is undulating, and the wind is blocked. The number of gale days in the plains area is greater than that in the middle and low mountainous areas, and the number of gale days in the Alashankou area in the western Junggar Basin is the highest. In short, due to the dry climate, less precipitation, sparse vegetation, loose soil quality and high winds, the annual ET 0 in northern Xinjiang is relatively large, and it is vulnerable to drought stress.
Xinjiang has a temperate continental arid climate, and the annual average temperature ranges from 1.96 °C to 9.06 °C. The annual precipitation ranges from 137.75 mm to 253.96 mm, and the annual evaporation ranges from 909.92 mm to 1328.76 mm (Table 1). In the eastern and western regions, the wind speed is high, and the number of windy days is high. The mountainous terrain is undulating, and the wind is blocked. The number of gale days in the plains area is greater than that in the middle and low mountainous areas, and the number of gale days in the Alashankou area in the western Junggar Basin is the highest. In short, due to the dry climate, less precipitation, sparse vegetation, loose soil quality and high winds, the annual ET0 in northern Xinjiang is relatively large, and it is vulnerable to drought stress.

Data Sources
Eight weather stations in northern Xinjiang were selected in this study (Figure 1), and daily meteorological data were collected from 2000 to 2020. The observation data include temperature, humidity, wind speed, rainfall, saturated water vapour pressure deficit and extraterrestrial radiation ( Table 1). The climatic data of historical days are obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn (access on 24 November 2021)). Among them, the Wujiaqu, Wusu, Shawan and Qinghe datasets span from 2000 to 2017, and the datasets of Hutubi, Fuyun, Hebukesel and Karamay span from 2000 to 2020.

Data Sources
Eight weather stations in northern Xinjiang were selected in this study (Figure 1), and daily meteorological data were collected from 2000 to 2020. The observation data include temperature, humidity, wind speed, rainfall, saturated water vapour pressure deficit and extraterrestrial radiation ( Table 1). The climatic data of historical days are obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn (access on 24 November 2021)). Among them, the Wujiaqu, Wusu, Shawan and Qinghe datasets span from 2000 to 2017, and the datasets of Hutubi, Fuyun, Hebukesel and Karamay span from 2000 to 2020.

FAO56 Penman-Monteith Model
The FAO56 P-M approach proposed by Allen et al. [1] served here as a reference method, it is expressed by Equation (1): where ET 0 is the reference evapotranspiration (mm day −1 ), Rn is net radiation (MJm −2 day −1 ), G is soil heat flux (MJm −2 day −1 ), γ is the psychometric constant (kPa • C −1 ), T is the mean air temperature ( • C), U 2 is the wind speed at 2 m height (ms −1 ), e s is the saturation vapour pressure (kPa), e a is the actual vapour pressure (kPa), and ∆ is the slope of the vapour pressure curve (kPa • C −1 ).

Empirical Models
Temperature-based models performed well in all subregions of the world. Due to the special climatic conditions in desert areas, it is necessary to select an empirical model that is more suitable for ET 0 in northern Xinjiang through model selection. Hargreaves and Allen, Berti and Dorji models were selected for comparison with FAO56 P-M models (Table 2). Similarly, four empirical models based on air mass transport, Dalton, Meyer, WMO and Albrecht, were selected to verify the applicability of the four models in northern Xinjiang to find an alternative model with higher applicability under fewer input meteorological parameters ( Table 2). Table 2. Empirical reference evapotranspiration model and parameter calculation.

Meteorological Inputs Equations Proposed by
Based-temperature ET 0 models R a , T max , T min , Mass transfer-based ET 0 models [27] Ra is the theoretical radiation, T max is the daily maximum temperature, T min is the daily minimum temperature, and other parameters are the same as above.

Random Forest-Based Reference Evapotranspiration (ET 0 ) Model
Random forest (RF) is a class of discriminant models that support classification, regression, and multi-classification [28]. It is based on bagging integration on decision trees and introduces random attribute selection of prediction in the training process of decision trees. It establishes a forest in a random way. The forest is composed of many decision trees. There is no correlation between each decision tree, and with the same distribution for all trees, multiple trees are used to train and predict the samples [29].
For a given classifier set, the numerical estimator is h 1 (X), h 2 (X), . . . , h k (X), and a random vector Θ is set; therefore, the tree predictor h (X, Θ) can take on numerical values.
As the number of trees in the forest increases, the training set drawn at random from the distribution of the random vector Y, X, and the marginal functions can be given as: For a large number of decision trees, the following two basic theorems are valid: We define the upper bound for generalization error of RF as follows: Theorem 2. As the number of decision trees increase, all sequences will finally converge to:

Least Square Support Vector Regression
A support vector machine (SVM) was introduced in the 1990s [30] and then extended for support vector regression (SVR), which is a new machine learning algorithm based on the structural risk minimization criterion. The biggest problem of traditional neural networks is that they can minimize the training error. However, it cannot minimize the generalization error in the learning process. SVM is successfully applied for the classification and prediction based on VC dimensions and structural risk minimization. Then, the least square support vector machine (LS-SVR) was proposed by Suykens and Vandewalle [31] of the United States on this basis. The difference between LS-SVR and SVM is that the relaxation variables in the optimization objective are changed from a penalty term to a quadratic term, and the constraint conditions are also changed into only equality constraints, which allows the solution to be realized by the least square method.
The linear regression process is as follows: the n-dimensional inputted vector and 1-dimensional outputted vector: (x k , y k ), x k ∈ R n , y k ∈ R, k = 1, . . . , N, are mapped from the original space to the high-dimensional space F, and the optimal linear regression function is constructed in F. That is: According to the principle of structural risk minimization, there are: where, ω 2 is the complexity of control model, γ is the determining penalty factor, R emp is the loss function, LS-SVM selects the two-norm of e as the loss function. At this time, the optimization problem becomes: The constraint condition is: The model is changed into dual space and Lagrange function is introduced: The partial derivation of the above formula is obtained: For k = 1, . . . , N, the system of linear equations is obtained by eliminating ω and e: Obtained: where, I is the identity matrix. a and b are calculated by the least square method: , the prediction model of LS-SVM is as follows where K(x, x k ) is the kernel function and the kernel function realizes the mapping from low-dimensional space to high-dimensional space and transforms the non-linear problem of low-dimensional space into a linear problem of high-dimensional space.

Bidirectional Long-Term and Short-Term Memory Network
The bidirectional long-term and short-term memory network (LSTM) cell, first proposed by Hochreiter and Schmidhuber [32], is an upgraded version of a recursive neural network that can overcome the problem of vanishing recurrent neural network (RNN) gradients. Standard LSTM networks deal with sequences chronologically, and they ignore future contexts. To overcome this shortcoming, the Bi-LSTM algorithm is introduced to improve the accuracy. The algorithm is a deformed structure of LSTM, and the forward and backward LSTM layers contained in the memory block are able to utilize past and future information. The hidden layer with opposite time series is obtained, and the two hidden layers are connected to obtain the same output result. The hidden layer of Bi-LSTM at time t includes forward → h t and backward ← h t bars: where T denotes the length of the time.

Back Propagation Neural Network Optimized by Genetic Algorithm
The genetic algorithm (GA) is a heuristic algorithm simulating biological evolution that has three genetic operators: selection, crossover and mutation [33]. BPNN consists of a three-layer network structure of an input layer, hidden layer and output layer [12]. With neurons as the basic unit, the BPNN solves the non-linear fitting problem mainly by setting activation functions in the hidden layer and output layer. In this study, a genetic algorithm was used to optimize the structural parameters of the BP neural network. To realize the dynamic updating of the value interval of the number of gene loci in the hidden layer, the crossover operator and mutation operator of the structural parameter chromosome in GA-BP were improved. First, the neuron gene calculation was carried out, and then the value interval of gene loci was derived according to the number of neurons to generate the value of gene loci. The BP structural parameter population was optimized by a genetic operator to obtain the optimal structural parameters ( Figure 2). neurons as the basic unit, the BPNN solves the non-linear fitting problem mainly by setting activation functions in the hidden layer and output layer. In this study, a genetic algorithm was used to optimize the structural parameters of the BP neural network. To realize the dynamic updating of the value interval of the number of gene loci in the hidden layer, the crossover operator and mutation operator of the structural parameter chromosome in GA-BP were improved. First, the neuron gene calculation was carried out, and then the value interval of gene loci was derived according to the number of neurons to generate the value of gene loci. The BP structural parameter population was optimized by a genetic operator to obtain the optimal structural parameters ( Figure 2). Through the correlation analysis of meteorological factors (Figure 3), it can be seen that under the condition of p < 0.05, all parameters are significantly correlated with ET0. Therefore, the following parameters (Table 3) are selected for ET0 prediction. We adopted four optimization algorithms, including RF, LS-SVR, BI-LSTM and GA-BP, and combined with different parameter combinations; a total of 28 reference evapotranspiration models were established. Through the correlation analysis of meteorological factors (Figure 3), it can be seen that under the condition of p < 0.05, all parameters are significantly correlated with ET 0 .
Water 2021, 13, x FOR PEER REVIEW 7 neurons as the basic unit, the BPNN solves the non-linear fitting problem mainly by ting activation functions in the hidden layer and output layer. In this study, a genet gorithm was used to optimize the structural parameters of the BP neural network. T alize the dynamic updating of the value interval of the number of gene loci in the hi layer, the crossover operator and mutation operator of the structural parameter chro some in GA-BP were improved. First, the neuron gene calculation was carried out then the value interval of gene loci was derived according to the number of neuro generate the value of gene loci. The BP structural parameter population was optimize a genetic operator to obtain the optimal structural parameters ( Figure 2). Through the correlation analysis of meteorological factors (Figure 3), it can be that under the condition of p < 0.05, all parameters are significantly correlated with E Therefore, the following parameters (Table 3) are selected for ET0 prediction adopted four optimization algorithms, including RF, LS-SVR, BI-LSTM and GA-BP combined with different parameter combinations; a total of 28 reference evapotrans Therefore, the following parameters (Table 3) are selected for ET 0 prediction. We adopted four optimization algorithms, including RF, LS-SVR, BI-LSTM and GA-BP, and combined with different parameter combinations; a total of 28 reference evapotranspiration models were established. Table 3. Summary of the inputs used for the implementation of each RF, LS-SVR, Bi-LSTM and GA-BP models. Figure 4 shows the flow chart of reference evapotranspiration model research. The data from 2000 to 2014 were taken as the training dataset, and the data from 2015 to 2020 were taken as the test dataset. Four different algorithms were combined to simulate ET 0 values in northern Xinjiang, and the simulation results were evaluated.
Water 2021, 13, x FOR PEER REVIEW 8 of 18 Table 3. Summary of the inputs used for the implementation of each RF, LS-SVR, Bi-LSTM and GA-BP models. Figure 4 shows the flow chart of reference evapotranspiration model research. The data from 2000 to 2014 were taken as the training dataset, and the data from 2015 to 2020 were taken as the test dataset. Four different algorithms were combined to simulate ET0 values in northern Xinjiang, and the simulation results were evaluated.

Performance Evaluation of Models
The assessment of the models was carried out using 5 statistical performance evaluations: the root mean square error (RMSE), the mean absolute error (MAE), the mean bias error (MBE), correlation of determination (R 2 ) and global performance indicator (GPI) expressions, which are as follows:

Performance Evaluation of Models
The assessment of the models was carried out using 5 statistical performance evaluations: the root mean square error (RMSE), the mean absolute error (MAE), the mean bias error (MBE), correlation of determination (R 2 ) and global performance indicator (GPI) expressions, which are as follows: where ET P-M, i is the FAO56 P-M model-observed reference evapotranspiration at the ith time step (mm day −1 ), ET E,i is the empirical or algorithm model (RF, LS-SV, Bi-LSTM and ANN models)-observed reference evapotranspiration at the ith time step (mm day −1 ), and n is the total observations. A i,new is the normalized value of the above four indicators, A i corresponds to the median of the indicators, and α i is equal to 1 for R 2 , equal to −1 for the RMSE and MAE and equal to ±1 for the MBE indicators (the sign is the same as MBE).

Performance Appraisal of Seven Empirical Models (Temperature-Based and Mass
Transfer-Based) for Estimating ET 0 Figure 5 shows that the three temperature-based empirical models have good applicability in northern Xinjiang, China, with R 2 values above 0.9. The Dorji model has the highest coefficient of determination, with R 2 up to 0.9275. However, the empirical model based on air mass transport is not ideal in northern Xinjiang, with an overall R 2 lower than 0.

Components Comparison of the Four Algorithm Models for Estimating ET 0
In this study, seven different parameter combination input modes (Table 3) and four optimization algorithms, including the RF, LS-SVR, Bi LSTM and GA-BP models, were used to rank and evaluate the models in the training period (Table 4) and test period ( Table 5) by introducing the global performance index (GPI). In the training period, the top 10% of the models were focused on the RF model, while the advantage of the RF model disappeared in the test period, which indicates that the RF model has overfitting in the training period, resulting in a significant reduction in the accuracy of the model in the test period. Therefore, the model evaluation in this paper is mainly based on the test period to reduce the impact of the overfitting algorithm on the simulation results. Table 5 shows that the GA-BP5 model has the best overall simulation effect during the test period, with an RMSE of 0.2542 mm·d −1 , MAE of 0.1706 mm·d −1 , MBE of −0.0039 and R 2 of 0.9918. The second best simulation is the LS-SVR6 model, wherein the RMSE is 0.2380 mm·d −1 , MAE is 0.1604 mm·d −1 , MBE is −0.0425 and R 2 is 0.9928. The worst model is the Bi-LSTM1 model, wherein the RMSE is 0.7811 mm·d −1 , MAE is 0.5374 mm·d −1 , MBE is −0.3776 and R 2 is 0.9227. The worst Bi-LSTM1 model has the least number of input parameters, and the best, the GA-BP5 model and LS-SVR6 model, have the most input parameters, which shows that the number of input parameters has a great impact on the model simulation accuracy. In the same algorithm with different parameter combinations, the best results of RF model and LS-SVR model appear in the sixth group of parameter inputs (T, T max , T min , Ra, Rs and U 2 ), and the best results of Bi LSTM model and GA-BP model appear in the fifth group of parameter inputs (T, T max , T min , Ra, Rs and U 2 ); this finding shows that the accuracy of the model cannot be improved by 100% only by increasing the number of input parameters. Table 5. Performance of RF, LS-SV, Bi-LSTM and GA-BP models during the testing period.

Models
Testing Period (2015-2020) Because the influence of meteorological factors on reference evapotranspiration is different, considering RH and U 2 , RMSE decreases by 0.5 mm·d −1 , MAE decreases by 0.32 mm·d −1 , MBE amplitude decreases by 0.3 and R 2 increases by 0.077. Considering the radiation terms Rs and U 2 based on the temperature model, RMSE decreases by 0.52 mm·d −1 , MAE decreases by 0.33 mm·d −1 , MBE amplitude decreases by 0.28 and R 2 increases by 0.076. The results show that the radiation term Rs and aerodynamic term U 2 have a significant effect on reference evapotranspiration, and the reference evapotranspiration model is mainly composed of three parts: the temperature term, radiation term and aerodynamic term.
We obtained the optimal model results under different input combinations (Table 6). Figure 6 shows a comparison between the optimal model and the P-M model under the conditions of various input parameters, which can also support the view that the accuracy of simulation results is improved after considering the radiation term Rs and aerodynamic term U 2 parameters.

Evaluation of Optimal Reference Evapotranspiration Model under Different Time-Scale Conditions
In the verification of the empirical model and optimization algorithm model, this paper takes the data from 2000 to 2020 as samples. A large amount of data and a long time span increase the model verification accuracy. Generally, the estimation of crop water requirements is only based on the reference crop evapotranspiration of 5-7 days. Therefore, the influence of different time scales on the simulation accuracy of the model is also included in the scope of this paper. We divided our data into four time scales, which were

Evaluation of Optimal Reference Evapotranspiration Model under Different Time-Scale Conditions
In the verification of the empirical model and optimization algorithm model, this paper takes the data from 2000 to 2020 as samples. A large amount of data and a long time span increase the model verification accuracy. Generally, the estimation of crop water requirements is only based on the reference crop evapotranspiration of 5-7 days. Therefore, the influence of different time scales on the simulation accuracy of the model is also included in the scope of this paper. We divided our data into four time scales, which were 7 days (1 July 2015 to 7 July 2015), one month (May 2019), one season (summer 2015) and one year (2017), and evaluated them with the Berti model, Hargreaves model and optimal model under seven different parameter input combinations. Figure 7 shows a comparison of Taylor diagrams of models under different time scales. It can be seen from Figure 7, when the time scale is 7 days, the range of standard deviation is 0.5~1.25, while, when the time scale is year, the range of standard deviation is reduced to 0.625~1, and the correlation is increased from 0.75 to 0.95. This result shows that the time scale has a significant effect on the accuracy of the model, and the accuracy of the model increases with increasing time scale. At the annual scale, the seven models based on the optimization algorithm are better. When meteorological data are missing, seven models can be used to estimate ET 0 in northern Xinjiang. Only the LS-SVR6 model and GA-BP7 model are recommended to estimate the local ET 0 at the daily, monthly and seasonal scales.

Discussion
The correlation between reference evapotranspiration model and meteorological factors is the key to determining the input parameters of the optimization model, and many studies have shown that air temperature is the primary factor affecting reference evapotranspiration [4,5,7]. This study shows that ET0 has a very significant positive correlation with the daily average temperature T, the daily minimum air temperature Tmin, the daily maximum air temperature Tmax, the net radiation Ra and the saturated water vapor pressure difference VPD at the level of p < 0.001 (Figure 3), and the saturated water air pressure difference VPD is a function of temperature and humidity. The relative humidity RH was negatively correlated with ET0 at the level of p < 0.001. The whole evapotranspiration process is dominated by temperature and net radiation Ra. This result is consistent with the conclusions of some scholars [22][23][24].
This study found that among the empirical models, the temperature-based model is more suitable for northern Xinjiang, this is consistent with the conclusion of Rachid et al. [34]. That is, the temperature-based model is more suitable for arid and semi-arid areas. The reason is that heat is the main driving factor of water transmission and the key step of the water cycle, especially in arid and semi-arid areas with limited water resources. It is feasible and reliable to apply the optimization algorithm to ET0 prediction. However, under the same parameter input conditions, the estimation accuracy of ET0 at different

Discussion
The correlation between reference evapotranspiration model and meteorological factors is the key to determining the input parameters of the optimization model, and many studies have shown that air temperature is the primary factor affecting reference evapotranspiration [4,5,7]. This study shows that ET 0 has a very significant positive correlation with the daily average temperature T, the daily minimum air temperature T min , the daily maximum air temperature T max , the net radiation Ra and the saturated water vapor pressure difference VPD at the level of p < 0.001 (Figure 3), and the saturated water air pressure difference VPD is a function of temperature and humidity. The relative humidity RH was negatively correlated with ET 0 at the level of p < 0.001. The whole evapotranspiration process is dominated by temperature and net radiation Ra. This result is consistent with the conclusions of some scholars [22][23][24].
This study found that among the empirical models, the temperature-based model is more suitable for northern Xinjiang, this is consistent with the conclusion of Rachid et al. [34]. That is, the temperature-based model is more suitable for arid and semi-arid areas. The reason is that heat is the main driving factor of water transmission and the key step of the water cycle, especially in arid and semi-arid areas with limited water resources. It is feasible and reliable to apply the optimization algorithm to ET 0 prediction. However, under the same parameter input conditions, the estimation accuracy of ET 0 at different scales is different, this conclusion is consistent with the views of some scholars [35][36][37]. The main reason is that the time scale itself is the factor restricting the accuracy of the algorithm. It can be seen from Figure 7 that the model based on the optimization algorithm is significantly better than the empirical, and this conclusion has also been verified in semi-arid areas of Spain [38]. Chen et al. [12] also pointed out that the model based on optimization algorithm is better than the empirical model, and the empirical model based on temperature factor is better than that based on radiation or humidity in the training and research area. It shows that the application of an optimization algorithm can greatly improve the estimation accuracy of the model and make up for the shortcomings of existing empirical formulas.

Conclusions
Based on the daily datasets of eight meteorological stations in northern Xinjiang, China, from 2000 to 2020, seven empirical models (Hargreaves, Berti, Dorji, Dalton, Meyer, WMO, and Albrecht) and 28 models based on optimization algorithms (RF, LS-SVR, Bi-LSTM and GA-BP) were compared. Through the evaluation of eight empirical models, the Berti model and Hargreaves model are better alternative models, while the WMO model and Meyer model are worse alternative models. The global performance index (GPI) is used to rank 28 models based on the optimization algorithm, and the optimal recommendation models of the annual scale and multi-year scale under different input parameter combinations are obtained as follows: LS-SVR1 (input: T, T max , T min , Ra), RF2 (input: T, T max , T min , Ra, Rs), LS-SVR3 (input: T, T max , T min , Ra, RH), LS-SVR4 (input: T, T max , T min , Ra, Rs, RH), GA-BP5 (input: T, T max , T min , Ra, RH, U 2 ), LS-SVR6 (input: T, T max , T min , Ra, Rs, U 2 ) and GA-BP7 (input: T, T max , T min , Ra, U 2 ). Only the LS-SVR6 model and GA-BP7 model are recommended as the best alternative models for ET 0 when the local climate dataset is incomplete at the daily, monthly and seasonal scales, which has high applicability in northern Xinjiang. This study can provide theoretical guidance and technical support for the determination of farmland irrigation systems and water resource planning and management in northern Xinjiang. This method is also applicable to other arid and semi-arid areas.