No single spatial interpolation method reigns supreme for modelling the precise spatial distribution of groundwater quality data. This study addresses this challenge by evaluating and comparing several commonly used geostatistical methods: Local Polynomial Interpolation (LPI), Ordinary Kriging (OK), Simple Kriging (SK), Universal Kriging (UK), and Empirical Bayesian Kriging (EBK). We applied these methods to a vast dataset of 3033 groundwater records encompassing a substantial area (11,100 km2) in the coastal lowlands of the western Netherlands. To our knowledge, no prior research has investigated these interpolation methods in this specific hydrogeological setting, exhibiting a range of groundwater qualities, from fresh to saline, often anoxic, with high natural concentrations of PO4 and NH4. The prediction performance of the interpolation methods was assessed through statistical indicators such as root means square error. The findings indicated that EBK outperforms the other geostatistical methods in forecasting groundwater quality for the five variables considered: Cl, SO4, Fe, PO4, and NH4. In contrast, SK performed worst for the species except for SO4. We recommend not using SK to interpolate groundwater quality species unless the data exhibit low spatial variation, high sample density, or evenly distributed sampling.