## Abstract

Reference evapotranspiration (ET_{0}) is an important parameter to characterize the hydrological water cycle and energy balance. An extremely heavy rainstorm occurred in Zhengzhou City, Henan Province on 20 July 2021, causing heavy casualties and economic losses. One of the important reasons for this rainstorm was abnormal water circulation. The purpose of this study is to estimate ET_{0} accurately and avoid extreme disasters caused by abnormal water cycles. This study compared and analyzed the accuracy and robustness of ET_{0} prediction based on the improved Levenberg–Marquardt (L-M) model based on artificial neural network and the genetic algorithm-backward neural network (GA-BP) model. The model uses seven weather stations in Zhengzhou, including mountain climate and plain climate. By utilizing the Pearson correlation analysis technique, six distinct input scenarios were identified, and the efficacy of the model was assessed using evaluation metrics, including RMSE, MAE, NSE, and SI. The results show that the estimation accuracy of the L-M model is better than that of the GA-BP model; when the number of input meteorological parameters is the same, the combined simulation effect including wind speed is the best; the *R*^{2} of L-M3 and L-M4 are 0.9285 and 0.9675, respectively; models can accurately estimate ET_{0} with limited data.

## HIGHLIGHTS

Estimation of ET

_{0}in Zhengzhou by improved L-M and GA-BP models based on neural networks.Six different input scenarios were introduced according to the Pearson correlation analysis.

Evaluate the robustness of different models in different input scenarios.

Temperature and wind speed provide more information for estimating ET

_{0}.This study provides methodological support for predicting ET

_{0}in regions lacking meteorological data.

## INTRODUCTION

Evapotranspiration is an important link in the hydrological cycle and energy conversion process (El-Shafie *et al.* 2014; Mattar 2018), which can transport surface water into the atmosphere while also carrying away net radiated energy from the surface, and is a key element in the global energy balance and water cycle. Since evapotranspiration is difficult to measure directly and is a relatively complex nonlinear natural process, reference evapotranspiration is often used to estimate actual evapotranspiration (Jovic *et al.* 2018). Therefore, accurate estimation of reference evapotranspiration is of great scientific significance for regional water resources planning and management, hydrological research, and prediction of drought and flood disasters (Shiri *et al.* 2014; Wu *et al.* 2019).

There are many methods to estimate reference evapotranspiration. Lysimeter is a direct method to measure evapotranspiration. However, due to the complex operation and high maintenance cost, this method is rarely used in practice. After that, scholars put forward empirical formulas based on meteorological data. For example, empirical formulas for estimating ET_{0} based only on meteorological data such as temperature and radiation (Valipour 2014; Mattar *et al.* 2016): Blaney-Criddle (BC) (Blaney & Criddle 1950), Priestley-Taylor (PT), Hargreaves-Samani (Hargreaves & Samani 1985), etc*.* Although these formulas only use a small amount of meteorological data in the calculation process, the simulation accuracy is greatly affected by the geographical location and climate change of the study area, and the prediction results are often overestimated or underestimated (Berti *et al.* 2014; Feng *et al.* 2017; Liu *et al.* 2017). The Penman–Monteith (PM) equation is a calculation method based on the water vapor diffusion theory and the energy balance theory, which takes into account both radiation and aerodynamic effects and is applicable to the global scale. The PM equation has been recommended by the FAO as a standard method for estimating ET_{0} and calibrating other models (Allen *et al.* 1998; Wang *et al.* 2019; Tikhamarine *et al.* 2020; Chen *et al.* 2020). However, the calculation process of this method is complex and requires a large amount of reliable and high-quality meteorological data (Roy *et al.* 2021), so the generalization of the PM equation in regions with incomplete meteorological data is limited (Valipour 2015; Yin *et al.* 2020; Wu *et al.* 2020a, 2020b).

