Comparative assessment of empirical and hybrid machine learning models for estimating daily reference evapotranspiration in sub-humid and semi-arid climates

Reference evapotranspiration (RET) is a crucial hydrological cycle element, responsible for a significant portion of water loss from continental surfaces^1,2,3. It accounts for approximately 62% of the rainfall contribution, equivalent to around 73,000 km³ per year⁴. RET serves as a powerful indicator for climate change studies^5,6,7,8 and plays a vital role in many fields, including hydrology, agriculture, ecology, and water resource management. Notably, RET can also be instrumental in addressing natural hazards like dry spells, heat waves, and flash droughts^9,10,11. While direct estimation methods like latent heat flux or real evapotranspiration from climate models offer accuracy, they often require complex modelling and high-resolution data, posing practical challenges. The intricate nature of modelling the soil-vegetation-atmosphere interaction further complicates accurate RET estimation. Therefore, understanding the processes and assessing RET is crucial for effective water resource management and planning, especially in semi-arid and dry locations where water availability is restricted. Raza et al.¹² presented potential evapotranspiration (PET) and RET differentiation and categorized their empirical equations based on different meteorological factors.

The Penman-Monteith technique (PM-FAO56), a revised version of the Penman Equation ¹³, is widely considered as the most accurate method for estimating RET, and has been endorsed by FAO¹⁴ and the Task Committee on Standardization¹⁵. Despite its global acceptance, PM-FAO56 relies dependent on meteorological parameters such as wind speed, relative humidity, and solar radiation, which may be unavailable in certain weather stations, particularly in developing countries. Further, this limitation hinders its application in regions with limited weather data. Consequently, researchers have sought alternative methods for estimating RET through comparison of PM-FAO56 method with other empirical methods or methods development through meteorological data obtained from remote sensing.

Numerous researches have been conducted to investigate the effectiveness of different RET methods across different regions and climatic conditions. These studies compare empirical models with PM-FAO56 using different approaches^16,17,18,19, including (a) temperature-based methods, (b) radiation-based methods, (c) mass transfer based methods, (d) methods combining radiation and temperature, and (e) methods integrating radiation, temperature, mass transfer, and other variables. For example, Almorox et al.¹⁸ assessed eleven temperature-based potential evapotranspiration (PET) estimation methods and determined that the Hargreaves and Samani model exhibited the most accurate performance on a global scale across diverse climatic regions. Additionally, comparisons have been made between empirical models and data obtained from lysimeters^20,21. In Morocco, few scientists have investigated RET performance with existing empirical models/methods^16,20,22,23. Er-raki et al.¹⁶ evaluated three empirical RET estimation methods for Tensift Basin (Morocco’s center) and Yaqui Valley (Northwest Mexico) during 2003–2004. In a semi-arid region, they suggest using Hargreaves and Samani model without calibration (as long as the wind remains low). They suggested that calibration is required for both Priestley-Taylor and Makkink parameters, particularly for dry periods. Several methods were examined by Bouhlassa and Paré²² to choose an appropriate solution to PM-FAO56 equation for 1989–2001 in Tafilalet, an arid region in southeastern Morocco. Their findings indicate that Jensen-Haise and Thornthwaite methods best-reflected evapotranspiration obtained by Penman-Monteith-FAO method. Similary, Hadria et al.²³ conducted a calibration and validation analysis of five temperature-based empirical models in 22 meteorological stations across Morocco. Their results demonstrated that Dorji’s estimate outperformed the other empirical models and they introduced a new fit version called RET-Hadria, specifically designed for assessing RET in arid and semi-arid areas. Zeggaf²⁰ compared lysimeter results with various empirical methods for Ouled Gnaou (Morocco’s semi-arid central region) during 1975, 1977, and 1978 years. They found that Priestley-Taylor method gave better results, followed by Penman-Monteith method. Besides, researchers like Liou and Kar²⁴ and Elfarkh et al.²⁵ have used remote sensing (RS) images and processes in geographical information system (GIS) to enhance evapotranspiration estimation.

