EvatCrop: a novel hybrid quasi-fuzzy artificial neural network (ANN) model for estimation of reference evapotranspiration

Study area

For the case study, the proposed model was applied in the Terai agro-climatic region of North Bengal, India. This area (9,614 sq. km) is located in the northern part of the state of West Bengal in India, from 25.94 ${ }^{\circ }$N to 27.00 ${ }^{\circ }$N and 88.41 ${ }^{\circ }$E to 89.87 ${ }^{\circ }$E (Mandal et al., 2022). Figure S1 (ESRI, 2009), presents a geographic map of the study area. This region is characterized by humid to hot sub-humid weather conditions. The minimum and maximum monthly average temperature during the rabi season (November to April) is about 13.9 °C and 27.2 °C, respectively, and during the kharif season (May to October), it is about 23.2 °C and 29.7 °C, respectively. The average monthly rainfall in this region is around 37.6 mm during the rabi season and 543 mm during the kharif season. The average monthly wind speed during the rabi season is between 1.0 and 1.2 m/s, whereas during the kharif season, it is between 0.88 and 1.97 m/s. The northern part is covered with forests, and the rest of the area (6,055 sq. km) contains agricultural lands with sandy clay loam to silty clay loam soil textures. The livelihood of the people of this area is primarily based on agriculture, resulting in intensive agricultural activities. Most cultivations in this area depend on groundwater irrigation during the rabi season. Rural farmers use traditional flood irrigation methods without precisely estimating the actual crop water requirement. Such indiscriminate use of groundwater resources threatens groundwater conservation. A recent study reported that a large part of this region is now subjected to moderate to extreme groundwater stress during the rabi season (Sahoo et al., 2021). Based on the availability of authenticated datasets, we selected three study locations dispersed across the Terai agro-climatic region of North Bengal, India. The three study locations are Berubari (latitude 26.38 °N, longitude 88.75 °E), Jayanti (latitude 26.70 °N, longitude 89.69 °E), and Tamaguri (latitude 26.07 °N, longitude 89.38 °E), as shown in Fig. S1 (ESRI, 2009).

Datasets