Evapotranspiration is a complex nonlinear dynamic process that requires methods beyond empirical models. In recent years, the rapid development of machine learning algorithms provides new ideas for estimating ET_{0} (Granata 2019; Roy *et al.* 2020; Ahmadi *et al.* 2021; Chia *et al.* 2021; Yan *et al.* 2021). Huo *et al.* (2012) estimated ET_{0} in northwest China using the ANN model, and the results showed that relative humidity, wind speed (at a 2 m height), maximum and minimum temperature were the key meteorological factors for estimating ET_{0} by the ANN model. Feng *et al.* (2016) used extreme learning machine (ELM), wavelet neural network (WNN), and generalized artificial neural network (GANN) models to estimate ET_{0}, and verified it in humid areas in southwest China. The results showed that the simulation accuracy of ELM and GANN models was better than that of the WNN model. Ferreira *et al.* (2019) evaluated the performance of artificial neural network (ANN) and support vector machine (SVM) for estimating daily ET_{0} across Brazil. The results showed that the performance of ANN and SVM for ET_{0} estimation is superior to other empirical formulas. Tang *et al.* (2018) used SVM and GANN to estimate the actual evapotranspiration of dryland maize with no and partial mulch coverage, respectively. The results show that the performance of the GANN model is better than the SVM model. Malik *et al.* (2019) predicted the monthly mean ET_{0} in the central hilly region of India based on multiple linear regression (MLR), radial basis neural network (RBNN), multi-layer perceptual neural network (MLPNN), and self-organizing mapping neural network (SOMNN) models. It is found that the prediction results were better when the input parameters included three (maximum air temperature, wind speed, and solar radiation) and five (minimum and maximum air temperatures, relative humidity, wind speed, and solar radiation) variables. Dong *et al.* (2022) used convolutional neural network (CNN), ELM, and multiple adaptive regression spline curve (MARS) to analyze the spatio-temporal variation of ET_{0} in China. The results showed that the performance of the CNN model in estimating ET_{0} was better than that of ELM and MARS models.

Although neural network models have been proven to be effective in estimating ET_{0}, there are still some problems in practical use. In the process of neural network training, the input data and output data establish a corresponding relationship, which mainly shows the generalization ability of the model. However, the general neural network model will have the problem of overfitting in the training process, which will reduce the generalization ability of the model and produce the suboptimal solution (EL-Shafie *et al.* 2014). In addition, due to the different sizes and structures of the available data input from different models (Granata 2019), problems such as slow operation speed, weak global search ability, and easy convergence to local extremum will also occur in the running process of the model (Sun *et al.* 2016), which will affect the overall performance of the model. To overcome these problems, an improved Levenberg–Marquardt (L-M) model based on neural networks is introduced in this study. The L-M model can be regarded as an improved form of the Gauss-Newton method, which has both the local characteristics of the Gauss-Newton method and the global characteristics of the gradient method (Sajedi *et al.* 2021). It can adjust the iteration speed and accuracy in the training process, so as to have more powerful generalization ability. Compared with the general neural network algorithm, it has a faster iterative convergence speed and higher accuracy (Ali *et al.* 2021). In recent years, the L-M model is usually used in artificial intelligence computing (Vidmar *et al.* 2020), electromagnetic simulation (Darisma & Marwan 2019), new energy battery development, and medical trials (Smirnova *et al.* 2019; Cheema *et al.* 2020). All these results show that the L-M model has good predictive ability. In addition, the GA-BP model is proposed as a competitive alternative to the L-M model, a BP neural network optimized by a genetic algorithm. Through the implementation of the genetic algorithm's global optimization capability, optimal initial weights and thresholds for the BP neural network were attained. These optimal values were then utilized as the initial weights and thresholds of the BP neural network to avoid falling into local minima during subsequent training (Zheng *et al.* 2018; Tu *et al.* 2020; Yin *et al.* 2022). The GA-BP model has the advantages of faster convergence speed, fewer iterations, and higher simulation accuracy, and is widely used in the prediction of long-term weather temperature (Dou & Sun 2021), local scour depth of bridge pier (Dong *et al.* 2020), concrete strength (Tu *et al.* 2020), flow and pressure of urban water supply network (Xia *et al.* 2021). It is worth noting that L-M and GA-BP models are rarely used in predicting ET_{0}.

In general, no specific model is suitable for all situations, and the performance of different models depends on how much effective information is provided by the input data and the structure of the input data (Granata 2019). The purpose of this study is to use the improved L-M algorithm and the GA-BP algorithm based on the traditional neural network to estimate ET_{0}, respectively. By inputting different combinations of meteorological variables and data structures, the accuracy and robustness of two different models to estimate daily ET_{0} under different input scenarios are compared and analyzed. In addition, this study is based on the data of seven weather stations in Zhengzhou City, in order to obtain a model that can use fewer meteorological parameters, it does not need to consider hydrological principles, and can accurately estimate ET_{0}.

## MATERIALS AND METHODS

### Case sites and data collection

In this study, the daily data of seven meteorological stations in the Zhengzhou Administrative region were used. The meteorological stations were evenly distributed in the area of Zhengzhou, including both the western mountain climate and the eastern plain climate, which was a good representation. Daily meteorological data of seven meteorological stations from 2000 to 2019 were selected, including wind speed (*U*_{2}, m/s), hours of sunshine (SS_{H}, h), relative humidity (*R*_{H}, %), maximum temperature (*T*_{max}, °C), and minimum temperature (*T*_{min}, °C). The daily meteorological data (2000–2019) collected were divided into two datasets: the training dataset (2000–2017) and the validation dataset (2018–2019). The model parameters are trained on the training dataset and the weights are estimated. The validation dataset further verifies the accuracy of the training model.

