We study spatial sampling design for prediction of stationary isotropic Gaussian processes with estimated parameters of the covariance function. The key issue is how to incorporate the parameter uncertainty into design criteria to correctly represent the uncertainty in prediction. Several possible design criteria are discussed that incorporate the parameter uncertainty. A simulated annealing algorithm is employed to search for the optimal design of small sample size and a two-step algorithm is proposed for moderately large sample sizes. Simulation results are presented for the Matérn class of covariance functions. An example of redesigning the air monitoring network in EPA Region 5 for monitoring sulfur dioxide is given to illustrate the possible differences our proposed design criterion can make in practice.