Annually calibrated k_Rs values are presented in Table 3, with values typically lower than the original k_Rs. All weather stations show high R² values (0.97–0.99), indicating strong correlations between ET_o–HS and ET_o–PM (Figure 10a and Figure 11a). The RMSE values are relatively low (0.45–0.73 mm d⁻¹), with MT having the lowest values (0.46 mm d⁻¹ and 0.45 mm d⁻¹ for calibration and validation, respectively). Annual calibration minimized the error by approximately 0.10 mm d⁻¹. The PBIAS values were generally lower for all stations, indicating low bias in the ET_o–HS estimates. The exception is DV, where the validation has a PBIAS of −1.51%, indicating a slight underestimation. Annual calibration reduced bias by approximately 6% compared to the original calibration of HS. The EF values are high for all stations (0.90–0.96), with values close to 1, indicating excellent ET_o–HS estimation accuracy. The DV weather station has a relatively higher RMSE (0.71–0.73 mm d⁻¹) compared to other stations, while the MT station has a relatively lower RMSE (0.46–0.45 mm d⁻¹) and NRMSE (12.79–13.32%), indicating the most accurate adjustment. The validation results support the robustness and reliability of the calibration. Although the accuracy changed slightly, the annual adjusted k_Rs were similar for most stations analyzed.

The results of the goodness-of-fit indicators are consistent with those in the literature [6,22,32]. For example, Paredes et al. [6] reported an average RMSE of 0.60 ± 0.15 mm d⁻¹ for the HS calibration for different climates. Moratiel et al. [18] present similar results for the Duero basin (0.69 mm d⁻¹), while Rodrigues and Braga [57] found higher values (0.83 mm d⁻¹) for the Alentejo region. Similarly, in Iran, Raziei and Pereira [32] found that the HS model performed well with RMSE values generally below 0.70 mm d⁻¹ across most locations. The model’s performance in the windy and humid Azorean climate, as reported by Paredes et al. [19], showed RMSE values ranging from 0.47 to 0.86 mm d⁻¹, with satisfactory accuracy. This suggests that in regions where the effects of wind are minimal, the HS model can reliably estimate ET_o with relatively low error rates. Conversely, in more arid and semiarid regions, the performance of the HS model deteriorates. Ren et al. [84] reported RMSE values ranging from 0.65 to 1.15 mm d⁻¹ in arid areas of Inner Mongolia. Under extreme conditions, the accuracy of the model further decreases, as Zhu et al. [85] found in their evaluation of 838 weather stations in China, particularly the HS model, which neglects relative humidity and wind speed. The limitations of the HS model become particularly evident during the dry season in the São Francisco River Basin, Brazil, as reported by Althoff [66], where lower correlations between meteorological variables and ET_o reduce the effectiveness of the model. Additionally, the high variability in wind speed within certain sub-regions further diminishes the model’s accuracy.

The K-Means algorithm defined the optimal number of clusters into which the studied months and years are divided. Therefore, with the same number of clusters (k = 2) for the annual and monthly data, the total sum of squares for monthly data decreases considerably across different temporal groupings, particularly highlighting seasonal and annual variations (Figure 9). The analysis showed that grouping by months produced better results than by years, as indicated by the higher average silhouette coefficient for monthly data (0.44) compared to annual data (0.34). This suggests that monthly clusters provided more coherent and separated groups and optimized the partitioning of meteorological data. Furthermore, the total within-cluster sum of squares decreased significantly for monthly data, reinforcing the notion that seasonal variations are crucial for improving model calibration.

The clusters created for monthly data aligned well with seasonal patterns, dividing the year into a wet season (clusters 11 and 21) and the dry season from May to September (clusters 12 and 22). This clustering coincides with periods of different evapotranspiration demands, especially since irrigation requirements are primarily relevant during the dry season. This reinforces the potential benefit of calibrating the Hargreaves–Samani equation for the wet and dry seasons, as suggested by Zanetti et al. [86], Aguilar and Polo [87], and Althoff et al. [66]. The results for the annual clusters can be found in the Appendix A (Table A3).

The cluster-based calibration using year and month clusters focused on four groups: years below the threshold (11 and 12), years above the threshold (21 and 22), wet months (11 and 21), and dry months (12 and 22). Most clusters had high R² values (0.94–0.99), indicating strong correlations between ET_o–HS and ET_o–PM (Figure 10b and Figure 11b). However, the accuracy varied depending on the cluster and weather station. For instance, the MR weather station consistently showed the lowest RMSE (0.32 mm d⁻¹) and NRMSE (9.2%).