### Data input combination

*et al*. 2020). The redundant information is removed to ensure that the input meteorological data provides more independent information for the model and improve the operation efficiency and accuracy of the model. The calculation formula is:where Sim(

*A*,

*B*) is the correlation between

*A*and

*B*variables;

*A*is the

_{i}*i*-th sample value of

*A*variable;

*B*is the

_{i}*i*-th sample value of

*B*variable; is the average value of variable

*A*; is the average value of variable

*B*;

*i*is the

*i*-th sample of the variable; and

*n*is the total sample size.

≥ 0.8 shows a high positive correlation between two variables; 0.5 ≤ < 0.8 indicating that the two variables are moderately correlated; < 0.5 indicating low correlation between the two variables. In this paper, variables with high or medium correlation were selected as the input of the prediction algorithm (Table 1).

Meteorological parameters . | T_{min}
. | T_{max}
. | R_{H}
. | SS_{H}
. | U_{2}
. |
---|---|---|---|---|---|

Correlation coefficient | 0.7424 | 0.8020 | − 0.5312 | 0.6953 | 0.5646 |

Meteorological parameters . | T_{min}
. | T_{max}
. | R_{H}
. | SS_{H}
. | U_{2}
. |
---|---|---|---|---|---|

Correlation coefficient | 0.7424 | 0.8020 | − 0.5312 | 0.6953 | 0.5646 |

Pearson correlation coefficients of *T*_{min}, *T*_{max}, *R*_{H}, SS_{H}, and *U*_{2} were all above 0.5, belonging to medium and high correlation levels. It is worth noting that the correlation coefficient of *R*_{H} is negative, indicating that *R*_{H} hinders the evapotranspiration process. In order to reflect the actual evapotranspiration process more truly, the influence of *R*_{H} on the actual evapotranspiration was also considered in this study. In addition, according to relevant research results (Traore *et al.* 2010; Yang *et al.* 2019), ET_{0} value is greatly affected by temperature, so temperature is a necessary input when selecting input meteorological factors. Therefore, the model input combination is divided into six scenarios (Table 2). Scenarios 1–3 have three input factors, scenarios 4–5 have four input factors, and scenario 6 has five input factors. The input meteorological factors in this study are all parameters easily obtained by meteorological stations, which to some extent solves the problems of backward monitoring facilities and difficulty in obtaining complex parameters in some regions.

. | T_{min}
. | T_{max}
. | R_{H}
. | SS_{H}
. | U_{2}
. |
---|---|---|---|---|---|

Model1 | √ | √ | √ | ||

Model2 | √ | √ | √ | ||

Model3 | √ | √ | √ | ||

Model4 | √ | √ | √ | √ | |

Model5 | √ | √ | √ | √ | |

Model6 | √ | √ | √ | √ | √ |

. | T_{min}
. | T_{max}
. | R_{H}
. | SS_{H}
. | U_{2}
. |
---|---|---|---|---|---|

Model1 | √ | √ | √ | ||

Model2 | √ | √ | √ | ||

Model3 | √ | √ | √ | ||

Model4 | √ | √ | √ | √ | |

Model5 | √ | √ | √ | √ | |

Model6 | √ | √ | √ | √ | √ |

### Research methods

#### Penman–Monteith formula

*et al.*1998) fully considers solar radiation, energy balance, aerodynamics, and other principles. The calculation results are accurate and do not require calibration (Roy

*et al.*2021). Therefore, FAO recommends the PM equation to calculate ET

_{0}value, and the calculated ET

_{0}value is generally used as the standard value. The equation is given as:where ET

_{0}is the reference evapotranspiration, mm/d;

*R*

_{n}is net radiation, MJ/(m

^{2}·d);

*G*is soil heat flux, MJ/(m

^{2}·d), when calculating evapotranspiration,

*G*= 0;

*U*

_{2}is the daily average wind speed 2 m above the ground, m/s;

*T*is the average temperature, °C;

*E*

_{s}is the saturated vapor pressure, kPa;

*E*

_{a}is the actual vapor pressure, kPa;

*γ*is the hygrometer constant, 0.065 kPa/°C; and

*δ*is the slope of the saturation steam pressure–temperature curve, kPa/°C.

#### L-M model