The study used five meteorological parameters as inputs to estimate $E{T}_0$ (mm/day). The five inputs were the minimum atmospheric temperature ( $T_{min}$) (°C), maximum atmospheric temperature ( $T_{max}$) (°C), wind speed ( $W_s$) (m/s; at 2 m above ground), relative humidity ( $R_h$) (%), and solar radiation ( $S_r$) (MJ m ${ }^{-2}$ d ${ }^{-1}$). The latest available meteorological datasets for three locations were obtained from National Centres for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) databases that provide benchmarked datasets (https://swat.tamu.edu/data/cfsr). The time variant data of each input parameter and that of the reference $E{T}_0$ of the three study locations are presented in Figs. S2 to S4.

No preprocessing was carried out as there were no missing values or outliers observed in the datasets. The statistical measures, such as minimum (Min), maximum (Max), standard deviation (Sd), coefficient of variation (Cv), and skewness coefficient (Sk), of each input parameter are presented in Table S1. It was observed that the skewness coefficient of wind speed ( $W_s$) has a positive value, while other parameters have both positive and negative skewness coefficients. The lowest coefficient of variation was recorded for the maximum temperature.

The acquired datasets were divided into two parts: 70% for training (from 1st January 2000 to 16th March 2010), consisting of 3,278 samples, and the rest 30% for model testing (from 17th March 2010 to 31st July 2014), having 1,598 samples. The training sets were used for model construction, while the test sets were used to evaluate the model’s performance.

Combinations of input parameters

Due to various constraints, the complete dataset of all input parameters may not always be available from all study locations. In such circumstances, limited input parameters are available instead of a complete set of five input parameters. Our model was designed to deal with various combinations of five meteorological parameters: $T_{min}$, $T_{max}$, $W_s$, $R_h$, and $S_r$. A total of eight possible combinations of five input parameters were considered to estimate $E{T}_0$ for three study locations. The eight possible input combinations ( $C1$ to $C8$) of meteorological parameters are presented in Table S2. The two input parameters $T_{min}$ and $T_{max}$ are measured using a common instrument (a thermometer) and are always available in pairs. This pair of input parameters was considered the base combination.

Reference values of $E{T}_0$

EvatCrop was designed based on the supervised learning algorithm. For proper training and testing of the model, the estimated values of $E{T}_0$ were compared with the reference (target) values of $E{T}_0$. The FAOPM56 method is recognized as the only standard method to compare the performances of any other models used for estimating $E{T}_0$ (Quinn, Parker & Rushton, 2018; Valipour, Gholami Sefidkouhi & Raeini-Sarjaz, 2020; Sarma & Bharadwaj, 2020). The reference values of $E{T}_0$ (mm/day) were calculated using the FAOPM56 method (Allen et al., 1998):

(1) $$E{T}_{0PM}=\frac{0.408\times \mathrm{\Delta }\times (R_a-G)+\gamma \times \frac{900}{T+273}\times u_2\times (e_s-e_a)}{\mathrm{\Delta }+\gamma \times (1+0.34\times u_2)}$$where the variables $R_a$, G, T, $u_2$, $e_s$, $e_a$, $\mathrm{\Delta }$, and $\gamma$ bear the same meaning and units as previously described in Banerjee, Sarkar & Ghosh (2022). All these variables were estimated using 24 additional equations involving different atmospheric inputs such as temperature, humidity, wind speed, etc., as suggested by Allen et al. (1998) to estimate $E{T}_0$.

Model architecture

Unlike a typical neuro-fuzzy system, our proposed EvatCrop model is a neuro-quasi-fuzzy system. The novelty in the architecture of this model is that an artificial neural network (ANN) is hybridized with a quasi-fuzzy inference system instead of a classical fuzzy inference system. It consists of two modules: a quasi-fuzzy module and an ANN module. The quasi-fuzzy module is coupled with an ANN module in such a way that it accepts the meteorological parameters as inputs and processes them to give the interim outputs. These interim outputs are then fed to the ANN module as inputs. After processing, the ANN module provides the estimated value of $E{T}_0$ as the final output. This model utilizes the merits of both fuzzy systems and neural networks. The model architecture and the process flow diagram are presented in Figs. S5 and S6, respectively. The model was implemented using the Scikit-learn library package in Python.

Quasi fuzzy module

A typical fuzzy inference system (FIS) consists of four functional components: a fuzzifier, a knowledgebase, an inference engine, and a defuzzifier (Odhiambo, Yoder & Yoder, 2001). The fuzzifier transforms the real numerical inputs into relevant fuzzy sets. The knowledge base contains the rule base (a group of control rules) and a database that defines the information about the membership of each fuzzy set. The inference engine fires some selected rules based on the control strategy and generates fuzzy outputs. Finally, the fuzzy output is transformed into numerical output by the defuzzyfier. Our proposed quasi-fuzzy module does not have any defuzzifier component. The defuzzifier component is replaced by a normalization component that normalizes the outputs produced by the inference engine before feeding them to the ANN module. An input interface unit is designed to accept any input combination out of the eight combinations ( $C1$ to $C8$).

Normalization unit

The quasi-fuzzy module is designed by replacing the defuzzifier of a typical FIS with a normalization unit. As the rules and their numbers are different in each rule group, the proportionate contribution of each rule should be considered. The proportional contribution is obtained by normalization. The normalization unit calculates the proportionate contribution of the output of each rule in a particular rule group. This normalized output of a rule is considered an interim output. For the $i-th$ rule with output $P_i$, the normalized output is:

(6) $$I_i=\frac{P_i}{{\sum }_{i=1}^N{P}_i}$$where N is the total number of rules in a rule group ( $R_G$).

After normalization, the set of rules in a rule group produces a set of normalized outputs. For convenience, the set of normalized outputs is called the interim output vector {I}. The interim output vector {I} is considered as input to the ANN module. The rule groups along with the number of rules in each rule group for each set of input parameters are presented in Table S3.

Performance metrics

Evaluation of the performance of a model and a comparative analysis with other contemporary models play a significant role in any model-building process. For evaluation and comparison of the performances of different models, several metrics have been suggested in the literature. Instead of a single metric, a combination of multiple metrics is more trustworthy for evaluating a model (Chai & Draxler, 2014). In this article, a well-accepted performance evaluation strategy was adopted where the performance of each model was evaluated in terms of four common metrics; coefficient of determination ( $R^2$), degree of agreement ( $d$), root mean squared error (RMSE), and root mean squared relative error (RMSRE). $R^2$ and $d$ are two goodness metrics that measure the model’s accuracy. The $R^2$ projects the relation between the actual and estimated values, and d measures the consistency between them. RMSE and RMSRE are two error metrics that project the errors in the estimated values. The mathematical equations of the evaluation metrics are defined below (Tao et al., 2018; Yaghoubi, Bannayan & Asadi, 2020).

(9) $$R^2=\frac{{\sum }_{i=1}^n{(E{T}_{0P{M}_i}-E{T}_{0_{Pi}})}^2}{{\sum }_{i=1}^n{(E{T}_{0P{M}_i}-{\overline{ET}}_{0_P})}^2}$$

(10) $$d=1-\frac{{\sum }_{i=1}^n{(E{T}_{0_{Pi}}-E{T}_{0P{M}_i})}^2}{{\sum }_{i=1}^n{(|E{T}_{0_{Pi}}-\overline{E{T}_{0PM}}|+|E{T}_{0P{M}_i}-\overline{E{T}_{0PM}}|)}^2}$$

(11) $$RMSE=\sqrt{\frac{1}{n}\sum _{i=1}^n{(E{T}_{0P{M}_i}-E{T}_{0_{Pi}})}^2}$$

(12) $$RMSRE=\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(\frac{E{T}_{0P{M}_i}-E{T}_{0_{Pi}}}{E{T}_{0P{M}_i}}\right)}^2}$$where $E{T}_{0P{M}_i}$ is the reference value of $E{T}_0$ obtained using the FAOPM56 method, $E{T}_{0_{Pi}}$ is the estimated value, $\overline{E{T}_{0PM}}$ is the mean value of reference evapotranspiration obtained using FAOPM56 method, $\overline{E{T}_{0_P}}$ is the mean predicted value of reference evapotranspiration and $n$ is the total number of input patterns.

