Hourly soil temperature prediction using integrated machine learning methods, GLUE uncertainty analysis, Taguchi search, and wavelet coherence analysis

In this study, hourly T s variations at 5, 10, and 30 cm soil depth were investigated and predicted for an arid site (Sirjan) and a semi-humid site (Sanandaj) in Iran. Standalone machine learning models (adaptive neuron fuzzy interface system (ANFIS), support vector machine model (SVM), radial basis function neural network (RBFNN), and multilayer perceptron (MLP)) were hybridized with four optimization algorithms (sunflower optimization (SFO), firefly algorithm (FFA), salp swarm algorithm (SSA), particle swarm optimization (PSO)) to improve prediction accuracy and reduce uncertainty. Uncertainty analysis was performed using generalized likelihood uncertainty estimation (GLUE), while wavelet coherence was used to assess interactions between T s and meteorological parameters. For the arid site, ANFIS-SFO (RMSE=1.18 o C, MAE=1.05 o C, NSE=0.93, PBIAS=7%, and R 2 =0.9998) produced the most accurate performance at 5 cm soil depth. At best, hybridization with SFO (ANFIS-SFO, MLP-SFO, RBFNN-SFO, SVM-SFO) decreased RMSE by 5.6, 18, 18.3, and 18.18 % compared with the respective standalone model. At the semi-humid site, all integrated models showed most accurate performance at 10 cm soil depth, with RMSE for the best model (ANFIS-SFO) increasing by 10.5%, and MAE by 10.1%, from 10 to 30 cm depth. GLUE analysis confirmed that integrating optimization algorithms with machine learning models decreased the uncertainty in T s predictions. Wavelet coherence analysis demonstrated that air temperature, relative humidity, and solar radiation, but not wind speed, had high coherence with T s at different soil depths at both sites, and meteorological parameters mostly influenced T s in upper soil layers.