The commonly used neural network algorithm is the gradient descent method. The parameters move in the opposite direction of the error gradient to reduce the error function until the minimum value is obtained. The complexity of the calculation is mainly caused by the calculation of partial derivatives. However, the linear convergence rate based on the gradient descent method is very slow, and the L-M algorithm is an improved form of the Gauss-Newton method, which has both the local characteristics of the Gauss-Newton method and the global property of the gradient method. The L-M algorithm is much faster than the gradient method because it uses the approximate second derivative information. The following is a brief description of the L-M algorithm.

*I*is the identity matrix;

*w*represents the vector composed of the weights and thresholds of the

^{k}*k*-th iteration, and

*w*

^{k}^{+1}represents the vector composed of the new weight and threshold.

*et al.*2016). For a given , if Δ

*w*can reduce the error exponential function

*V*(

*w*), then is divided by the factor ; if the error index function

*V*(

*w*) increases, is multiplied by the factor . The flowchart of the L-M algorithm is shown in Figure 2. The steps of the L-M algorithm:

- (a)
Given the allowable training error

*ɛ*, , and constants, and initializing the weight and threshold vectors, set*k*= 0, = ; - (b)
Calculate the network output and the error exponential function

*V*(*w*);^{k} - (c)
Calculate the Jacobian matrix

*J*(*w*);^{k} - (d)
Calculate the Δ

*w*; - (e)
If

*V*(*w*) <^{k}*ɛ*, go to (g); - (f)
If

*V*(*w*) >^{k}*ɛ*, take*w*^{k}^{+1}=*w*+ Δ^{k}*w*, as the weight and threshold vector, calculate the error index function*V*(*w*^{k}^{+1}), if*V*(*w*^{k}^{+1}) <*V*(*w*), set^{k}*k*=*k*+ 1, = , go to (b), otherwise = , go to (d); - (g)
Operation is complete.

_{0}. The normalized equation is:where

*x*is the normalized value,

_{s}*x*is the measured value of a factor in the sample,

*x*

_{max}is the maximum value of the sample data, and

*x*

_{min}is the minimum value of the sample data.

#### Genetic algorithm-back propagation neural network

Genetic algorithm is an iterative optimization model based on Darwinian evolution theory and genetic evolution theory. Each species undergoes natural evolution through successive generations, inheriting genetic traits from their parents while also experiencing certain genetic changes. In each generation, individuals with genetic traits better suited for the environment have a higher chance of survival. Through the process of natural selection and evolution over multiple generations, the offspring that remain possess the genetic makeup that best fits their environment. As a chromosome, it is a survival of the fittest problem-solving process, through the evolution of chromosome generations, to obtain the optimal or satisfactory solution to the problem. The BP neural network algorithm is a local search optimization method, the essence of which is the gradient descent method. The advantage lies in strong adaptive ability, good generalization ability, and fault tolerance ability. The disadvantage is that the algorithm is extremely sensitive to the initial weight and threshold, its results are easy to converge to local minima, and the convergence speed is slow (Liu *et al.* 2019). Combined with BP neural network and genetic algorithm, the population search method is used to optimize the weight and threshold of the neural network, which can better overcome the local defect of the BP neural network which tends to the local optimal solution (Saleh *et al.* 2016).

*x*is the number of neuron nodes in the input layer, and is the number of input meteorological factors;

*z*is the number of neuron nodes in the output layer. Because the output layer only has reference evapotranspiration,

*z*= 1;

*y*is the number of hidden layer neuron nodes.

*ξ*, which are the sum of the absolute errors between the estimated value and the calculated value of the training sample, as shown in the following formula:where

*k*is the calculated value of ET

_{j}_{0}in the

*j*-th day,

*o*is the estimated ET

_{j}_{0}value in the

*j*-th day.

In this study, both L-M and GA-BP models are improved algorithms based on neural network, including input layer, hidden layer, and output layer in the process of machine learning. The number of input layer and output layer is relatively easy to determine, while the number of hidden layer is a very important factor in the operation process (Faris *et al.* 2019). Improper selection of hidden layer will not only affect the efficiency of machine learning but also directly affect the simulation results. However, in many pieces of literature (El-Shafie *et al.* 2014; Maroufpoor *et al.* 2020), there is no specific rule to determine the number of hidden layers. In this study, the number of hidden layer neurons is determined by trial and error to ensure the simulation effect, operation efficiency, and robustness of the model. The L-M model and GA-BP model adopt the same input structure. It is worth noting that the iteration types of the L-M model and GA-BP model are different, so the number of iterations is also different. It is found (Maroufpoor *et al.* 2020) that when the number of iterations reaches a certain number, increasing the number of iterations has little impact on improving the simulation accuracy of the model, and even increases the running burden of the model. Through experiments, the iteration times of L-M model and GA-BP are 1,500 times and 10 times, respectively.

