Analysis of geostatistical data through spatial process model necessitates computation and storage that becomes prohibitive with increment in number of spatial locations. Inference on spatial parameters is critical to understand the structural dependence in the data. In spatial Gaussian process regression, finite sample inference proceeds either using posterior distributions (Bayesian approach) or via resampling/subsampling techniques (frequentist). Resampling methods have become more attractive for big data as, unlike Bayesian models which require sequential sampling from MCMC, they naturally lend themselves to parallel computing resources. However, spatial boostrap involves an expensive Cholesky decomposition to decorrelate the data. In this manuscript, we develop a highly scalable parametric spatial bootstrap that uses sparse Cholesky factors for parameter estimation and decorrelation. The proposed BRISC algorithm (https://github.com/ArkajyotiSaha/BRISC) requires linear oder memory and computations, and is embarrassingly parallel, thereby delivering massive scalability. Simulation studies and real data analysis highlight the accuracy and computational efficiency of our approach.