On the other DV, the weather station experienced high errors in certain clusters (e.g., cluster 21 with NRMSE of 27.46%) and biases (e.g., cluster 11 with PBIAS of −3.06%) during validation. The MT weather station shows consistently high performance, although minor declines were observed in some clusters during validation. In contrast, VA shows moderate performance, with certain clusters experiencing declines in validation (cluster 11 with NRMSE of 24.7% and PBIAS of 4.87%) (Table 4).

In the context of climate change studies, capturing long-term trends and variability in climate patterns is critical. The analysis focused on two sets of clusters (Table 4): clusters representing years below the threshold (11 and 12) and clusters representing years above the threshold (21 and 22) (Figure 10b and Figure 11b).

The analysis of the monthly clusters showed better performance during dry month clusters 12 and 22 (Table 4), with lower NRMSE values (varying between 9.2 and 16.96% for calibration), suggesting lower ET_o−HS estimation errors with little influence of the annual clusters. The higher ET_o values during the dry months naturally lead to higher RMSE values compared to the wet months, but when normalized, the NRMSE values were lower, which suggests that the model was able to capture the higher evapotranspiration rates more effectively in drier conditions. The results also suggest that the model was more reliable for some weather stations, such as MR, than for others, like DV, underscoring the importance of localized calibration.

These results suggest better performance for dry-month clusters, which is consistent with the results of Moratiel et al. [18]. However, these authors also present acceptable results for annual calibration (RMSE < 0.69 mm d⁻¹) despite poorer performance in winter. These results are also coherent with the work of Teixeira et al. [88], in which they showed that the relative error of the Hargreaves–Samani method for the Alentejo region decreased significantly to about 12% in the summer. Similar conclusions were reached by Shahidian et al. [44], for two different regions, California and Bolivia, where they showed that the annual correlations between Hargreaves–Samani and Penman–Monteith are not the most correct method to perform calibration because this correlation tends to increase in the dry season compared to the wet season, resulting in significant variation in model performance for these two periods.

Rodrigues and Braga [56] calibrated the HS equation for the irrigation season and presented an RMSE of 0.88 mm d⁻¹. In addition to the previous work [56] that focused on k_Rs calibration specifically for irrigation periods, the present study uses calibration approaches across all months and meteorological conditions (dry and wet seasons). This is important as it allows the application of HS for the entire cropping pattern, which includes winter crops (e.g., fodder crops and wheat). Furthermore, it is also intended to include inter-annual climate variability in this study to assess whether the calibrated k_Rs accurately reproduce ET_o for dry and wet years, which is very relevant for climate change studies. This broader scope increases the applicability of the methodology to other climates. Another contribution of our study that was not considered in the study [56] is the quality and consistency of meteorological data across the stations, which may affect the model results.

An accuracy improvement was expected when data were clustered by years and months, compared to annual calibration. Wet season clusters performed only slightly (1%) or even worse (5% below the annual calibration), while dry season clusters also performed slightly better (up to 3%) (Figure 12). The performance of the HS equation is often reported in the literature. For humid conditions, some studies [53,89,90] showed that the HS equation performs poorly, as observed in the present study, while other studies reported favorable results for humid conditions [22,32,53,84,89,90].

Note that solar radiation and relative humidity have a stronger correlation with ET_o in spring (Figure 8), and these variables are not considered in the HS equation. This could explain the worse performance of HS during wet season clusters. The maximum RMSE improvement for the cluster approach over the annual approach was 0.13 mm d⁻¹. According to Althoff [66], calibration of the HS equation for specific climatic periods or regions can improve its accuracy. However, this does not guarantee that predictions will be considerably better for every period or region. They found that minimal improvement was achieved by cluster calibration in tropical regions including the Cerrado, Caatinga, and Atlantic Forest with semi-arid to humid subtropical climates. Additionally, they suggested that calibration results using annual k_Rs provided good results, consistent with our observations that seasonal cluster-based calibration provided minimal improvements compared to simpler methods (annual k_Rs). The annual calibration approach reduces the error and bias of ET_o−HS estimates for all weather stations compared to the original HS equation. However, further refinement using cluster-based calibration, which accounts for annual and monthly patterns in the meteorological variables, produced unclear results. While for dry season clusters results showed small improvements compared to the annual approach, the gains do not justify the increased complexity of the cluster-based approach.