### Model performance evaluation

_{0}(Ferreira

*et al.*2019; Ahmadi

*et al.*2021). RMSE is used to measure the deviation between the simulated value and the standard value. The smaller the RMSE, the more accurate the model prediction. MAE is the average value of the absolute difference between the simulated value and the standard value, which can avoid the problem of the error canceling each other, and does not need to consider the sign, so it can accurately reflect the size of the model error. As the primary estimation parameter in this study, NSE quantifies the goodness-of-fit between simulated and observed values. Serving as an indicator of the model's simulation accuracy, NSE values range from (−∞, 1), with values closer to 1 indicating higher levels of accuracy. Conversely, NSE values less than 0 suggest poor performance. SI is a parameter used to evaluate the performance of the model. When SI < 0.1, the model performance is excellent, when 0.1 < SI < 0.2, the model performance is good, when 0.2 < SI < 0.3, the model performance is average, and when SI > 0.3, the model performance is poor. The expression is as follows:where

*N*is the number of samples,

*S*is the simulated value of ET

_{i}_{0}in the

*i*-th day,

*O*is the calculated value of ET

_{i}_{0}in the

*i*-th day, and is the average value of the calculated value of ET

_{0}.

## RESULTS AND ANALYSIS

### Model evaluation

In this study, the performance of the model was evaluated according to the different numbers of input data, including scenarios 1–3 with three meteorological data input, scenarios 4–5 with four meteorological data input, and scenario 6 with five meteorological data input.

_{0}estimation results of L-M and GA-BP models under different input scenarios. When the three meteorological parameters are input, and the input parameters are

*T*

_{min},

*T*

_{max}, and

*U*

_{2}, the simulation accuracy of the two models is the highest, indicating that wind speed is the main factor affecting the evapotranspiration process, which is consistent with the practical significance that wind speed accelerates air flow and facilitates evapotranspiration. When the input meteorological parameters are

*T*

_{min},

*T*

_{max}, and

*R*

_{H}, the simulation accuracy of the model is the worst. It is consistent with the negative correlation between

*R*

_{H}and ET

_{0}in the Pearson correlation analysis mentioned above, indicating that

*R*

_{H}plays an obstructive role in the natural evapotranspiration process (Maroufpoor

*et al.*2020; Ahmadi

*et al.*2021). For different types of models, the simulation results of the L-M model are better than the GA-BP model. The RMSE, MAE, and NSE of L-M3 are 0.4209 mm, 0.3300 mm, and 0.9271, respectively. Compared with GA-BP3, RMSE and MAE are reduced by 7.11 and 3.45%, respectively. NSE increased by 1.27%. For the same model, the simulation accuracy varies with the change of input parameters. When the input parameters include wind speed, the simulation accuracy of the model reaches the maximum, and the simulation effect of scenario 3 is the best, followed by scenario 2, and scenario 1 is the worst. RMSE and MAE of L-M3 decreased by 20.99 and 17.04%, respectively, while NSE increased by 5.21%. Compared with GA-BP1, RMSE and MAE of GA-BP3 decreased by 19.05 and 15.12%, respectively, while NSE increased by 4.15%. It can be seen that under different input scenarios of the same model, the L-M model is more sensitive to meteorological parameters, and the simulation effect is better than the GA-BP model. Figures 3 and 4 show the comparison of predicted values, residuals, and standard values of L-M and GA-BP models in scenarios 1–3. It can be clearly seen that

*R*

^{2}of L-M1 and GA-BP1 are 0.8868 and 0.8772, respectively. Although

*R*

^{2}performs well, both of them are less than 0.9. When the input meteorological parameters include SS

_{H}and

*U*

_{2}, respectively, the

*R*

^{2}of the two models are above 0.9, and in the same model, the input combination of scenario 2 and scenario 3 produces similar simulation results. In general, the simulation accuracy of the L-M model is better than the GA-BP model, and the

*R*

^{2}of L-M3 is 0.9285, which is 1.28% higher than GA-BP3 (

*R*

^{2}). In addition, the residual figure shows that the residual value of the L-M model is smaller than that of the GA-BP model. The error range of L-M1 (b) is −1.5 ∼ 2.0 mm, and the error range of GA-BP1 (b) is −1.5 ∼ 3.0 mm, with the largest error range. The error ranges of L-M2 (b) and GA-BP2 (b) are basically the same as those of L-M3 (b) and GA-BP3 (b), but the L-M model is superior to the GA-BP model. It shows that model type and input structure are important factors affecting ET

_{0}.

Model . | Structure . | Training (2000–2017) . | . | Validating (2018–2019) . | . | ||||
---|---|---|---|---|---|---|---|---|---|