In recent years, machine learning models have garnered considerable attention in various fields^{26,27,28,29,30,31,32}. For instance, a model’s capacity to represent intricate nonlinear relationships has a significant impact on RET estimation. Goyal et al.³³ highlighted the promising findings of ML models in various climates and environments, emphasizing their ability to improve accuracy above standard empirical models. Besides, researchers have applied various ML models, such as artificial neural networks (ANN)^34,35, support vector regression (SVR)^36,37, M5 model tree^38,39, random forests (RF)^40,41,42,43, reduced error pruning tree (REPTree)^44,45, extreme gradient boosting (XGBoost)^46,47,48, light gradient boosting machine (LightGBM)^34,49,50 and decision trees (DT)^51,52,53 to estimate daily RET uing restricted meteorological data. For instance, Granata³⁸ conducted a comparative investigation of M5P Regression Tree, Bagging, RF and SVR with differing input combination (T_mean, RH_mean, R_s, U₂) for RET estimation in a humid subtropical climate region of Central Florida. They concluded that the M5P models exhibited good performance, while RF proved to be the least accurate. In China, Fan et al.⁴⁰ examined limited meteorological data to investigate four empirical models and three ML models (LightGBM, M5Tree, and RF) to estimate daily RET. Their findings suggested that LightGBM outperformed the other models, with input combinations comprising T_max, T_min, U₂, R_s, and RH_mean. Similary, Fan et al.⁴⁶ compared six ML models, including SVM, gradient boosting decision tree (GBDT), M5Tree, XGBoost, ELM and RF, and using meteorological data from eight Chinese stations. Their findings revealed that the GBDT and XGBoost models exhibited performance on par with the SVM and ELM models, while offering advantages in terms of simplicity, accuracy, stability, and reduced computational costs, making them recommended options for daily RET estimation. Additionally, Yong et al.³⁴ evaluated performance of LightGBM, ANN and decision forest regression (DFR) in five Malaysian meteorological stations and reported that LightGBM and ANN have proven stable and accurate in determining daily RET. It is worth mentioning that the ML models’ performance in daily RET estimation is influenced by various factors, including the selection of input climatic variables, model structure, basic parameters, and performance criteria. Careful consideration of these factors and the correlation between individual input variables and RET is crucial for optimizing the ML models’ accuracy and efficiency. Additionally, effective tuning of ML model parameters further enhances their performance and efficiency for correct estimation^{33,34,35,36,38,40,44,46}.

Nowadays, recent studies^33,54 in RET estimation have emphasized the use of hybrid models, combining multiple ML algorithms through blending or stacking techniques, to address the challenges posed by highly complex meteorological data. Goyal et al.³³ noted that standalone ML models may not achieve sufficient accuracy in such cases. For example, Elbeltagi et al.⁵⁴ assessed five hybrid models (additive regression (AR) with bagging, M5tree, random subspace, REPTree, ANN, and RF) for a semi-arid area in Pakistan. They concluded that AR-M5tree model is the most appropriate hybrid model for estimating RET. Hence, this highlights the effectiveness of hybrid models in improving performance while maintaining interpretability. However, their RET estimation’s utilization is currently limited, and the available information on this subject is incomplete and fragmented. There is a need for further research and investigation to fully explore the hybrid models’ potential and effectiveness in addressing evapotranspiration estimation challenges.

In Morocco, few studies^51,55,56 have focused on evaluating ML models for estimating RET. Recently, Lachgar et al.⁵¹ studied the performance of five ML models, like RF, Linear Regression (LR), SVR, k-Nearest Neighbor (k-NN), and DT for estimating RET between 2011 and 2019 in Fez. They highlighted the ML models’ ability to capture the variance in RET. In Marrakesh, El Hachimi et al.⁵⁵ investigated the performance of SVM, RF, DT, k-NN, LR, XGboost for estimating RET during 2013–2020 and found that XGboost surpasses the other models, followed by RF. To our best knowledge, there is currently no existing comparative research on the use of hybrid models for estimating daily RET in Morocco. Thus, the novelty of this research lies in the comprehensive comparison of empirical models, ML models, and their hybrid models for estimating daily RET in subhumid and semi-arid climates.

The primary objectives of the present research are as follows: (i) estimating RET using eight empirical models, namely, Valiantzas 2013 (VAL2013a and VAL2013b); Dalton 1802 (Dal1802); Trabert 1896 (Trab1896) Hargreaves, 1975 (Harg1975); Irmak and Haman, 2003 (Irs2003); Hargreaves and Samani, 1985 (HargS1985); and Allen and Pruitt, 1986 (BC1986) and comparing their performance with standard FAO-PM56 method, (ii) development of ML models (RF, M5P, XGBoost and LightGBM) and their hybrid models (RF-M5P, RF-XGBoost, RF-LightGBM and XGBoost-LightGBM) using different meteorological input combinations (based on maximum and minimum air temperature (T_max and T_min), relative humidity (RH_mean), solar radiation (R_s), and wind speed (U₂)), (iii) evaluating performance of developed ML and hybrid models using different statistical indices to determine the best for RET estimation. The findings of this research will contribute to assessing the performance and suitability of selected machine learning models for reference evapotranspiration (RET) estimation under the specific climatic conditions of Morocco. This will support improved planning for water resource management and irrigation practices.

