Abstarct: 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‎. ‎Complex structures of (spatial or spatio-temporal) data‎, ‎their sizes and available incomplete models‎, ‎have been pushed users to develope new and dynamic models‎. ‎In this thesis‎, ‎we use low-rank models‎, ‎particularly predictive process models‎, ‎to analyze Gaussian geostatistical data‎. ‎This models improve MCMC sampler convergence rate and decrease sampler run-time by reducing parameter space‎. ‎This improvement is due to the avoidance of expensive matrix computations that also can be used for large data sets‎. ‎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.