RMSE . | MAE . | NSE . | SI . | RMSE . | MAE . | NSE . | SI . | ||

L-M1 | 3-6-1 | 0.4843 | 0.3543 | 0.8977 | 0.1472 | 0.5327 | 0.3978 | 0.8812 | 0.1567 |

L-M2 | 3-6-1 | 0.4741 | 0.3529 | 0.9128 | 0.1359 | 0.4257 | 0.3352 | 0.9240 | 0.1264 |

L-M3 | 3-6-1 | 0.4728 | 0.3516 | 0.9148 | 0.1345 | 0.4209 | 0.3300 | 0.9271 | 0.1238 |

L-M4 | 4-8-1 | 0.3064 | 0.2167 | 0.9604 | 0.0917 | 0.2897 | 0.2175 | 0.9655 | 0.0852 |

L-M5 | 4-8-1 | 0.4009 | 0.2904 | 0.9321 | 0.1200 | 0.3647 | 0.2910 | 0.9452 | 0.1073 |

L-M6 | 5-10-1 | 0.1857 | 0.1970 | 0.9854 | 0.0556 | 0.1397 | 0.1845 | 0.9920 | 0.0411 |

GA-BP1 | 3-6-1 | 0.4960 | 0.3626 | 0.8953 | 0.1490 | 0.5597 | 0.4027 | 0.8790 | 0.1646 |

GA-BP2 | 3-6-1 | 0.4803 | 0.3560 | 0.9048 | 0.1421 | 0.4557 | 0.3431 | 0.9145 | 0.1340 |

GA-BP3 | 3-6-1 | 0.4800 | 0.3552 | 0.9049 | 0.1419 | 0.4531 | 0.3418 | 0.9155 | 0.1333 |

GA-BP4 | 4-8-1 | 0.3580 | 0.2622 | 0.9459 | 0.1072 | 0.3397 | 0.2540 | 0.9525 | 0.0999 |

GA-BP5 | 4-8-1 | 0.4010 | 0.2995 | 0.9321 | 0.1201 | 0.4099 | 0.3161 | 0.9308 | 0.1206 |

GA-BP6 | 5-10-1 | 0.3013 | 0.2278 | 0.9617 | 0.0902 | 0.2622 | 0.2087 | 0.9717 | 0.0771 |

Model . | Structure . | Training (2000–2017) . | . | Validating (2018–2019) . | . | ||||
---|---|---|---|---|---|---|---|---|---|

RMSE . | MAE . | NSE . | SI . | RMSE . | MAE . | NSE . | SI . | ||

L-M1 | 3-6-1 | 0.4843 | 0.3543 | 0.8977 | 0.1472 | 0.5327 | 0.3978 | 0.8812 | 0.1567 |

L-M2 | 3-6-1 | 0.4741 | 0.3529 | 0.9128 | 0.1359 | 0.4257 | 0.3352 | 0.9240 | 0.1264 |

L-M3 | 3-6-1 | 0.4728 | 0.3516 | 0.9148 | 0.1345 | 0.4209 | 0.3300 | 0.9271 | 0.1238 |

L-M4 | 4-8-1 | 0.3064 | 0.2167 | 0.9604 | 0.0917 | 0.2897 | 0.2175 | 0.9655 | 0.0852 |

L-M5 | 4-8-1 | 0.4009 | 0.2904 | 0.9321 | 0.1200 | 0.3647 | 0.2910 | 0.9452 | 0.1073 |

L-M6 | 5-10-1 | 0.1857 | 0.1970 | 0.9854 | 0.0556 | 0.1397 | 0.1845 | 0.9920 | 0.0411 |

GA-BP1 | 3-6-1 | 0.4960 | 0.3626 | 0.8953 | 0.1490 | 0.5597 | 0.4027 | 0.8790 | 0.1646 |

GA-BP2 | 3-6-1 | 0.4803 | 0.3560 | 0.9048 | 0.1421 | 0.4557 | 0.3431 | 0.9145 | 0.1340 |

GA-BP3 | 3-6-1 | 0.4800 | 0.3552 | 0.9049 | 0.1419 | 0.4531 | 0.3418 | 0.9155 | 0.1333 |

GA-BP4 | 4-8-1 | 0.3580 | 0.2622 | 0.9459 | 0.1072 | 0.3397 | 0.2540 | 0.9525 | 0.0999 |

GA-BP5 | 4-8-1 | 0.4010 | 0.2995 | 0.9321 | 0.1201 | 0.4099 | 0.3161 | 0.9308 | 0.1206 |

