Confronting uncertainty in model-based geostatistics using Markov Chain Monte Carlo simulation

This paper demonstrates for the first time the use of Markov Chain Monte Carlo (MCMC) simulation for parameter inference in model-based soil geostatistics. We implemented the recently developed DiffeRential Evolution Adaptive Metropolis (DREAM) algorithm to jointly summarize the posterior distribution of variogram parameters and the coefficients of a linear spatial model, and derive estimates of predictive uncertainty. The DREAM method runs multiple different Markov chains in parallel and jumps in each chain are generated from a discrete proposal distribution containing a fixed multiple of the difference of the states of randomly chosen pairs of other chains. This approach automatically scales the orientation and scale of the proposal distribution, and is especially designed to maintain detailed balance and ergodicity, thereby generating an exact approximation of the posterior probability density function (pdf) of the parameters of the linear model and variogram. This approach is tested using three different data sets from Australia involving variogram estimation of soil thickness, kriging of soil pH, and spatial prediction of soil organic carbon content. The results showed some advantages of MCMC over the conventional method of moments and residual maximum likelihood (REML) estimation. The posterior pdf derived with MCMC conveys important information about parameter uncertainty, multi-dimensional parameter correlation, and thus how many significant parameters are warranted by the calibration data. Parameter uncertainty constitutes only a small part of total prediction uncertainty for the case studies considered here. The prediction accuracies using MCMC and REML are similar. The variogram estimated using conventional approaches (method of moments, and without simulation) lies within the 95% prediction uncertainty interval of the posterior distribution derived with DREAM. Altogether our results show that conventional kriging and regression-kriging still remain a viable option for production mapping.