Study area and data collection

Figure 1 illustrates the research region’s geographic distribution. This study region consists of two perimeters which are among Morocco’s most important irrigated perimeters. It covers an area of 6,007.14 km², which is 57.2% for Gharb and 42.8% for Loukkos perimeter. It has a Mediterranean climate (Csa) with an oceanic impact, according to Köppen’s classification. The differences between the two perimeters are shown in Table 1.

Observed daily minimum and maximum daily air temperature (T_min and T_max), minimum and maximum relative humidity (RH_max and RH_min), solar radiation (R_s), mean wind speed at 2 m height (U₂) and precipitation (P) were obtained from five weather stations to estimate reference evapotranspiration. These stations were strategically chosen, with three situated within the Gharb perimeter, representing three out of the five districts, and the remaining two located within the Loukkos perimeter (Figs. 1 and 2). The selection criteria considered the availability of reliable data and practical constraints, allowing for a focused analysis on representative stations within the study area. Meteorological data were provided by two Regional Offices for Agricultural Development, ORMVAG and ORMVAL, with varying data collection periods covering both perimeters. Detailed information on the studied station locations, statistical characteristics of observed meteorological data, and missing data percentage are described in Table 2.

Geographical ___location of study area, meteorological stations distribution and digital elevation model extracted from Shuttle Radar Topography Mission (SRTM) ) (created by ArcGIS 10.8.2).

Full size image

Digital elevation model extracted from Shuttle Radar Topography Mission (SRTM) (created by ArcGIS 10.8.2).

Full size image

Missing data imputation

In general, the database contains missing data due to weather station malfunctions. To address this issue, we employed either deletion or imputation. For the imputation process, we utilized the Multivariate Chain Equations Imputation Method (MICE) developed by Van Buuren and Groothuis-Oudshoorn⁵⁷. This approach enabled us to assess whether imputing missing data would impact the selection of estimation methods. Otherwise, the other missing data was completely removed. The MICE method comprises three main steps. Firstly, a regression model is chosen for the variable being studied. Then, missing data values are iteratively assigned random values based on observed data. Finally, imputed values are estimated using the regression coefficients obtained for each dataset. We opted for this method due to its practicality and effectiveness in analysing precipitation data⁵⁸.

Table 2 reveals that the solar radiation series for the Sidi Allal Tazi station had the highest proportion of missing data, accounting for 10.77% of the observed series. Through the MICE method, we imputed 9.4% of the solar radiation series for the Sidi Allal Tazi station. Overall, the deleted data represented between 1.37% and 4.48% of the database for each station studied.

Estimating evapotranspiration via FAO-56 Penman-Monteith and empirical models

A series of models were designed by researchers to estimate reference evapotranspiration^59,60,61,62. In this research, eight empirical models, divided into four groups (combination, mass transfer, radiation, and temperature), were selected, as indicated in Table 3. Model selection focused on variables’ availability, usage extent, and simplicity. Subsequently, we compared these models with Penman-Monteith model [PM-FAO56, Eq. 1¹⁴]. Table 4 summarizes the climate parameters for each empirical model.

Machine learning models and hybrid models description

Different equations and models to estimate daily RET were compared. Linear Regression (LR) was used to compare eight empirical equations with the PMFAO56 model. Additionally, we explored the impact of various meteorological variables on RET estimation using ML algorithms like Random Forest (RF), M5 Pruned (M5P), eXtreme Gradient Boosting (XGBoost), and Light Gradient Boosting Machine (LightGBM). Additionally, we combined these models to create hybrid models like RF-M5P, RF-XGBoost, RF-LightGBM and XGBoost-LightGBM.

Linear regression (LR)

Linear regression is a well-known approach for modelling a dependent variable’s value via one or more independent variables. In this research, The LR equation is written as follows.

\(\:{\text{R}\text{E}\text{T}\:}_{\text{P}{\text{M}}_{\text{i}}}\) represents observed values (estimated by PM-FAO 56); \(\:{\text{R}\text{E}\text{T}\:}_{\text{c}\text{a}{\text{l}}_{\text{i}}}\) represents values estimated by different empirical models.

Random forests (RF)

RF model, introduced by Breiman⁶⁶, is an ensemble approach that integrates numerous decision trees to create a powerful prediction model. It is widely utilized for regression and forecasting tasks due to its ability to capture complex, non-linear interactions between features. It works by generating a collection of random binary trees through bootstrapping, where each tree is trained on a randomly sampled subset of observations from the training dataset. The remaining data, known to as “out-of-bag” (OOB) data, is used for evaluating the model’s performance. RF exhibits several advantages, including strong generalization capabilities, robustness to outliers, and the ability to tune hyperparameters easily. By aggregating the results from individual trees, RF produces a final prediction, often using methods like majority voting. This ensemble approach helps mitigate overfitting and minimize variance by training on various data samples. To further control overfitting, the minimum leaf size parameter can be adjusted, requiring a minimum number of observations to generate child nodes.