It is difficult, costly, and time-consuming to measure hourly Ts directly (Feng et al. 2019).
Therefore, Ts measurements are not easily available for developing countries, such as Iran (Mehdizadeh et al. 2020a). In addition, Ts varies on the land surface and at different soil depths (Moazenzadeh and Mohammadi 2019), but stations for recording Ts are often sparse and limited.
In recent years, indirect methods, including analytical, numerical, and data-driven models, have been suggested for accurate prediction of Ts. However, analytical models cannot reproduce the complex relationships between Ts and climatological variables, especially for different climates and global scales (Kang et al. 2000;Plauborg 2002). Numerical models use finite volume and finite element methods to solve heat transfer laws, and are highly complex (Badache et al. 2016;Samadianfard et al. 2018b). As an alternative, intelligence machine learning models can be used for predicting Ts and describing the interactions between Ts and other predictor variables, such as meteorological parameters (Moazenzadeh and Mohammadi 2019;Feng et al. 2019). Current advances in intelligence machine learning models are improving their capability for modeling nonlinear relationships (Padarian et al. 2020).
To overcome the existing drawbacks of standalone machine learning approaches and obtain higher prediction accuracy, hybrid models can be used (Riahi-Madvar et al. 2020;. Hybrid methods apply nature-inspired meta-heuristic algorithms to optimize the performance of standalone models. Quick data-processing speed, achieving a global optimum, higher precision than standalone models, and robust generalization are some advantages of integrated machine learning-optimization algorithm techniques (Baydaroğlu and Koçak 2014; Moazenzadeh and Mohammadi 2019). For example, Nahvi et al. (2016) found that self-adaptive evolutionary (SaE) enhanced the estimation accuracy of ELM in predicting daily Ts at different soil depths based on air temperature (Ta), compared with the standalone ELM, GP, and ANN models. Samadianfard et al. (2018b) developed a coupled MLP-firefly algorithm (MLP-FFA) model for monthly Ts estimation in the 0-100 cm soil layer in Adana, Turkey, using the climate variables Ta, atmospheric pressure (Pa), and solar radiation (Rs), and found that it performed better than the standalone MLP model, especially at 20 cm soil depth. Moazenzadeh and Mohammadi (2019) estimated daily Ts at six soil depths in Maragheh, Iran, using two standalone models, elman neural network (ENN) and support vector regression (SVR), and four coupled models, (SVR-FFA), SVR-krill herd algorithm (SVR-KHA), ENN-FFA, and ENN-KHA. They found that all models had their highest accuracy at 10 cm soil depth and that SVR-KHA showed the best estimation accuracy based on Ta, relative humidity (RH), and sunshine hours. Mehdizadeh et al. (2020b) hybridized the ENN model with two algorithms, ant colony optimization (ACO) and gravitational search algorithm (GSA), to predict daily Ts at four soil depths in Rasht and Isfahan, Iran, and found that ENN-GSA produced the most accurate predictions at all depths.
Hourly evaluation of Ts is essential for high-resolution modeling within hydrology, geosciences, ecology, and crop models as a decision-support system to achieve better crop yield (Qi et al. 2016;Sanikhani et al. 2018). Since hourly Ts measurement with instruments is difficult, previous studies have generally used daily or monthly datasets in Ts prediction. To our knowledge, few studies have investigated the precision of conventional machine learning models in predicting Ts at hourly scale (Araghi et al. 2017;Heddam 2019;Feng et al. 2019), while none has investigated model hybridization with integrated optimization algorithms to improve hourly Ts estimation in different climates.
Uncertainty analysis is a robust step to assess the reliability of different models .
Previous studies on Ts prediction are mostly based on standalone machine learning models and consider straightforward prediction, with no explicit investigation of the associated uncertainty and coherence analysis.
In this study, four machine learning models (ANFIS, MLP, SVM, and radial basis function neural network (RBFNN)) were used to predict hourly Ts at three soil depths in two climate zones. The models were then hybridized with four bio-inspired meta-heuristic algorithms (sunflower optimization (SFO), salp swarm algorithm (SSA), FFA, particle swarm optimization (PSO)) and model parameter values optimizing the precision of the models were analyzed. In a further novel approach, generalized likelihood uncertainty estimation (GLUE) was applied in uncertainty analysis to compare the validity and reliability of the models.
Since inferring the relationships between climate variables (Ta, RH, wind speed (U), Rs) and hourly Ts data can lead to a better understanding of soil temperature response to model structure, wavelet coherence analysis was used to assess the behavior of hourly Ts at different soil depths relative to meteorological variables and to identify the most significant meteorological variable for these predictions. To our knowledge based on a literature review, no previous study has applied wavelet coherence analysis to Ts and other climate variables.
Despite the popularity of hybrid optimization models for daily and monthly Ts prediction, few studies have compared hybrid model precision in predicting hourly Ts, especially by analyzing related uncertainty and investigating significant meteorological factors. There have also been few accuracy evaluations of hybrid ANFIS and SVM models in estimating hourly Ts in Iran. The primary aims of this study were therefore to: 1) develop and compare ANFIS, MLP, RBFNN, and SVM hybridized with the SFO, SSA, FFA, and PSO algorithms in predicting hourly Ts at three soil depths in two climate zones; 2) select best values of optimization algorithm parameters, using the Taguchi method; 3) identify model uncertainty at different soil depths using the GLUE approach; and 4) determine the relationships between hourly Ts and Ta, U, RH, and Rs at different soil depths, using wavelet coherence analysis.

2-1-Study area and dataset
Two datasets of measured hourly Ts were chosen to evaluate the hybrid and standalone models.
The first dataset was obtained from an automatic monitoring station in a pistachio orchard (55˚82̍ N, 29˚30̍ E) in Sirjan city, Iran (Fig. 1), where information about Ts at soil depth is essential in disease control, evaporation modeling, irrigation management, and frost protection. Sirjan has an arid climate, with mean annual air temperature of 27˚C and mean annual rainfall of 40 mm. Hourly Ts and meteorological parameters (Ta, Rs, U, RH) in the field were monitored from 12 September 2012 to 17 January 2012. From a total of 1966 data obtained, 1500 hourly data were randomly selected for model learning and the remaining 466 hourly data were used for model testing.
The second Ts dataset was obtained from a synoptic station in Sanandaj (35.33°N, 47.00°E, 1373.4 m height), western Iran (Fig. 1). This region has a semi-humid climate, based on the calculated aridity index using 30-year meteorological data (Zolfaghari et al. 2016), with hot conditions in summer and very cold winters. Mean annual air temperature measured at Sanandaj station is 13.4˚C and mean annual rainfall is 450 mm. Hourly Ts and meteorological parameters (Ta, Rs, U, RH) were recorded from 9 July 2009 to 12 October 2009. A set of 800 hourly Ts data was randomly selected and used for learning the models, while 200 hourly Ts data were used for testing.
Time series plots of Ta and hourly Ts at three depths (5, 10, and 30 cm) in the soil at Sirjan and Sanandaj revealed that temperature variations showed a decreasing trend at all soil depths (Fig. 2).
The greatest difference in temperature between initial and final points was at 5 cm soil depth (15˚C) in Sirjan, while in Sanandaj station, it was 21 ˚C at 10 cm and 19˚C at 30 cm soil depth. Statistical characteristics on the measured meteorological and Ts data are given in Table 1.

