sampling Gaussian process (GP) models are widely used in the Analysis of computer experiments. However, two issues have not been solved satisfactorily. One is the computational issue that hinders GP from broader application, especially for massive data observed on irregular grids. The other is the underestimation of the GP prediction uncertainty. To tackle these problems, we introduce a new version of bootstrap using design-based subsampling and proposed two methods to construct bootstrap predictive distributions. The new bootstrap procedure borrows the strength of space-filling designs to provide an efficient subsample and reduce computational complexity. It is shown that this procedure not only alleviates the computational difficulty in GP estimation and prediction, but also provides an unbiased predictor and an efficient and accurate quantification of uncertainty. Theoretical properties of the bootstrap predictors are discussed. The finite-sample performance is examined through simulations. We illustrate the proposed method by an ice sheet computer experiment.