M5 pruned (M5P)

The M5P model, also known as M5 Pruned, is a decision tree algorithm introduced by Quinlan⁶⁷. This model runs in two steps, providing a novel approach to regression challenges. It separates the input data into subgroups in the first phase and applies linear regression models to each subset based on their partial attribute values. This enables the model to record variable relationships and construct regression equations at each node. Furthermore, the M5P model is typically built as a tree, starting with a root node and branching out into subnodes that reflect the regression equations. It can easily handle huge datasets with high dimensionality⁶⁸, making it useful for investigating complicated systems like evapotranspiration estimates. Additionally, it does not require parameter adjustment, which simplifies the process.

Extreme gradient boosting (XGBoost)

XGBoost model is an improved Gradient Boosting Machines (GBMs) version presented by Chen and Guestrin⁶⁹ that expands on the notion of “boosting” weak learners. Through additive training procedures, it combines numerous weak models to generate a powerful learner. By simplifying goal functions and providing parallel calculations during training, XGBoost tries to minimize overfitting while decreasing computational costs. It offers a scalable and effective solution for both regression and classification workloads²⁷. With features such as distributed computing, pruning strategies, and management of missing data, the approach is meant to maximize speed. Because of its efficacy, versatility, and capacity to handle big datasets, XGBoost has become a popular choice.

Light gradient boosting machine (LightGBM)

LightGBM model is a gradient-boosting architecture that improves model performance while using less memory than conventional models. LightGBM is distinguished by its novel leaf-wise development, which develops trees by focusing on individual leaves rather than sequentially expanding the branches⁷⁰. This approach enables efficient tree formation and increased computing performance. LightGBM also incorporates two approaches: gradient-based one-sided sampling and exclusive feature bundling (EFB). These approaches improve model performance by allowing for more efficient feature sampling and grouping. Overall, LightGBM is an efficient technique for dealing with massive datasets and producing accurate predictions while conserving resources. One of the most essential variables impacting the accuracy of a given model is the selection of proper hyper-parameters. The hyper-parameters employed in this study, shown in Table 5, were chosen through grid search optimization, supported by prior ___domain knowledge. This approach enabled systematic evaluation of parameter values to minimize prediction error and optimize performance. The resulting tuning achieved a balance between computational efficiency and accuracy across all models.

Weighted hybridization for ML algorithms

Weighted hybridization is an approach that combines different algorithms to increase the accuracy of evapotranspiration estimate. Individual algorithms are given varied weights based on their performance, as proposed by Nourani et al.³⁶. More information can be found in^33,71. The hybrid model provides more accurate and resilient predictions by harnessing the capabilities of each algorithm and assessing their relative relevance. This method is critical for reducing the constraints of individual algorithms and increasing the dependability of evapotranspiration estimations.

To ensure consistent scale and improve modelling capabilities, the input data for the ML models were normalised. This process, described by Eq. (11), transformed the data to a range between 0 and 1.

where, x_norm represents the normalised value, x₀ is the real value, and x_min and x_max are the minimum and maximum values respectively.

It should be noted that the XGBoost model does not require input variable normalization since it is not sensitive to monotonic input variable normalization.

Input combinations

To investigate the impact of different meteorological variables on RET estimation, six combinations were utilized, as outlined in Table 6. The inputs for the ML models included air temperature (T_max and T_min), relative humidity (RH_mean), solar radiation (R_s), and wind speed (U₂). Likewise, the observed meteorological data were divided into two sets: a training set comprising 70% of the data and a separate testing set for evaluating the performance of the models. This division ensured that the models thoroughly evaluated on a substantial amount of data from each ___location and were rigorously assessed on an independent dataset. All the chosen ML models and simulations were implemented using R software (version 4.2.2).

Evaluation performance

To evaluate the model’s performance in comparison to the standard Penman-Monteith model (FAO-PM56), four statistical metric parameters were employed. These metrics encompassed the Kling Gupta Efficiency index (KGE, Eq. 12), Coefficient of determination (R², Eq. 13), Mean Squared Error (RMSE, Eq. 14), and Root relative squared error (RRSE, Eq. 15). The selection of these metrics aimed to assess the precision, accuracy, under/overestimation, and provide a means for model comparison^72,73,74. To rank the models, the RMSE and RRSE values were arranged in ascending order, while the R² and KGE values were ordered in descending order. The specific formulas for these metrics can be found in Table 7.

