Study site

This study was carried out using data from automated weather stations (AWS) located in the state of Campeche, Mexico (Figure 1). The predominant climate is warm sub-humid, which is seen in 92 % of the territory and 7.75 % presents warm humid climate, localized in the eastern part of the state, and in the northern part a percentage of 0.05 % with semi-dry climate. The highest temperature is over 30 °C and the lowest 18 °C. The mean annual temperature is 26 to 27 °C. The rainfall is abundant to very abundant during the summer. The total annual precipitation ranges between 1 200 and 2 000 mm, and in the northern region, of semi-dry climate, it is around 800 mm annually (^{INEGI, 2017}).

Figure 1 Location of the weather stations in the state of Campeche (Mexico) used in this study. 

The databases for every 10 min were obtained from the automated weather stations of the National Water Commission (Comisión Nacional del Agua, Conagua) in Mexico. Table 1 shows the geographic information of the weather stations used in this study, as well as the periods of registry time from each weather station.

Management of missing data and quality

The databases were processed every 10 minutes, detecting missing series of time through algorithms implemented in the Microsoft Excel® software, which were later refilled using the interpolation technique called Piecewise Cubic Hermite Interpolating Polynomial (PCHIP). For a more detailed description, see ^{Salazar, Ureña and Gallego (2010)}, and ^{Torrente-Cantó (2018)} . Once the data series were completed, daily databases were constructed. Likewise, the data were analyzed to identify atypical values, where the values above three standard deviations of the mean were marked as atypical. The data signaled as atypical were analyzed, and if an atypical extreme inferior value was associated to a rainfall event it was not eliminated; on the contrary case, it was eliminated, with the objective of having functional models even during rainy season. In the case of the atypical extreme superior values, they were eliminated; in both cases the technique of PCHIP interpolation was used to complete them. Figure 2 shows the box and whisker plots of the weather stations implied in the models of this study. The whisker represents the minimum and maximum of the variables. Concerning the maximum temperature, the mean varies between 33 and 34 °C with maximum values between 42 and 43 °C (except the station of Cd. del Carmen), the minimum value of the maximum temperature ranges between 22 and 25 °C. The minimum temperature has a mean between 20 and 24 °C, with maximum values between 25 and 28 °C and minimum between 12 and 17 °C. Figure 2 also shows the atypical values, generally associated to the minimum temperature.

Figure 2 Box and whisker plots of the temperature variable in the weather stations analyzed. 

FAO 56 PM (ETo-FAO56PM) equation

The FAO56PM equation is the standard model used to estimate accurately the ETo, proposed by the Food and Agriculture Organization of the United Nations (FAO). It incorporates thermodynamic and aerodynamic aspects, taking into consideration many meteorological parameters related to the process of evapotranspiration such as net radiation, air temperature, vapor pressure deficit, and wind speed. It has been shown to be a relatively accurate method under different conditions or regions (^{Allen et al., 1998}). These aspects have been incorporated into the following equation:

where R_n = net radiation on the surface (MJ m^-2 day^-1); G = heat flow of the soil (MJ m^-2 day^-1); T_med = mean temperature of the air at 2 m of height (°C); u₂ = wind speed at 2 m of height (m s^-1); e_s = vapor pressure at saturation (kPa); e_a = real vapor pressure (kPa); ∆ = slope of the vapor pressure curve (kPa °C^-1); γ = psychrometric constant (kPa °C^-1).

In this study, the ETo-FAO56PM method was used to evaluate the conventional and artificial intelligence methods.

Hargreaves and Samani Equation

The HS model is considered an alternative model to estimate the ETo when only the temperature records are available in the study site; it is one of the methods that has been used consecutively because of its simple implementation and the accuracy of its results (^{Gong et al., 2016}; ^{Shiri, 2017}). Equation (2) of the Hargreaves and Samani model is structured as follows:

where ETo = Reference evapotranspiration (mm day^-1); K_HG = is an empirical coefficient, which was initially established at 0.0023 but has been recalibrated according to the place used; T_med = mean temperature; T_max = maximum temperature; T_min = minimum temperature; R_a = extraterrestrial solar radiation. R_a was calculated in function of the day of the year, the site’s latitude, and the solar angle according to the equation proposed by ^{Allen et al. (1998)}.

Camargo equation

The Camargo model is a modification of the Thornthwaite (TH) equation; it is a model based on the climate variable of temperature. Camargo substituted the value of the mean temperature of the Thornthwaite equation by the mean effective temperature (T_ef) (^{Camargo, Marin, Sentelhas, & Picini, 1999}). Equation (3) of the Camargo model is structured in the following way:

where ETo = reference evapotranspiration (mm day^-1); K_CA1 and K_CA2 = empirical coefficients, where their original values are 16 and 0.36, respectively, and should be calibrated according to the place of use; I = annual heat index; a= empirical exponent in function of I, N = maximum hours of sunshine; (K_CA2 * (3T_max-T_min) = effective temperature, replacing the mean temperature in the Thornthwaite equation. The value of I is defined as the sum of 12 values of monthly heat indexes, as shown in the following equation:

where T_medj = mean monthly temperature (°C).

and:

The value of a ranges from 0 to 4.25, while the annual heat index I varies from 0 to 160.

Parameter adjustment in conventional methods

The conventional methods based on temperature should be adjusted to the local conditions before being used (^{Almorox et al., 2018}) to obtain good estimations of ETo; therefore, the original coefficients of the equations were calibrated by using non-linear regression techniques through the Levenberg-Marquardt algorithm.

Support Vector Machines (SVM)

The technique of support vector machines (SVM) was developed by ^{Vapnik (2000)} and is one of the approaches based on automated learning. It is a robust supervised learning technique to solve classification and regression problems applied to large sets of complex data with noise; it selects a unique hyperplane of separation of each class and the basic idea is to map the data x in a space of characteristics of high dimension through non-linear mapping and to make a linear regression in this space. During training, only the examples found in the margin of separation are considered, called support vectors. They are applied successfully in regression problems, generally called SVR (support vectors regression), using SVM for a set of data {(X_i , Y_i )} N/i = 1, where X_i is the entry vector, Y_i is the exit value, and N is the total number of sets of data through mapping of X in a characteristic space through a non-linear function φ(x) to later find a regression function ( ^{Fan et al., 2018a}; ^{Mehdizadeh, Behmanesh, & Khalili, 2017}; ^{Quej, Almorox, Arnaldo, & Saito, 2017}; ^{Topi & Vanita, 2017}; ^{Wen et al., 2015}):

where φ (x) is the function of non-linear mapping. ω is a weight vector and b is a bias value, they are the parameters of the regression function, which can be calculated by minimizing the following function of regularized risk:

where the term 12 ω² improves the generalization of the SVM model, normalizing the degree of complexity of the model; C is a parameter of positive compensation that determines the degree of error in the problem of optimization chosen by the user; (ε) is the loss function by Vapnik (size of the tube of the SVM model) and is defined as:

That is, if the difference between predicted and measured values is less than ε, then the loss is equal to 0. If the predicted values are inside the tube, the loss error is equal to 0. For the rest of the predicted points found outside the tube, the loss is equal to the difference between the predicted value and the radius ε of the tube. For the detection of atypical values, the width values ξ and ξ ⃰, measure from top to bottom in the tube of ε.

Because both variables acquire positive values, the risk should be minimized with the following equation:

Subject to yi - ωϕxi - bi ≤ ε + ξi ωϕxi+ bi - yi ≤ε + ξi*ξ, ξi*≥0 [/p]

where C∑i=1n(ξi + ξi*) control the degrees of empirical risk.

The SVM model estimates the regression in function of a series of Kernel functions, which convert original entry data of smaller dimensions to a space of characteristics of larger size in an implicit way. Among the most used Kernels, there is polynomial SVM (SVM-Poly) and the radial base function (SVM-RBF), whose Kernel parameters ought to be adjusted previously through an algorithm. For example, the optimal Kernel parameters and the SVM model are generally obtained using the grid search method (^{Mehdizadeh et al., 2017}):

Implementation of the SVM Model in Estimating the ETo

In this study, the SVM model to estimate the ETo was built by using the R software (^{RDevelopment, 2009}). As entry variables, the meteorological data of T_max, T_min and R_a , were used, and as target variable the values of ETo-FAO56PM (Eq. 1). For training and validation of the SVM model, the R software was used together with the LIBSVM 3.1 package (^{Chang, Lin, & Tieleman, 2013}). The Kernel Radial Base Function (RBF) was used (Eq. 10) to resize the data. With the aim of avoiding the over-adjustment and increasing the performance of the SVM model to estimate the ETo, the parameters ε, C and ɣ of the SVM, and the Kernel Radial Base Function were optimized through the genetic algorithm (GA), using cross-validation (CV= 5 folders) (^{Quej et al., 2017}; ^{Shrestha & Shukla, 2015}), and varying the parameters ε = 0.002 to ε = 2, C = 0.0001 to C = 10, and ɣ = 0.0001 to ɣ = 2. The GA was implemented using the R software together with the e1079 and Caret library; 60 % of the data were used during the stage of training and 40 % in the stage of validation.

The parameters optimized by the GA used in the training of SVM are shown in Table 2.

Genetic Expression Programming (GEP)

Genetic expression programming (GEP) was presented by ^{Ferreira (2001)} . It is a branch of the evolutionary algorithms that has the capacity of modelling dynamic and nonlinear processes. It is an algorithm that belongs to the family of traditional genetic algorithm (GA) and genetic programming (GP). It can emulate biological evolution based on computer programming to solve a problem defined by the user.

GEPs are considered a hybrid between GA and GP. They use genetic programming for the solution of the problem in tree shape, where there are two types of nodes:

In GEP, individuals are codified first as linear chains of fixed length as in GA. Then, they are expressed as nonlinear entities of different sizes and shapes, as in GP (^{Ferreira, 2001}). In addition, a set of terminals (coefficients and predictors), functions and mathematical operators are used in the GEP to estimate the dependent variable (^{Mehdizadeh et al., 2017}), creating functions randomly and selecting those that present a better adjustment to the experimental results, allowing the generation of algorithms and mathematical expressions automatically for the solution of problems (^{Mattar, 2018}; ^{Shiri, 2017}).

Implementation of the GEP model in the estimation of the ETo

In this study the implementation of the GEP technique was carried out by using the GenexprooTols V. 5.0 software, the entry variables are values of meteorological data of T_max, T_min, R_a and ETo-FAO56PM values as target variable. The arithmetic operators and mathematical functions implemented in the program were {+,-, ×, ÷, x, x3,  x2, x3, Inx, eᵡ,sin⁡x,cos⁡x, Arctanx}, recommended for hydrological studies ((^{Mattar, 2018}; ^{Mehdizadeh et al., 2017}; ^{Shiri, 2017}). Of the data, 70 % were used for the stage of training and 30 % for the validation, using cross-validation (CV =5 folders) to avoid over-adjustment. The GEP parameters used in this study are shown in Table 3 (^{Shiri et al., 2014}).

XGBoost

It is one of the most important and potent algorithms of “Machine Learning” created by ^{Chen and Guestrin (2016)} , used for the analysis of regression problems and statistical classification, which produces a complex prediction model from the assembly of decision trees (simple models), in a context of supervised learning.

The model is based on the theory of slope increase, which is why predictions of several “weak” apprentices (models whose predictions are slightly better than random assumptions), are combined to develop a “strong” apprentice. These “weak” apprentices are combined by following a strategy of gradual learning, avoiding over-adjustment and optimizing computer resources. This is obtained by simplifying all the functions that allow combining predictive and regularization terms, but which at the same time maintain an optimal computational speed during the whole processing. At the beginning of the calibration process, a “weak” apprentice is adjusted to the entire data space, and then, a second apprentice is adjusted to the residues of the first. This adjustment process of a model to the residues of the previous one continues until some stoppage criterion is reached (minimization of the root-mean-square error). The result is a type of weighed mean of individual predictions of each weak student. Traditionally, the regression trees are selected as “weak” apprentices (^{Fan et al., 2018a}). Under this context the XGBoost model is based in the following objective function: loss + regulator:

where l is the predictive term and Ω the regularization term. The loss function of the predictive term can be specified by the user. The regularization term is obtained as an analytical expression based on the number of tree leaves and the punctuation of each leaf. The key point of the XGBoost calibration process is that both terms are reordered in the last instance in the following expression:

where G and H are obtained from the expansion of the Taylor series of the loss function, λ is the regularization parameter, and T is the number of leaves in a tree. This analytical expression of the target function allows a fast scanning from left to right of the possible divisions of the tree, but always taking into account the complexity.

XGBoost has a broad range of adjustment parameters. In addition, the flexibility of the algorithm improves when giving the user the chance of including some auto-defined parameters, such as the loss function or the measurement used for validation and trial ^{Urraca, Antonanzas, Antonanzas-Torres, & Martinez-De-Pison, 2017}).

Implementation of the XGBoost model in the estimation of the ETo

For the implementation of the XGBoost model, as the first step, the hyper-parameters were optimized: nrounds, max_depth, eta, gamma, colsample bytree, min_child_weight, and subsample (Table 4) using the Caret library of Software R (^{Fan et al., 2018a}); second, the XGBoost model was adjusted using the “Xgboost” library of Software R, using cross-validation (CV = 5 folders) to avoid over-adjustment. Of the data, 70 % were used for the stage of training and 30 % for validation. The entry variables to the model are the values of meteorological data of T_max, T_min and R_a , the ETo-FAO56PM values as target variable.

Statistical analysis

Four statistical indicators are used in this study to evaluate the performance of conventional and artificial intelligence models; these indicators are: coefficient of determination (R²; Eq. 13), root-mean-square error (RMSE; Eq. 14), mean absolute error (MAE; Eq. 15), mean bias error (MBE; Eq. 16):

where n is the number of comparisons, P_i and O_i are estimated and observed values of ETo-FAO56PM, respectively. P_avg is the average of estimated values of ETo, O_avg is the average of observed values of ETo. The units of ETo are shown in mm d^-1.

R² is commonly used to estimate the performance of hydrological models; it represents the fraction of estimated values that are the closest to the line of measurement data. Values of the coefficient of determination close to 1 indicate more efficient models and the regression line adjusts better to the data. The RMSE is a measurement used frequently to compare errors of prediction in different models, and the lower its value, the better the predictive capacity of a model will be in terms of its absolute deviation. MAE is the sum of absolute values of the errors divided by the number of observations, and it is used frequently to measure how close the estimated values are to the observed values. The MBE provides information about the tendency of the model to overestimate or underestimate the variable, quantifying the systematic error of the model.

The evaluation of uncertainty of statistical indicators R², RMSE, MAE and MBE was carried out through the construction of bootstrap Intervals of Confidence (BIC) at 95 % of the level of trust, and for such a purpose the non-parametric percentile bootstrap method was used as the technique of resampling using B = 10 000 replicas with replacement with the purpose of inducing more accuracy in the estimations (^{Efron, 1992}).

The BICs offer a way of estimating with high probability a range of values where the value of the parameter is found (statistical indicator).

The standard error of the distribution (SE) and the value of each statistical bootstrap indicator were also computed, calculating the standard deviation and the mean of the B replicas.