This study is the first attempt to characterize the quality status of El Fahs aquifer by combining graphical tools, multivariate statistical techniques and traditional geostatistical methods. Water samples are collected from thirty-six observation wells during April 2016 to characterize the physicochemical properties of the aquifer. Subsequently, these samples are partitioned into three hydrochemically distinct water classes (i.e., C1, C2, and C3) using the K-means clustering method. Principal Component Analysis is used to reduce the dimensionality of the dataset prior performing the clustering computations, resulting in clusters of higher quality than the non-reduced case in terms of Silhouette coefficient. Piper diagram is used to display the chemical composition of the samples, revealing the dominant role of Mg–Ca–Cl water type for all three classes, whereas Sodium and Sulfate were found to be the second most important cations and anions respectively. Indicator kriging (IK) is used to identify the probability of occurrence of the hydrochemical classes beyond the sampling locations. It is found that Class 1, associated with fresh groundwater component, is most probable to occur at the central part of the plain, mainly due to the presence of a dense hydrological network, whereas Classes 2 (agricultural activities) and 3 (dissolution of evaporate geological formations) are expected to occur at the southern and northern regions respectively. IK also identified the regions associated with high levels of uncertainty, mostly occurring in a large portion of the northern area due to the absence of available hydrochemical information. The results showed that integration of graphical methods, multivariate statistical techniques and geostatistical modeling, is an efficient approach for characterizing the hydrochemical status of the aquifer system, to spatially optimize the groundwater monitoring well networks and quantify the uncertainty levels of the water classes in a systematic way.