Proposed model development for RET estimation

Figure 3 illustrates the flowchart outlining the suggested empirical, machine learning (ML), and hybrid models for estimating reference evapotranspiration (RET). LR, RF, M5P, XGBoost, LightGBM, RF-M5P, RF-XGBoost, RF-LightGBM, and XGBoost-LightGBM are the models used for RET estimation. As input variables, these models use six distinct climatic combinations (Comb1 - Comb6). The performance of each model was extensively explored by assessing the statistical metrics parameters listed above.

Proposed empirical, ML and hybrid models flowchart for RET estimation.

Full size image

Correlation between PM-FAO 56 daily RET and meteorological variables

Figure 4 shows a correlation between meteorological variables and daily reference evapotranspiration estimated by PM-FAO 56 model at each station studied. Results shows the RET is primarily affected by solar radiation where Pearson’s coefficients were above 78%, indicating a good correlation. This suggests that models that incorporate radiation-based variables may perform better in estimating RET compared to those that rely on temperature-based or mass transfer-based variables. For all Gharb stations, R_n obtained the highest correlation, with values ranging from 0.96 to 0.97, followed by R_s and R_a. However, R_s is the most correlated to R_n and R_a for all Loukkos stations. One of these three variables is explicitly contained in all combination models and those based on radiation-based (Table 3). Temperatures occupy the second place as r values vary from 0.53 to 0.79. Moreover, maximum temperature for Loukkos stations is more correlated with reference evapotranspiration than with mean and minimal temperatures. This leads to saying that a model containing maximum temperature coupled with solar radiation could better estimate RET for these stations. Recent research by Chia et al.⁷¹ pointed out that temperature and radiation are indispensable for estimating RET in semi-arid regions. Conversely, in sub-humid regions, the RET estimation requires the inclusion of evaporation in addition to temperature and radiation. This review underscores the significance of considering distinct climatic conditions when estimating RET in various regions.

Figure 4 further indicate that air vapor pressure deficit (VPD) is moderately correlated with reference evapotranspiration ranging from 0.56 to 0.74. Moreover, mass transfer models usually use this VPD variable as input and lack other variables that could improve RET estimation. On the other hand, relative humidity is negatively correlated with reference evapotranspiration, with r values varying from − 0.63 to -0.27. In the Gharb stations, it is noteworthy that RH_max correlation coefficient is higher than RH_min, regardless of the study period. In line with previous studies^33,44,71, our findings support the positive correlation of T_mean, T_max, T_min, and R_s with RET. Wind speed shows a slight correlation, while relative humidity (RH_mean) exhibits a negative correlation with RET. The correlation coefficients for wind speed (U₂) are relatively low, ranging from 0.21 to 0.40, with the exception of the MB station, which exhibits a correlation coefficient of 0.50. The consideration of wind impact in RET estimation varies among researchers, with some arguing that wind is a significant factor due to potential data inaccuracies⁷⁵, while others suggest that wind has minimal influence⁷⁶ except in areas with high wind conditions. These results align with the findings of other studies^46,54, providing further evidence of the relationship between meteorological variables and RET.

Correlation of Pearson coefficients between daily RET calculated by PM-FAO 56 model and meteorological variables for different weather stations.

Full size image

Empirical models’ comparison for daily reference evapotranspiration estimates

The statistical results of the eight empirical models (Dal1802, Trab1896, Harg1975, Irs2003, HargS1985, BC1986, VAL2013a and VAL2013b) for estimating daily RET at the five meteorological stations within the Gharb and Loukkos perimeters are presented in Tables 8 and 9. As above mentioned, the computed statistical indicator values (KGE, R², RMSE and RRSE) were obtained using model performance equations [Eqs. 12–15], evaluated against the PM - FAO 56 model during the training and testing phases. The models were prioritized based on their statistical errors, and the best model was identified accordingly. KGE, R², RMSE, and RRSE were determined to be 0.291–0.974, 0.249–0.964, 0.182–1.287 mm/day and 8.572–51.913%, respectively, during training; and 0.386–0.972, 0.366–0.966, 0.172–1.302 mm/day and 8.137–47.993% during testing. As seen in the table, the empirical models exhibited minimal differences (with an RRSE gap from − 1.57 to 4.63%) between the training and testing phases.

Notably, the combination models (VAL2013a, VAL2013b) outperformed other models across all stations, with KGE R², RMSE, and RRSE ranging 0.947–0.974, 0.942–0.966, 0.318–0.402, 8.137–10.739 respectively, during testing phase. Except for Mensara station, Trab1896 performed successfully at the training and testing phases. It was noticed that The VAL13b model performs noticeably better than the VAL13a model, owing to the distinct variable requirements of each model. This is because it includes all contributing factors affecting RET.

