‎Analysis of large geostatistical data sets‎, ‎usually‎, ‎entail the expensive matrix computations‎. ‎This problem creates challenges in implementing statistical inferences of traditional Bayesian models‎. ‎In addition,researchers often face with multiple spatial data sets with complex spatial dependence structures that their analysis is difficult‎. ‎This is a problem for MCMC sampling algorithms that are commonly used in Bayesian analysis of spatial models‎, ‎causing serious problems such as slowing down and chain integration‎. ‎To escape from such computational problems‎, ‎we use low-rank models‎, ‎to analyze Gaussian geostatistical data‎. ‎This models improve MCMC sampler convergence rate and decrease sampler run-time by reducing parameter space‎. ‎The idea here is to assume‎, ‎quite reasonably‎, ‎that the spatial information available from the entire set of observed locations can be summarized in terms of a smaller‎, ‎but representative‎, ‎sets of locations‎, ‎or ‘knots’‎. ‎That is‎, ‎we still use all of the data but we represent the spatial structure through a dimension reduction‎. ‎So‎, ‎again‎, ‎in implementing the reduction‎, ‎we need to design the knots‎. ‎Consideration of this issue forms the balance of the article‎. ‎To evaluate the performance of this class of models‎, ‎we conduct a simulation study as well as analysis of a real data set regarding the quality of underground mineral water of a large area in Golestan province‎, ‎Iran‎.