Water shortage and quality are major issues in many places, particularly arid and semi-arid regions such as Makkah Al-Mukarramah province, Saudi Arabia. The current work was conducted to examine the geochemical mechanisms influencing the chemistry of groundwater and assess groundwater resources through several water quality indices (WQIs), GIS methods, and the partial least squares regression model (PLSR). For that, 59 groundwater wells were tested for different physical and chemical parameters using conventional analytical procedures. The results showed that the average content of ions was as follows: Na+ > Ca2+ > Mg 2+ > K+ and Cl− > SO42− > HCO32− > NO3− > CO3−. Under the stress of evaporation and saltwater intrusion associated with the reverse ion exchange process, the predominant hydrochemical facies were Ca-HCO3, Na-Cl, mixed Ca-Mg-Cl-SO4, and Na-Ca-HCO3. The drinking water quality index (DWQI) has indicated that only 5% of the wells were categorized under good to excellent for drinking while the majority (95%) were poor to unsuitable for drinking, and required appropriate treatment. Furthermore, the irrigation water quality index (IWQI) has indicated that 45.5% of the wells were classified under high to severe restriction for agriculture, and can be utilized only for high salt tolerant plants. The majority (54.5%) were deemed moderate to no restriction for irrigation, with no toxicity concern for most plants. Agriculture indicators such as total dissolved solids (TDS), potential salinity (PS), sodium absorption ratio (SAR), and residual sodium carbonate (RSC) had mean values of 2572.30, 33.32, 4.84, and −21.14, respectively. However, the quality of the groundwater in the study area improves with increased rainfall and thus recharging the Quaternary aquifer. The PLSR models, which are based on physicochemical characteristics, have been shown to be the most efficient as alternative techniques for determining the six WQIs. For instance, the PLSR models of all IWQs had determination coefficients values of R2 ranging between 0.848 and 0.999 in the Cal., and between 0.848 and 0.999 in the Val. datasets, and had model accuracy varying from 0.824 to 0.999 in the Cal., and from 0.817 to 0.989 in the Val. datasets. In conclusion, the combination of physicochemical parameters, WQIs, and multivariate modeling with statistical analysis and GIS tools is a successful and adaptable methodology that provides a comprehensive picture of groundwater quality and governing mechanisms.