Notably, VAL13b incorporates additional variables, such as relative humidity and wind speed, which account for the aerodynamic effects on RET. Similary, Kisi⁷⁷ evaluated seven empirical models to the PM-FAO 56 in Mediterranean climate and Valiantzas (2013b) was found to be the best model.

The analysis indicates that temperature-based models, specifically Hargreaves and Samani⁶² and Brutsaert and Chen⁷⁸, demonstrate superior performance compared to radiation-based models such as Hargreaves⁶¹ and Irmak et al.⁷⁹, as well as mass transfer models like Dalton⁵⁹ and Trabert⁶⁰. These findings are consistent with the conclusions reported by multiple researchers in the field¹⁷. Conversely, it is important to note that some researchers^20,76,80 have suggested that radiation-based models may perform better than temperature-based models. The higher effectiveness of temperature-based models can be related to temperature’s relative stability compared to solar radiation, which changes depending on conditions like as cloud cover, meteorological conditions, and time of day. Consequently, variations in solar radiation create uncertainty in radiation-based models. Generally, temperature appears to be a more powerful element than solar radiation in promoting evapotranspiration in dry or semi-arid regions where water supply is limited^72,73,74.

Specifically, the HargS1985 model demonstrates superior performance compared to the BC1986 model across all stations. This finding is congruent with Er-raki et al.¹⁶ results, who found that the HargS1985 model provides more accurate estimations in semi-arid conditions of the Tensift basin. Similarly, Almorox et al.¹⁸ reported that the HargS1985 model had the greatest overall performance after analyzing eleven temperature-based models on a worldwide scale. In contrast, other studies have reported that the Hargreaves-Samani model⁶² tends to underestimate reference evapotranspiration in arid regions and overestimate it in wetland environments^81,82.

In term of mass transfer models, Trab1896 performed better than Dal1802 at most stations. The mass transfer models used in this study were unsuitable for estimating daily RET in this region due to substantial statistical errors, with the exception of the Menasra station. In term of radiation models, Harg1975 was generally superior to Irs2003. On the other hand, we found that RET ranking is generally like SAT station data and data supplemented by Mice imputation method (SAT mice). This means that the filling of solar radiation (R_s) gaps did not affect the choice of RET estimation method.

The scatter plots in Fig. 5 illustrate the comparison between the estimated RET values obtained from the best empirical models and the FAO56-PM values during testing at the two meteorological stations. When the data points in the scatter plot were closely matched with the diagonal 1:1 trend line, a strong fit was noticed, showing high agreement between RET estimated by PM-FAO 56 and by empirical model. When the data points deviated considerably from the trend line, it indicated a poor fit, suggesting a lack of connectivity between the two previously estimated models. Overall, the data points in the plots demonstrate a strong correlation, aligning closely to the 1:1 line for models Val2013a, Val2013b, Harg195, and Harg1985.

Comparison of standalone and hybrid ML models using various input combinations

Table 10 presents the averaged statistical performance indicators (KGE, R², RMSE, and RRSE) values for estimating RET across five meteorological stations in the Loukkos and Gharb perimeters, categorized by model and input combination used during testing phase. Additionally, Fig. 6 displays the KGE and RRSE values for each meteorological station during training and testing phases, comparing the four ML models and four hybrid models with six different input combinations. Overall, the statistical indicators values were found to varied substantially based on input combination, model types, and phase employed.

It can be seen from Table 10 that, independent of the perimeter and input combination, KGE, R², RMSE, and RRSE values varied from 0.557 to 0.982, 0.585–0.979, 0.015–0.108, and 6.925–35.360 respectively during testing phase. These results clearly surpass those produced by the empirical model. The RMSE values achieved by these models were smaller than those reported by other researchers using different ML models in various regions^34,40,46,51. It’s worth noting that discrepancies in RMSE values between studies might be caused by various factors such as ML model type and input variable selection, time periods chosen, climate conditions, and data quality^83,84,85.

Scatter plots depicting the relationship between RET values estimated by different empirical models and FAO-56 PM values during testing phase at Sidi Allal Tazi and Aouamra stations.

Full size image