GA-BP6 | 5-10-1 | 0.3013 | 0.2278 | 0.9617 | 0.0902 | 0.2622 | 0.2087 | 0.9717 | 0.0771 |

The same results were obtained when four meteorological parameters were input. For different types of models, the simulation effect of the L-M model is still better than that of the GA-BP model in the same input scenario. The RMSE, MAE, and NSE of L-M4 were 0.2897 mm, 0.2175 mm, and 0.9655, respectively. Compared with GA-BP4, RMSE and MAE decreased by 14.72 and 14.37%, respectively, while NSE increased by 1.36%. For the same model, compared with L-M5, RMSE and MAE of L-M4 decreased by 20.56 and 25.25%, respectively, while NSE increased by 2.10%. Compared with GA-BP5, RMSE and MAE of GA-BP4 decreased by 17.13 and 19.65%, respectively, while NSE increased by 2.33%. This may be related to the inclusion of both temperature and wind speed in the input scenario (Maroufpoor *et al.* 2020).

*R*

^{2}of the two models are above 0.94 when the four meteorological parameters are input. For the same model, the

*R*

^{2}of L-M4 is 0.9675, and that of L-M5 is 0.9603. The error range of L-M4 is −1.0 ∼ 1.0 mm, and the maximum error value of L-M4 is reduced by 50% compared with L-M5. The GA-BP model yielded similar results. For different types of models, the

*R*

^{2}of GA-BP4 and GA-BP5 were 0.9542 and 0.9603, respectively. Compared with GA-BP4 and GA-BP5, L-M4 and L-M5 increased by 1.37 and 1.71%, respectively. The residual range of the GA-BP model is basically concentrated in the range of −1.5 ∼ 2.0 mm. It indicates that when four meteorological factors are input, both models show good prediction performance, but the L-M model is better than the GA-BP model on the whole, and the two models are more sensitive to the input combination including wind speed.

*y*=

*x*) (Maroufpoor

*et al.*2020), with the highest fitting degree and correlation coefficient

*R*

^{2}of 0.9925. It is worth noting that the scatter distribution of GA-BP model simulation results (Figure 8) showed an obvious sag trend to the lower right, which is related to the simulation accuracy of the model, indicating that under the same input of meteorological parameters, GA-BP6 simulation results overestimate ET

_{0}(Huo

*et al.*2012). In addition, it can be seen from the residual figure that the error range of L-M6 is −0.6 ∼ 0.3 mm, which is closer to zero compared with GA-BP6 (−1.0 ∼ 1.0 mm). It shows that the prediction performance of the L-M model is better than that of the GA-BP model.

### Comparison of the models

It can be seen from the above analysis that scenarios 3, 4, and 6 are the best input combinations with different numbers of meteorological parameters for L-M and GA-BP models, respectively. Therefore, this study mainly conducts model comparison and analysis of the above three input scenarios.

_{0}value predicted by different models and the daily ET

_{0}value calculated by the PM method had the same intra-year variation trend. From April to September, the temperature is high, the sunshine duration is sufficient, and the ET

_{0}value is large. From October to December and January to March of the next year, the temperature is low, the sunshine hours are short, and the ET

_{0}value is small. For the same model, the coincidence degree of predicted value and specific collimation is L-M3 < L-M4 < L-M6, and the GA-BP model has the same conclusion. For different models, the predicted values and standard values of the three scenarios conform to the L-M > GA-BP model.

*et al.*2020), the first quartile (Q1), the third quartile (Q3), the quartile interval (IQR), and the rectangular portion of the display median. The most important thing in the box diagram is the calculation of relevant statistical points, which can be achieved by the calculation of percentiles. Among them, Q3 is more important than Q1 in the error judgment of data distribution, because the error of Q3 accounts for 75%, while the error of Q1 only accounts for 25%. The smaller the interquartile range (IQR), the more concentrated the dataset is, the better the model simulation effect is; otherwise, the model simulation effect is poor. The residual boxplots of ET

_{0}predicted by the model in scenarios 3, 4, and 6 are shown in Figure 10. The Q1 of L-M3, L-M4, and L-M6 is −0.3383, −0.2104, and −0.0917, respectively, which is better than that of GA-BP3 (−0.3463), GA-BP4 (−0.1830), and GA-BP6 (−0.1110). For ΔQ3, compared with the GA-BP model, ΔQ3 is 0.0022, 0.0778, and 0.1758, respectively. Similarly, for IQR, the L-M model spacing range is small, and L-M6 error distribution is close to zero. It can be seen that the performance of the L-M model is better than that of the GA-BP model, and the L-M6 model has the best simulation performance.

*et al.*2012; Maroufpoor