For a comprehensive study of the performance of our proposed model in comparison to other prevailing models, two other metrics, average goodness (Ag) and average error (Ae) were used. Where Ag is the average value of two goodness metrics: coefficient of determination ( $R^2$), and degree of agreement (d), and is defined as:

(13) $$Ag=\frac{R^2+d}{2}$$

Similarly, Ae is the average value of two error metrics: root mean squared error (RMSE) and root mean squared relative error (RMSRE), and is defined as:

(14) $$Ae=\frac{RMSE+RMSRE}{2}$$

A higher value of Ag indicates a better goodness of fit, while a lower value of Ae indicates a better performance of a model.

Decision tree (DT)

A DT model uses a greedy technique to build a decision tree from root to leaf nodes. A feature that best divides the data into two or more parts is used as the root node of the tree. Thereafter, until a stopping criterion is met, the tree is developed by recursively splitting the data into smaller subsets depending on the feature values. The aim of a decision tree is to find the feature that yields the greatest information gain for each split, resulting in a tree that illustrates the relationships between characteristics and the target variable (Duda, Hart & Stork, 2006; Murty & Devi, 2015). For the present study, the DT model was developed using the scikit-learn library package in Python. The two hyperparameters: the maximum depth and the minimal cost-complexity pruning were adjusted using a fivefold grided search strategy using the training dataset. The range of maximum depth was set between 2 and 16 and the minimal cost-complexity pruning was set between 0.0 and 1.0.

Artificial neural network (ANN)

Artificial neural network (ANN) models try to mimic the information processing techniques of a biological neuron (Fausett, 2006). A classical ANN model consists of processing units called artificial neurons, which are interconnected using synaptic weights, and it learns patterns from data by adjusting these weights. Various architectures were proposed to design an ANN using a variety of activation functions. For the comparative study, the ANN model was implemented using the scikit-learn library package in Python. A two-layer ANN model with one hidden layer and one output layer was developed. The activation functions, the values of hyperparameters, and the number of hidden layer neurons were adjusted using a five-fold grided search strategy using the training dataset. The outperforming activation functions were selected by trial and error. The number of hidden layer neurons was set between 2 and 16. Three learning algorithms: the adaptive moment estimation, the L-BFGS, and the stochastic gradient descent were used to train the model, and the best one was selected for the present study.

Adaptive neuro fuzzy inference system (ANFIS)