Among the different input combinations evaluated across all eight models, it was observed that the ML and hybrid models displayed the poorest performance when utilizing T_max, T_min, and RH_mean as input (combination 2). On average, the R² values ranged from 0.585 to 0.669, while the RRSE values varied between 26.724% and 35.360%, indicating relatively lower accuracy compared to other input combinations. This can be attributed to the negative correlation observed between RH_mean and RET estimated, as depicted in Fig. 4. Furthermore, the limited information provided by these variables may not fully capture the complex relationships involved in accurately estimating RET. When comparing combination 1 (T_max, T_min, R_s) with combination 3 (T_max, T_min, R_s, U₂), it was observed that the inclusion of U₂ improve slightly the R² (difference < 0.12) and reduced RMSE values (difference < 0.01). However, the most significant difference was found in the RRSE values, particularly in the Gharb stations, where the improvements ranged from 6.597 to 9.862%.

In the Loukkos stations, the differences in RRSE values were between 0.910% and 1.886%. These disparities suggest that the Loukkos stations experience higher wind speeds compared to those in the Gharb (Table 2), and the correlation between U₂ and RET is slightly weaker in the Loukkos (as shown in Fig. 4). Similarly, when comparing combination 1 (T_max, T_min, R_s) with combination 5 (T_max, T_min, R_s, RH_mean), improvements were observed, although they were lower than in the previous case. The RRSE differences were 3.871–8.923% for Gharb and 0.499–1.537% for Loukkos.

Goyal et al.³³ suggested that incorporating either U₂, RH_mean, or both can improve model performance. Although relative humidity is considered the least significant parameter, its addition to the combination of T_max, T_min, and R_s results in a decrease in RMSE values. Besides, the combination 6, which included T_max, T_min, RH_mean, R_s, and U₂ as input meteorological variables, demonstrated the best performance with R² and RRSE values ranging from 0.955 to 0.979 and from 6.925 to 11.272%, respectively. This improved performance can be attributed to the inclusion of additional variables capture the complex interactions and dynamics involved in estimating RET accurately.

Overall, it’s worth noting that temperature data is a foundational requirement for the models presented in Tables 9 and 10. Without sufficient temperature data, model predictions become unreliable, thereby limiting accuracy of RET estimations. Alternative gridded datasets, such as those from reanalysis products like ERA5-Land, can be used effectively to run models like Penman-Monteith FAO-56⁸⁶. Nouri et al.⁸⁶ demonstrated that the ERA5-Land dataset provides reliable RET estimates, especially in data-limited and windy regions. Their findings revealed that while some models, like recalibrated Hargreaves-Samani and Penman-Monteith with localized wind speed, performed well, others struggled due to wind speed variation.

In term of standalone ML models, the testing phase revealed the following ranking for Gharb: LightGBM6 > XGBoost6 > RF6 > LightGBM3, with average RMSE of 0.015–0.017 mm/day. For Loukkos, the ranking was: LightGBM6 > XGBoost6 > LightGBM3 > XGBoost3 with average RMSE of 0.025–0.027 mm/day. These findings support previous studies by Fan et al.⁴⁰ and Yong et al.³⁴, where LightGBM consistently outperformed other standalone ML models with an RMSE of 0.08–0.58 mm/day and 0.041–0.315 mm/day, respectively. Further, there was a minor difference between LightGBM6 model and XGBoost6 model, with LightGBM6 having a little higher RRSE value (0.34–0.5%). It should be pointed out that while XGBoost had the highest KGE value among all models, our ranking also considered the low RMSE and RSSE values, which placed XGBoost in the second position.

KGE and RRSE results for each meteorological station throughout training and testing.

Full size image

In both study area, the worst performance were obtained by M5P2 < XGBoost2 < RF2 < LightGBM2, where on average KGE < 0.772, R² < 0.661, RMSE > 0.067 mm/day and RRSE > 27.097% (Table 10). These findings contradict those of Granata³⁸, who claimed that the M5P models performed well while the RF models were the least accurate. From Fig. 6, it can be observed that RF demonstrated strong performance during the training phase across all stations, likely due to its ensemble nature and robustness to noise.

However, LightGBM’s performance was often inferior to that of the RF and XGBoost models. This could be due to its different optimization approach and less effective capturing of complex relationships. This finding aligns with the study by Chia et al.⁷¹, which reported that LightGBM exhibited relatively weaker performance than RF and the M5 tree model during the training phase. The authors explained that LightGBM’s leaf-wise optimization required sufficient training data for effective performance⁷¹. Interestingly, the situation reversed during the testing phase, demonstrating that LightGBM performed better when properly trained. It is noteworthy to mention that for all ML model, the model’s performance difference was triggered by a difference in the training and testing datasets, such as temporal differences in meteorological data patterns during both phases.