2-2-Adaptive neuron fuzzy inference system (ANFIS)
ANFIS is a type of ANN model combined with Takagi-Sugeno fuzzy, i.e., it uses both fuzzy logic and neural networks and can capture the advantages of both (Najafi-Ghiri et al. 2019). Two rules in the "if-then" method for the Takagi-Sugeno model were considered here (Njafi-Ghir et al.
The ANFIS structure contains five layers, where layers 1 and 4 include adaptive parameters and layers 2, 3, and 5 consist of non-adaptive nodes. In layer 1, each node adjusts to a function parameter and produces a value of membership degree ( ( ) −2 ( )). In layer 2, each node is defined and shows a firing strength for each rule (wi). The outcome of multiplying the signals coming into the node is the output of each node. Normalized firing strength (NFS, ̄) is created in the third layer as the ratio between the ith rule's firing strength (FS) and the total sum of FS for all rules. In the adaptive nodes of layer 4, the node function is defined as 4, =̄. In layer 5, the overall output is calculated through combining all incoming signals as = ∑̄.
In this study, the bell membership function was used, due to its smoothness and brevity: where ai, ci, and bi are premise parameters of membership functions.

2-3-Multilayer perceptron (MLP) model
MLP is a feedforward supervised neural network that has been applied successfully to complex problems. The backpropagation learning algorithm is commonly used for training MLP models, but it may get trapped in a local optimum (Pouladi et al. 2019). The MLP model uses multiple layers with a non-linear activation function to learn the relationship between input and output datasets. The first layer (input layer) contains the inputs to the MLP model, while the middle layers (hidden layers) contain a number of neurons. Each neuron performs a weighted summation of inputs. The activation functions are used to calculate the inner product of input parameters and adjustable weight vectors of synapses (Pouladi et al. 2019). More details about MLP can be found in Kisi et al. (2015).

2-4-Radial basis function neural network (RBFNN)
The RBFNN model is a particular type of ANN model that has been widely used for modeling hydrological variables such as streamflow, runoff, temperature, drought, and groundwater level.
The main difference between RBFNN and feedforward multilayer ANN is the transfer function properties employed in the hidden layers (Walczak and Massart 2000). RBFNN contains three layers in its structure. The first layer receives the input vector data, where the hidden layers include several nodes and a nonlinear transfer function (Tayebi et al. 2019). The Euclidean norm (distance) associated with the hidden layer's input vector and center is computed by the activation function related to each hidden neuron. Commonly, the hidden layers contain the nonlinear radial basis Gaussian function and another activation function placed in the last layer. More details about RBFNN can be found in Kisi et al. (2015).

2-5-Support vector machine (SVM)
SVM was initially introduced by Cortes and Vapnik (1995) and has been widely used for both regression and classification analysis, due to its advantage in minimizing model complexity and estimation error simultaneously (Vapnik 1999;Zheng et al. 2020). The SVM model uses different kernel functions to estimate the regressions, implicitly converting inputs to high-dimensional feature space using a hyper-plane . The SVM equation is: where ( ) is a deterministic function, ( ) is a non-linear transferring the input vector, is a vector of weight coefficients, and b is bias.
The b and values are determined by minimizing risk structure-function. To reduce the complexity and obtain a more robust model, slack variables , * can be used in the risk structurefunction equation: where is intensive loss value, is the ith output, and 0 is the cost constant.

2-6-Sunflower optimization algorithm (SFO)
SFO is a new optimization algorithm proposed by Yang (2012) where P is the power of source and ri is the distance between the current and the best i plant.
The redirection of position by the sunflower is computed as: where is the inertial displacement of the sunflower and is constant, * is the best location of the sunflower, and is the current location of the plant, and (‖ + −1 ‖) is the pollution probability, where i th sunflower cross-pollinates with another near i-1 th sunflower to generate a new individual in an updated location.
It is essential to limit the maximum step of individuals by: where dmax is the maximum step, is the upper bound value, Xmin is the lower bound value, and Npop is the number of sunflowers in the overall population.
Finally, the next population is updated as: where ⃗ +1 is the position of sunflower i+1.

2-7-Firefly algorithm (FFA)
FFA was first developed by Yang (2010), inspired by the light emission capability of fireflies. The attractiveness of one firefly to another is related to its brightness to neighboring fireflies, where a less bright firefly is attracted to a brighter one (Alor et al. 2019;Riahi-Madvar et al. 2020). The excellent information-sharing mechanism is one of the advantages of FFA. The firefly mechanism is formulated as: where ( ) is the attractiveness of a firefly, 0 is the attractiveness of firefly at r=0, is the absorption coefficient in the range of 0 and 1, ( ) is the light intensity, 0 is the light intensity at r=0, and the distance between two fireflies i, j at locations of xi and xj can be defined as = ‖ − ‖. The firefly i is attracted by firefly j and updates its location as: where is randomization parameter, and is a vector of random parameters.

2-8-Salp swarm algorithm (SSA)
SSA is a new meta-heuristic optimization algorithm, inspired by the collective behavior of salps where is the leader's position, yi is the food source, and are the lower and upper bounds, respectively, r3 and r2 are random numbers, and r1 is computed as: where l is the existing iteration, and L is the maximum number of iterations. SSA uses the parameter r1 to increase stability in exploration and exploitation capability. For each follower, the position is updated as: where is the location of jth salp in the ith dimension.

2-9-Particle swarm optimization (PSO)
PSO is a well-established stochastic/random search approach related to the swarm and inspired by the particle's social behavior. Each particle in the swarm protects the updating of search behavior in accordance with learning experiences of all other particles. In each generation, information is integrated by particles to set the velocity on every dimension. The position and velocity of particles are updated as: , +1 = * , + 1 * * ( , − , ) + 2 * * ( , − , ) , +1 = , + , +1 where , +1 is the velocity of i th particle at iteration t+1, is the weight inertia, 1 and 2 are acceleration coefficients, , +1 is the location of i th particle at iteration t+1, , is the optimal location experienced by particles, and , is the optimal location experienced by any particle.

2-10-Hybridizing ANFIS, MLP, RBFNN, and SVM with optimization algorithms
In this study, SFO, FFA, SSA, and PSO were used to fine-tune the parameters of the ANFIS, MLP, RBFNN, and SVM models. The input data were entered for training. The optimization algorithms then started with random selection of agents (particles, fireflies, salps, and sunflowers) and finally Robust design of random parameters included in the optimization algorithms is important to enhance model performance. In a novel approach, Taguchi search was used in this study for selecting the best values of these random parameters. It is a powerful advanced technique that uses orthogonal array and signal to noise ratio (S/N) to minimize the number of experiments and greatly decrease the time, cost, and effort in finding optimal parameters of algorithms (Zhang et al. 2015).
An orthogonal array table is created by calculating total degree of freedom (DOF) based on the combined degree of freedom of all parameters (Canbolat et al. 2019;Bademlioglu et al. 2020).
Based on number of levels (L) and number of parameters (NV), the total number of experiments needed to find the optimum level of parameters is computed. The minimum number of experiments (N) is calculated as = 1 + ( − 1). Thus e.g., for the four parameters with four levels in the PSO algorithm, at least 13 experiments should be conducted to discover the optimal values of PSO parameters, while without the Taguchi method, the total number of experiments would be 4 4 . Thus, the Taguchi orthogonal array of L16 (4 3 ) was established.
The S/N ratio is calculated as: where n is the number of the case, and y is the performance characteristic performance value (the higher the S/N value, the better the ratio). S/N was used as a performance characteristic of each parameter in this study.
The framework and flowchart of the work in this study is presented in Fig. 3, which also shows the general framework of the 16 hybrid models used to predict hourly Ts (ANFIS-SFO, ANFIS- Level 1: Random sampling is used to create a large number of sampling sets from a prior distribution of input data. Level 2: Likelihood value is calculated from model runs and compared with a certain threshold value to evaluate each input parameter as behavioral (value above the threshold likelihood) or nonbehavioral (value below the threshold likelihood). Behavioral parameters are retained to judge the models. Likelihood is calculated as: where N is adjustable parameter, 2 is the error variance for the i th model, is the variance of observations, is the parameter set, and ( | ) is the likelihood measure for the i th model calculated with the observations Y.
Level 3: Simulation weights for behavioral parameter sets are rescaled and the cumulative weighted distribution of estimations is used in quantile estimation for uncertainty prediction: where is a likelihood weight and n is the number of data.
Two indices, p and r, are used to quantify model uncertainty, where p is percentage of bracketed observations at 95% prediction uncertainty (95PPU) and r is mean 95PPU range, separated using the standard deviation of observations : where k is the number of observations, is the standard deviation of observations, and and are the upper and lower boundary of the 95% prediction uncertainty, respectively. Lower values of r and higher values for p indicate lower uncertainty. In this study, the spectral responses of models in the GLUE technique were used to assess variations in r and p.

2-12-Wavelet coherence analysis in interaction analysis
In reality, Ts is tele-connected to different meteorological variables such as Ta, RH, U, and Rs.
Wavelet coherence analysis was used here to assess these interactions. It calculates the coherence value of cross-wavelet transform among two time series in the time-frequency domain (Lee and Kim 2019). Detailed theory on wavelet transform can be found in Grinstead et al. (2004). In the present case, wavelet coherence between Ts and meteorological variables was calculated as (Lee and Kim 2019): where x and y are two time series, s is the smoothing over both scale and time, S is smoothing operator, is the wavelet coefficient of soil temperature and meteorological parameters in the time series x, is the wavelet coefficient of soil temperature and meteorological parameters in the time series y, ( ) is the cross-wavelet power spectrum of x and y, is the occasion, Sscale is smoothing along with the wavelet, and is smoothing in time (see Fig. 4).

3-Results and discussion
The best optimization algorithm parameters, extracted with Taguchi search, for training the standalone and hybrid models to predict Ts based on climate variables (Ta, Rs, U, RH) are presented in section 3.1. In section 3.2, model performance and the results obtained at different sites and soil depths are compared. The uncertainty associated with the simulations, determined using the GLUE approach, is presented in section 3.3. Wavelet coherence results for the interconnections between Ts and climate variables are summarized in section 3.4.

3-1-Optimal model parameters derivation by Taguchi search
The optimization parameters of the meta-heuristic algorithms combined with SFO, SSA, FFA, and PSO at different levels, and the associated S/N values, are shown in Table 2. The values of population size, c1, c2, inertia weight, and γ in optimization algorithms were derived based on the Taguchi method. After constructing an orthogonal array and calculating S/N ratio for each algorithm parameter, the optimal values were chosen (highest S/N ratio) ( Table 2). The effect of each optimization algorithm process parameter was divided into four levels. For the PSO algorithm, the mean value of RMSE for population size varied from 1.12 to 1.54, inertia weight from 1.23 to 1.76, c1 from 1.12 to 1.45, and c2 from 1.37 to 1.45, and the best conditions were obtained at a value of 2, 4, 1, and 1, respectively. The greatest variation in S/N values was seen for inertia weight, i.e. PSO was mainly affected by inertia weight, which should be the main focus in model development. The most critical parameter for SSA, SFO, and FFA was r3, population size, and γ, respectively. The optimal value of r3, r2, and population size for SSA, determined by Taguchi search, was 0.5, 0.8, and 400 respectively. For SFO, population size of 400 was most effective, while for FFA γ had an optimal value of 0.8. These results indicate that Taguchi search can effectively and systematically provide robustness values in selection of optimization parameters, replacing the trial-and-error approach. It can be used in future studies of metaheuristics to improve the ability and applicability of models.

3-2-Performance analysis of models
The performance of standalone ANFIS, MLP, RBFNN, and SVM and the hybrid models was evaluated based on statistical criteria and by comparison of predictions against measured hourly Ts at 5, 10, and 30 cm soil depth (Tables 3 and 4). Soil temperature was used as the target parameter, and corresponding Ta, U, Rs, and RH as input data for all models. At the training level, integration of SFO with ANFIS, SVM, RBFNN, and MLP led to better results. The standalone ANFIS model outperformed the other standalone SVM, RBFNN, and MLP models. All models provided the highest accuracy at 10 cm soil depth, and the lowest at 5 cm soil depth. The prediction accuracy in terms of RMSE, MAE, NSE, and PBIAS for the training level was lower for SVM-PSO than for the other hybrid ANFIS, SVM, MLP, and RBFNN models.
It should be noted that the measured data at 5 cm depth were from the arid Sirjan site and the measured data at 10 and 30 cm soil depths from the semi-humid Sanandaj site. Table 4 shows the performance of the models in predicting hourly Ts in terms of RMSE (°C), MAE (°C), PBIAS (%) (optimal value=0), and NSE (optimal value=1) in the testing phase. Scatter plots and correlation comparisons between simulated and measured Ts' values at 5, 10, and 30 cm soil depths are presented in Fig. 5, 6, and Nanda et al. (2020) also found that SVM had the lowest precision of all machine learning models tested in hourly and half-hourly Ts prediction. The main reason for the low accuracy of SVM may be non-linearity between several analyzed parameters.
Comparison of the results at 10 and 30 cm depth indicated that the error in Ts estimation increased with depth, indicating that meteorological parameters mostly influenced soil temperature in the upper layers (Kazemi et al. 2018;Behmanesh and Mehdizadeh, 2017). Similarly, Nahvi et al. (2016) found that the precision of machine learning models in estimation of Ts in soils in Turkey and Iran declined from 10 cm to 100 cm depth, but increased from 5 to 10 cm depth. Behmanesh and Mehdizadeh (2017) also found that the accuracy of machine learning models increased between 5 and 10 cm depth. In the present study, the highest accuracy was observed at 10 cm depth for all hybrid and standalone models, as also found by e.g., Behmanesh and Mehdizadeh (2017) and Kisi et al. (2017). Among the models they studied, Kisi et al. (2017) found that ANFIS achieved the highest precision at 10 cm depth (RMSE=1.29) based on 25-year monthly dataset at 10, 50, and 100 cm depth.
Based on the results in  Samadianfard et al. (2018) found that integrated MLP-FFA and SVM-FFA models had significantly higher Ts prediction accuracy compared with the standalone MLP and SVM models at 5, 10, and 20 cm soil depth.
In the scatterplots of all test dataset, all estimations fell on the 1:1 line, with R 2 values >0.99 for all models (Figs. 5-7). The R 2 value increased from 0.9958 to 0.9998, 0.9934 to 0.9993, and 0.9950 to 0.9998 for all hybrid and standalone models at 5, 10, and 30 cm depth, respectively. All models had the most accurate performance for 10 cm soil depth and the worst for 5 cm soil depth . This reflects the complexity and rapid changes in Ts at shallow depths. Overall, the ANFIS-SFO predictions were closest to the observed data, with the highest R 2 and lowest error at each soil depth, confirming the suitability of ANFIS-SFO for estimating hourly Ts in regions with arid (Sirjan) and semi-humid (Sanandaj) climates.
From Table 4 and Figures 5-7, it is evident that use of the all meta-heuristic algorithms (SFO, SSA, FFA, PSO) improved the accuracy of the corresponding standalone models. This may be due to more accurate searching and finding the best solution in the local and global spaces. This is in agreement with findings by Mehdizade et al. (2020), Moazenzadeh and Mohammadi (2019) and Shamshirband et al. (2020) that hybrid models show higher accuracy than conventional models.
Based on results from the present study, the ANFIS model integrated with a meta-heuristic optimization algorithm is highly recommended for estimating Ts at different soil depths in different climates.

3-3-Uncertainty analysis of models
Spectral representations of uncertainty values of the models at the testing level for 5, 10, and 30 cm soil depth, determined using the GLUE approach, are shown in Fig. 8. and RBFNN models. This indicates that using meta-heuristic algorithms such as SFO results in certain parameter adjustment in AI models and reduces the case-specificity, improving the generality of the models. In the present case, hybridization of the model decreased the uncertainty in soil temperature predictions and provided more accurate and reliable depth analysis of temperature.
In the uncertainty analysis, an increase in the number of data points increased the r values for the different models (Fig. 8). In general, SVM had the highest r and lowest p, indicating a high level of uncertainty. The standalone and hybrid MLP models produced more accurate predictions than the standalone and hybrid RBFNN models. Overall, the optimization algorithms improved the performance and reliability of the standalone models in terms of r and p values. SFO outperformed the other optimization algorithms in reducing uncertainty in estimation of hourly Ts at all soil depths.
In general, the evaluation criteria values and GLUE uncertainty analysis of models confirmed that the hybrid ANFIS-SFO model had most accurate performance in predicting hourly Ts based on Ta, RH, U, and Rs for semi-humid and arid climates of Iran. The results demonstrated the suitability of SFO for hybridizing machine learning models.

3-4-Wavelet coherence analysis
To assess the importance of input vector selection in Ts simulation, the effect of four meteorological parameters was considered using wavelet coherence analysis. The wavelet coherence between Ts and other climate variables (Ta, U, Rs, and RH) at 5, 10, and 30 soil depths is shown in Figures 9 and 10 for the period 0-1944 hours at 5 cm soil depth and 0-1000 hours at 10 and 30 cm soil depth. There was high coherence between Ta and Ts at periodicity of 4-8 hours in the time domain 486-709 hours at 5 cm soil depth at Sirjan (Fig. 9). Overall, Ta showed the highest correlation with Ts, confirming findings by Kisi et al. (2015) and Nanda et al. (2020).
Analysis of the diagram indicated no significant coherence between U and Ts at different time periodicities, but there was a significant relationship between Ts and Rs at periodicity of 28-34 hours, and between Ts and RH at periodicity of 28-34 hours, in the time domain 1252-1944 hours ( Fig. 9) To further analyze the results of wavelet coherence, the values were plotted against periodicity of Ts. Fig. 10 shows the coherence at 10 and 30 cm soil depth for Sanandaj station. The maximum coherence between Ts at 10 and 30 cm soil depth and Ta was 0.8 and 0.7, respectively, at periodicity 8-16 hours in the time domain 0-1000 hours. There were no significant coherences between Ts at different periodicities and U at 10 and 30 cm soil depth. The maximum coherence between Ts at periodicity 4-8 hours and Rs in the time domain 0-1000 hours was 0.9 and 0.7 at 10 and 30 cm depth, respectively. The highest coherence found (0.9) was at 10 cm soil depth between RH and Ts at periodicity 8-16 hours in the time domain 0-1000 hours. These results indicate that wavelet coherence analysis can provide more quantitative evidence on periodicity and interaction of Ts with climate parameters at different soil depths.

4-Conclusions
Soil temperature strongly affects soil biological processes, irrigation scheduling, plant growth, the environment, and water resource management, but measuring soil temperature is time-consuming and costly. Accurate, reliable and certain prediction of soil temperature using models is thus necessary in agricultural, environmental and geosciences practices.
In this study, the machine learning models ANFIS, MLP, RBFNN, and SVM were used for estimating hourly soil temperature at various soil depths, based on meteorological data from an arid site (Sirjan) and a semi-humid site (Sanandaj) in Iran. Four meta-heuristic optimization algorithms (SFO, FFA, PSO, SSA) were used to hybridize and optimize the standalone models.
Parameter selection of the meta-heuristic models was optimized by Taguchi search, while uncertainty in model results was analyzed by the GLUE approach. Wavelet coherence analysis was used to infer the interactional periodicity of variables.
ANFIS was the most accurate standalone model, while SVM was the least accurate, for both sites and all soil depths (5, 10, 30 cm). The optimization algorithms SFO and SSA were best in enhancing performance for all standalone models. Uncertainty results indicated that the bracketed observed values of the best model, ANFIS-SFO, were in the range 94-96%.
Coherence wavelet analysis indicated weak effects of wind speed, but strong effects of air temperature, relative humidity, and solar radiation, on temperature in upper layers. Periodicity of Ts with other effective parameters was inferred.
The approach developed in this study proved reliable and accurate for hybrid meta-heuristic model development and assessment. It is applicable in regions with similar climate conditions to the study sites. Future studies should consider the effects of climate scenarios on soil temperature for the base period and future period. ANIFS, SVM, RBFNN, and MLP can also be hybridized with multiobjective optimization algorithms to determine the optimal values of hyperparameters and appropriate inputs to the models.

Ethics declarations Ethical approval
Not applicable.

Consent to publish
All authors have agreed with the content and all have given explicit consent to publish.

Data availability
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Code availability
The code is available from the corresponding author on reasonable request.

Funding
Not applicable.