The ANFIS model was first proposed by Jang (1993). This hybrid model combines the benefits of fuzzy logic and artificial neural networks using a three-stage architecture. In the first stage, the input variables are transformed into fuzzy sets using membership functions. In the second stage, a decision is achieved by applying a set of IF-THEN rules to the fuzzy inputs. The final stage reduces the outputs of the previous stage to a single crisp value by using defuzzification. The parameters for the rule antecedents and membership functions are weighted to reduce the difference between what happened and what was predicted. Due to the interpretability of fuzzy logic and the capacity for learning and adaptation of an ANN, the ANFIS excels at solving complicated, non-linear issues that traditional rule-base systems find challenging (Du & Swamy, 2013; Patel & Gianchandani, 2011). For the present study, the ANFIS model was developed using the ANFIS toolbox (https://in.mathworks.com/help/fuzzy/anfis.html) of MATLAB (version R2022a). The total number of fuzzy membership functions was set to three, and the parameters associated with membership functions (pre-parameters) of the developed ANFIS model were adjusted by a backpropagation algorithm on the training datasets.

Estimating $E{T}_0$ of the Berubari location

The values of the four primary metrics ( $R^2$, $d$, RMSE, and RMSRE) along with two combination metrics; average goodness ( $Ag$) and average error ( $Ae$) obtained for the four models (DT, ANN, ANFIS, and EvatCrop) at the training and testing phases for the Berubari location are summarised in Tables S5 and S6, respectively. All four models were compared in terms of $Ag$ and $Ae$ against different input combinations (from $C1$ to $C8$). The results presented in Table S5, depict that in the training phase, the EvatCrop was the most accurate model for all input combinations. The ANN was the worst model, whatever the input combinations. Table S5 depicts that for combinations $C2$, $C5$, and $C8$, ANFIS has the same accuracy compared to EvatCrop in the training phase. For other input combinations, DT and ANFIS fought neck-to-neck to secure second place. The results in Table S5 show that DT, ANN, and ANFIS were not quite as good as EvatCrop. They were about the same for input combinations $C2$, $C5$, $C6$, and $C8$ when it came to $Ag$ (differences ranging from 0.10% to 1.83%) and $Ae$ (differences ranging from 2.66% to 68.69%). The performances of all the models are very poor and unacceptable for $C1$, $C3$, $C4$, and $C7$. In conclusion, in the training phase, for five input combinations, the performances of DT, ANFIS, and the EvatCrop are comparable to each other. For the four input combinations $C6$, $C7$, and a three-input combination $C2$, the accuracy differences between all four models are also marginal, while the ANN is the worst model in all the cases.

In the testing phase, all the models were compared against different input combinations. Table S6 shows that EvatCrop did better than all the other models in terms of having the highest $Ag$ (from 0.675 for $C1$ to 0.992 for $C8$) and the lowest $Ae$ (from 0.116 for $C8$ to 0.737 for $C1$), no matter what combinations of inputs were used. It was the best performer with all input parameters (combination $C8$). However, the performances of other models vary with input combinations. In Table S6, we can see that DT, ANN, and ANFIS are not as good as EvatCrop when it comes to $Ag$ (differences range from 0.50% to 1.65%) and $Ae$ (differences range from 20.53% to 40.23%) for input combinations $C2$, $C5$, and $C6$. The DT model was the least accurate model for $C1$ (with $Ag$ = 0.650 and $Ae$ = 0.764), $C2$ (with $Ag$ = 0.971 and $Ae$ = 0.228), $C3$ (with $Ag$ = 0.672 and $Ae$ = 0.737), $C5$ (with $Ag$ = 0.967 and = 0.242), and $C6$ (with $Ag$ = 0.977 and $Ae$ = 0.201). For $C4$, ANN is the worst performing (with $Ag$ = 0.710 and $Ae$ = 0.705). Whereas, ANFIS was the least accurate model for $C7$ (with $Ag$ = 0.639 and $Ae$ = 0.814) and $C8$ (with $Ag$ = 0.658 and $Ae$ = 0.821). A possible reason is that the ANFIS model is developed using the grid partition approach, where the number of rules increases exponentially with the increasing number of input parameters, leading to increased computation time and reduced performance. A similar observation was reported by Yeom & Kwak (2018). Additionally, the ANFIS model suffers from an out-of-range problem in the test set that further degrades its performance, as described by Fu et al. (2020).

To study the impact of the meteorological inputs on model accuracy at the Berubari location, we provide a detailed analysis of the results during the testing phase using EvatCrop as the best optimum model. As reported in Table S5, using only two inputs, $T_{min}$ and $T_{max}$ ( $C1$), the EvatCrop model provided the lowest accuracy ( $Ag$ = 0.675 and $Ae$ = 0.737). Using three input variables, when wind speed $W_s$ was added with $T_{min}$ and $T_{max}$ ( $C2$), the performances of the EvatCrop model increased significantly. The value of $Ag$ was increased by about 31.00% from 0.675 to 0.981, and the value of $Ae$ was decreased by almost 75.00% from 0.737 to 0.184. When relative humidity ( $R_h$) or solar radiation ( $S_r$) was added with two inputs ( $T_{min}$ and $T_{max}$) instead of $W_s$ ( $C3$ and $C4$), the estimation accuracy of the EvatCrop model dropped significantly. The $Ag$ value was cut down by 26.60% for combination $C3$ and 36.18% for combination $C4$. Correspondingly, the $Ae$ values were increased by 73.00% for combination $C3$ and 72.00% for combination $C4$. When four parameters were used as inputs (combinations $C5$, $C6$, and $C7$), two combinations with $W_s$ ( $C5$ and $C6$) produced reasonably good results. The average value of $Ag$ for EvatCrop against two combinations, $C5$ and $C6$, was obtained as 0.988, which is almost 27.00 and 23.00% better than that of combinations $C3$ and $C4$, respectively. Correspondingly, the average value of $Ae$ for EvatCrop against two combinations $C5$ and $C6$ was observed to be 0.147, which represents a huge drop of 78.40% and 76.10% for combinations $C5$ and $C6$, respectively. However, when $W_s$ was not present in the four-input combination ( $C7$), the accuracy of all the models dropped significantly. For combination $C7$ instead of $C5$, the value of $Ag$ was cut down by 19.35% for EvatCrop. Similarly, when combination $C7$ was taken instead of $C6$, the values of $Ag$ were decreased by 20.00% for EvatCrop.

The results presented in Tables S5 and S6 indicate that better accuracy was achieved when the wind speed ( $W_s$) was available as one of the input parameters. This observation specifies that the meteorological parameter wind speed ( $W_s$) significantly contributes to estimating the $E{T}_0$ across the targeted study region. Similar observations were reported by Traore, Wang & Kerh (2010) in the Sudano-Sahelian zone of Burkina Faso, Yang et al. (2019) in China, and the literature survey conducted by Amani & Shafizadeh-Moghadam (2023). Although the complete set of five input parameters provided the highest estimation accuracy for $E{T}_0$, fairly good performance was achieved with the other three input combinations, $C2$, $C5$, and $C6$. Based on estimation accuracy, the input combinations are ranked as $C8$, $C6$, $C5$, $C2$, $C7$, $C4$, $C3$, and $C1$. The results presented in Tables S5 and S6 summarise that (i) the EvatCrop was the most accurate model for estimating $E{T}_0$ at the Berubari location; (ii) the ANFIS model performed worst in the testing phase for combinations $C7$ and $C8$ due to the out-of-range problem; and (iii) the combinations that include $W_s$ as one of the input parameters (combinations $C2$, $C5$, $C6$, and $C8$) produce more accurate results than those without $W_s$. Figure S10 presents the scatter plots of reference values $E{T}_{0PM}$ versus the values of $E{T}_0$ estimated by the EvatCrop against eight input combinations. This figure shows that the highest performance for $E{T}_0$ estimation was achieved when all the input parameters were available. Figures S11 and S12 display the boxplot and Tailor diagrams, respectively of the four models for a better comparative analysis of the results obtained in the testing phase. Additionally, the distribution plots for the estimated and reference $E{T}_0$ are provided in Fig. S13.

Estimating $E{T}_0$ of the Jayanti location

Results obtained in the training phase at the Jayanti location are reported in Table S7. From Table S7, we observe that in the training phase, the EvatCrop model comparatively outperformed all other models for input combinations $C1$, $C3$, $C4$, $C5$, and $C7$, and for the other combinations, the DT model was the best. In the case of combinations $C2$, $C6$, and $C8$, DT performed 0.20% (2.84%), 0.30% (14.43%), and 0.50% (62.03%) better than the EvatCrop in terms of $Ag$ ( $Ae$), respectively. It is to be noted that, in the training phase, the ANFIS model exhibited the same performance as EvatCrop in terms of $Ag$ for combinations $C5$, $C6$, and $C8$. For the other five combinations ( $C1$, $C2$, $C3$, $C4$, and $C7$), the ANFIS and DT models secured the second position with nearly equal performances. For the combinations $C2$, $C5$, $C6$, and $C8$, the accuracy of all the models was similar in terms of $Ag$ (difference ranging from 0.10% to 2.10%) and $Ae$ (difference ranging from 3.53% to 80.64%). In conclusion, for combination $C8$, the DT, ANFIS, and EvatCrop models produced their best results, and their performances were comparable to each other. Whereas, the ANN model performed worst.

Table S8 summarizes the results obtained in the testing phase for the Jayanti location. The results indicate that EvatCrop performed best with the highest (lowest) average values of $Ag$ ( $Ae$), irrespective of the input combinations. The accuracies of the other three models, DT, ANN, and ANFIS, are also comparable with EvatCrop for input combinations $C1$, $C4$, $C5$, and $C6$ in terms of $Ag$ (differences ranging from 0.20% to 2.05%) and $Ae$ (differences ranging from 0.79% to 53.28%). Except for input combinations $C7$ and $C8$, the DT model exhibited the poorest performance ( $Ag$ ranging from 0.731 to 0.956 and $Ae$ ranging from 0.210 to 0.519).

To reveal the impact of the meteorological inputs on model accuracy, an in-depth analysis was performed based on the results summarised in Table S8 in the testing phase. The EvatCrop exhibited the highest level of accuracy at the Jayanti location when all five inputs ( $T_{min}$, $T_{max}$, $W_s$, $R_h$, and $S_r$) were considered as inputs (combination $C8$). The input combination $C1$ with only two input variables ( $T_{min}$ and $T_{max}$) produced the worst results ( $Ag$ = 0.784 and $Ae$ = 0.504). From Table S8, it is observed that the differences in accuracy in terms of $Ag$ ( $Ae$) in comparing $C1$ to $C8$ are 31.28% (73.41%). The value of $Ag$ ( $Ae$) decreased (increased) by 12.83% (65.01%) when we removed $W_s$ from combination $C8$ (combination $C7$). But when we removed $R_h$ instead of $W_s$ (combination $C6$), the performance of EvatCrop marginally tarnished in terms of $Ag$ (decreased by 0.10%) and $Ae$ (increased by 2.19%). By removing $S_r$ from $C8$ (combination $C5$), $Ag$ was reduced by 3.76%, and $Ae$ was increased by 42.73%. It is to be noted that among the three four-input combinations, $C6$ produced the comparatively best accuracy, as observed in Table S8. For three input combinations, when we removed $W_s$ from combination $C6$ (combination $C4$), the performance of EvatCrop decreased further by 17.53% in terms of $Ag$ and increased by 62.78% with respect to $Ae$. The drop in $Ag$ and rise in $Ae$ were reasonably small when we replaced $S_r$ by $R_h$ in $C4$ (combination $C3$). The $Ag$ decreased by 2.71% and the $Ae$ increased by 7.36%. However, when $R_h$ was replaced by $W_s$ in combination $C3$ (combination $C2$), the overall performance of EvatCrop increased significantly. The $Ag$ increased by 16.56%, and the $Ae$ dropped by 82.00%. The results presented in Tables S7 and S6 summarize that (i) EvatCrop is the most accurate model for estimating $E{T}_0$ in the testing set at Jayanti location; (ii) for combinations $C2$, $C6$, and $C8$, the DT model suffered from overfitting of data in the training phase, whereas the performance degraded in the testing phase; (iii) the ANFIS model failed to provide a good estimation for combinations $C7$ and $C8$ due to an out-of-range problem; and (iv) the input combinations with $W_s$ (combinations $C2$, $C5$, $C6$, and $C8$) yielded more accurate results than the other four combinations without $W_s$.

The scatterplots of the reference values of $E{T}_0$ ( $E{T}_{0PM}$) and the predicted values of $E{T}_0$ are presented in Fig. S14 for the training and testing phases. Additionally, boxplots, Taylor diagrams. and distribution plots are provided for the testing phase in Figs. S15–S17, respectively.

Estimating ET₀ of the Tamaguri location

From the results reported in Table S9, the EvatCrop proved to be the most accurate model during the training phase for all input combinations at the Tamaguri location. The ANN model was the least accurate model, whatever the input combinations. The ANFIS model displayed the same accuracy as that of EvatCrop for input combination $C8$ in terms of $Ag$. However, ANFIS was slightly inferior to EvatCrop for the combinations $C1$ to $C7$ in terms of $Ag$ (differences ranging from 1.50% to 4.00%) and $Ae$ (differences ranging from 3.31% to 3.96%), respectively. For combinations $C1$, $C3$, $C4$, and $C7$, the accuracies of the DT, ANFIS, and EvatCrop were found to be very low. In conclusion, the performance of DT, ANFIS, and EvatCrop varied marginally from one another in the training phase for almost all the input combinations.

The results obtained in the testing phase at the Tamaguri location are reported in Table S10. From Table S10, it is observed that, regardless of the input combinations, the EvatCrop was the most accurate model in terms of the highest values of $Ag$ and the lowest values of $Ae$. The performances of the ANFIS model were very poor against the combinations $C7$ and $C8$ due to the out-of-range problem ( $Ag$ ranging from 0.548 to 0.948 and $Ae$ ranging from 0.110 to 0.558). But compared to the performance of ANFIS at Berubari and Jayanti, its performance was relatively acceptable ( $Ag$ = 0.948) at Tamaguri for combination $C8$. The main reason for this is the lesser number of outliers in the testing set, as observed from the boxplot (Fig. S19H). Contrary to the results obtained in the training phase, the ANN model performed comparatively better in the testing phase in terms of $Ag$ (difference ranging from 1.00% to 3.61%) and $Ae$ (difference ranging from 0.780% to 45.000%). The DT model ranked in third place, whereas the ANFIS was the worst.

To study the impact of the meteorological inputs on model accuracy at the Tamaguri location, the results obtained in the testing phase by using EvatCrop are summarised in Table S10. Using two inputs ( $T_{min}$ and $T_{max}$) (combination $C1$), the EvatCrop exhibited its lowest accuracy ( $Ag$ = 0.661 and $Ae$ = 0.641). When wind speed $W_s$ was considered with $T_{min}$ and $T_{max}$ (combination $C2$), the performance of EvatCrop was remarkably improved. $Ag$ was increased by around 47.66% from 0.661 to 0.976, whereas $Ae$ was declined by 72.54% from 0.641 to 0.176. However, it is to be pointed out that the accuracy of EvatCrop decreased significantly when relative humidity ( $R_h$) or solar radiation ( $S_r$) was considered instead of $W_s$, with two inputs ( $T_{min}$ and $T_{max}$) ( $C3$ and $C4$). $Ag$ values were reduced by 28.68% for combination $C3$ and 27.46% for combination $C4$. Accordingly, for combination $C3$, the $Ae$ values increased by 70.96%, and for combination $C4$, the increase was 70.71%. In the case of four input combinations (combinations $C5$, $C6$, and $C7$), two combinations with $W_s$ ( $C5$ and $C6$) yielded reasonably acceptable results. The EvatCrop scored an average $Ag$ value of 0.984 against two combinations, $C5$ and $C6$, which was 41.34 and 38.98% better than combinations $C3$ and $C4$, respectively. Accordingly, it was found that the average value of $Ae$ for EvatCrop against two combinations $C5$ and $C6$ was 0.141, indicating a significant decrease of 76.73% and 76.54%, respectively, for combinations $C5$ and $C6$. However, the accuracy of all the models was significantly decreased when $W_s$ was not included in the four-input combination ( $C7$). The value of $Ag$ was reduced by 23.82% when combination $C7$ was used instead of combination $C5$. Comparably, the value of $Ag$ decreased by 24.75% when combination $C7$ was used in place of $C6$. From the results reported in Tables S9 and S10, it is concluded that (i) out of the four models, the EvatCrop turned out to be the most accurate model, (ii) the out-of-range problem prevents the ANFIS model from providing an acceptable accuracy for combinations $C7$, and (iii) input combinations with $W_s$ (combinations $C2$, $C5$, $C6$, and $C8$) yielded more accurate results than the other four combinations without $W_s$.

Figure S18 provides the scatterplots for estimated values against the reference values of $E{T}_0$ in the training and testing phases at the Tamaguri location. For a more intensive analysis of the results, the boxplots, Taylor diagrams, and distribution plots for the testing phase are also presented in Figs. S19–S21, respectively.