In term of hybrid ML models, the Gharb’s performance ranking was XGBoost-LightGBM6 > RF-LightGBM6 > RF-XGBoost6 > RF-M5P6. The R² (RRSE) values were 0.979 (7.598), 0.978 (7.813), 0.974 (8.685) and 0.973 (9.211) respectively. For Loukkos, the ranking models were XGBoost-LightGBM6 > RF-LightGBM6 > RF-XGBoost6 > RF-LightGBM3, with R² (RRSE) values of 0.975 (6.925), 0.974 (7.038), 0.971 (7.619) and 0.969 (7.635) respectively. From Table 10,, the result shows that XGBoost-LightGBM perform the best. This might be because XGBoost is recognized for its regularization approaches and successful handling of complicated relationships, whereas LightGBM excels in efficient computation and handling huge datasets. Consequently, their hybridization could detect a wider range of patterns and improve overall performance. In contrast, the lowest performance occurred in RF-M5P2 < RF-LightGBM2 < XGBoost-LightGBM2 < RF-XGBoost2. Nonetheless, during training phase, RF-XGBoost model gave good performance, as shown in Fig. 6.

When comparing standoalone ML models with hybrid ML models, it was found that XGBoost-LightGBM was highly close to LightGBM in term of all statistical performance indicator. For instance, compared with XGBoost-LightGBM6, RRSE of LightGBM6, XGBoost6 and RF6 was in difference of 0.318% (0.096%), 0.655% (0.592%), 2.468% (1.258%), respectively for Gharb (Loukkos). This suggests that while the hybridization approach slightly improved model performance, the improvement was not significant. It’s worth noting that factors such as computational efficiency and implementation flexibility might impact the selection between standalone ML models and hybrid models. Therefore, LightGBM and XGBoost are the recommended ML model for estimating RET in our research study.

Collectively, these studies highlight the effectiveness of standalone and hybrid machine learning models for improved accuracy in RET estimation^{33,34,38,40,51,54,55,56,71}. Similary, our findings indicate that standalone ML models exhibited better performance and accuracy compared to the empirical models for estimating daily RET at Gharb and Loukkos stations, using T_max, T_min, RH_mean, R_s, U₂ as input variables.

This research has some notable limitations that need to be addressed. The study’s period is relatively limited, spanning from 2011 to 2017. Hence, expanding the study’s period would provide a more comprehensive understanding of the ML models’ performance across diverse climatic fluctuations over a longer period. Moreover, the investigation of current study limited from subhumid to semi-arid climatic conditions only, further investigation should incorporate more climatic conditions to examine climatic variability thoroughly. Additionally, the accuracy and performance of ML and hybrid models heavily relied on the availability and quality of meteorological data, which could impact their effectiveness in areas with limited or incomplete data. To address this issue, future studies could utilise reanalysis data products, such as ERA5 ERA5-Land, and MERRA-2, which provide continuous, high-resolution data across extended temporal and spatial scales. Such data sources could prove invaluable in regions where in situ measurements are scarce. Moreover, ML models lack physical mechanisms, making it challenging to comprehend their inner workings and create accurate models without knowledge of functional specifications⁷⁴. The issue of over-fitting and under-fitting during the training/testing phases of ML models, due to the dataset random division, may also affect model accuracy. Employing advanced validation techniques could help enhance model reliability and generalizability.

Estimating reference evapotranspiration (RET) accurately has remained a major focus across a wide range of applications, including water resource management, agricultural water requirements, irrigation scheduling, and climate change assessments. In this research, we investigated the ability of four machine learning (ML) models, and their hybrid models along with eight empirical models (grouped in mass transfer-based, temperature-based, radiation-based and combination models) to estimate daily RET in subhumid and semi-arid irrigated perimeters in Morocco. For six input combinations, RF, M5P, XGBoost, LightGBM, RF-M5P, RF-XGBoost, RF-LightGBM, and XGBoost-LightGBM were all thoroughly evaluated. The results showed that combination models (particularly, Valiantzas 2013 (VAL2013b)) were the best empirical models and that temperature-based models generally outperformed radiation-based models. Compared with empirical models, ML models gave more accurate RET estimation, and the hybrid XGBoost-LightGBM models provided the highest statistical indicator values (KGE, R², RMSE, RRSE). Interestingly, the standalone ML model LightGBM also demonstrated acceptable accuracy across all stations and input combinations, indicating its potential as a promising model for RET estimation with limited data. Moreover, the XGBoost model is also an intriguing alternative ML model. Overall, models with the input variables T_max, T_min, RH_mean, R_s, and U₂ performed better for daily RET estimation.

The current research highlights the ML and hybrid models’ efficiency in estimating daily RET within two irrigated Moroccan perimeters. Ultimately, further investigations could explore additional ML algorithms, hybrid model configurations, and their relevance to long-term datasets at different time scales and various climate regions.