*et al.*2020). Although the SI value of the GA-BP model is also in the excellent range, the overall simulation effect is not as good as that of the L-M model. In addition, the SI values of the two models decrease with the increase of the number of meteorological parameters, indicating that with a certain model, the more meteorological parameters involved in the estimation, the better the simulation effect of the model. However, when the input parameters are limited, L-M3, L-M4, or L-M6 can be selected as the estimation model according to the actual accuracy requirements.

According to the above analysis, the simulation effect of the L-M model is better than that of the GA-BP model. The main reasons are as follows: on the one hand, the operation structure of the model is different, which may be because the GA-BP model only optimized the selection of initial weights and thresholds, and didn't consider the problem of overfitting in the iterative process. However, the L-M model combines Gauss-Newton and the steepest descent method in the operation process and adaptively adjusts the damping factor to achieve the convergence characteristics, which has a high iterative convergence speed. At the same time, the problems such as overfitting in the operation process are solved (El-Shafie *et al.* 2014), and the performance of estimating ET_{0} value is improved. On the other hand, models have different sensitivities to input data. Under the same input combination, the L-M model is more sensitive to meteorological parameters. As can be seen from Table 3 and Figures 3–8, all evaluation indexes under the same input combination are better than the GA-BP model, and simulation accuracy improvement of different input combination models is greater in the L-M model than in the GA-BP model.

## CONCLUSION

For the planning and management of the regional hydrological cycle and water resources, accurate estimation of reference evapotranspiration is crucial. If given the right meteorological information, machine learning algorithms are a potent tool that can reliably estimate ET_{0}. To estimate the reference evapotranspiration of Zhengzhou in the usual temperate monsoon climate, respectively, better L-M and GA-BP models based on machine learning algorithms were introduced in this study. Six alternative input scenarios made up of typical meteorological factors were presented using the Pearson correlation analysis. It compares, analyzes, and assesses how well different models execute ET_{0} estimate under various input conditions. The results lead to the following conclusions:

- (1)
The L-M model combines Gauss-Newton and the steepest descent method. In all input scenarios, the estimation accuracy, applicability, and robustness of the L-M model are better than those of the GA-BP model.

- (2)
The L-M model is more sensitive to the input meteorological combination. When the three meteorological parameters are input, the combination of

*T*_{max},*T*_{min}, and*U*_{2}has the best simulation effect, and the correlation coefficient*R*^{2}of L-M3 and GA-BP3 is 0.9285 and 0.9168, respectively. When four meteorological parameters are input, the combination of*T*_{max},*T*_{min}, SS_{H}, and*U*_{2}has the best simulation effect, and the correlation coefficient*R*^{2}of L-M4 and GA-BP4 is 0.9675 and 0.9542, respectively. This indicates that temperature and wind speed can provide more effective information when the model estimates ET_{0}. - (3)
When the number of input meteorological parameters is different, the GA-BP model also shows good prediction performance, but the L-M model shows the optimal prediction ability. Among them, L-M3 (

*T*_{min},*T*_{max},*U*_{2}), L-M4 (*T*_{min},*T*_{max}, SS_{H},*U*_{2}), and L-M6 (*T*_{min},*T*_{max}, SS_{H},*U*_{2},*R*_{H}) have better simulation results, with*R*^{2}values of 0.9285, 0.9675, and 0.9925, respectively. Therefore, within the allowable error range, the L-M model can be used as a reliable model to estimate ET_{0}with a small number of meteorological parameters.

Due to data limitations, further research is needed to determine the applicability of the L-M model in other regions of the world. However, this study has shown that the L-M model has good prediction ability, and its accuracy and robustness are better than those of the GA-BP model. The effects of human activities and climate change on the model's performance were also taken into consideration. Therefore, while the applicability of the L-M model in other regions of the world should be studied, it is also recommended to improve the model's performance by combining it with a deep learning model to increase its reliability.

## AUTHOR CONTRIBUTIONS

For this research paper with several authors, a short paragraph specifying their individual contributions was provided. C.N. developed the original idea and contributed to the research design for the study. S.J. gave some revised opinions for the article. S.L. and C.L. were responsible for data collection and charting. S.S. provided some guidance for the writing of the article. C.H. provided guidance and improved suggestions. All authors have read and approved the final manuscript.

## FUNDING

This work was supported by the National Natural Science Foundation of China, Project No. U2243219, Key projects of National Natural Science Foundation of China (51979250, 51739009), and Key Research and Promotion Projects (technological development) in Henan Province, grant number 222102320455.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*Journal of Irrigation and Drainage*

**35**(S1), 112–115. DOI: 10.13522/j.cnki.ggps.2016.z1.028.