Gaussian process models have been widely used in spatial statistics but face tremendous computational challenges for very large data sets.The model fittiug and spatial prediction of such models typically require O(n3) operations for a data set of size n.Various approximations of the covariance functions have been introduced to reduce the computational cost.However, most existing approximations can not simultaneously capture both the large and small scale spatial dependence.