Geostatistical tools have been increasingly applied to compile spatial distribution maps of groundwater quality. Results of geostatistics can be helpful for decision-makers to carry out appropriate remedial actions to sustain the quality of groundwater sources. The main purpose of this paper is to assess the groundwater hydrogeochemistry of Tehran-Karaj plain aquifer, Iran, by a forward geostatistical method and an inverse geostatistical method called sequential Gaussian simulation. Seven important hydrogeochemical properties of the aquifer including total dissolved solids, sodium adsorption ratio, electrical conductivity, sodium, total hardness, chloride, and sulfate were analyzed and compiled geostatistically. Data were taken from 137 well samples in 2016. After data normalization, variography was compiled, and experimental variograms were plotted; then, the best theoretical model was fitted on each variogram based on the minimum residual sum of squares. Cross-validation was used to determine the accuracy of parameters related to the variograms. Estimation maps of the groundwater hydrogeochemistry were prepared, and the estimation variance map was drawn to assess the accuracy of estimation in each estimated point. Forward geostatistical methods are subjected to smoothing, while inverse geostatistical methods are not subjected to this problem. The results of this study revealed that the utilized inverse geostatistical methods called simulation algorithms are more accurate than forward methods. Eventually, estimation maps of each parameter, as well as error maps, were compiled, and critical regions have been proposed according to the simulated maps.