Accurate representation of reservoir heterogeneity using stochastic modeling techniques requires careful synthesis of the multivariate probability distribution characterizing the spatial variations of petrophysical properties. A well designed scheme for sampling from that distribution is necessary in order to accurately portray the uncertainty in estimating reservoir properties arising from the incomplete knowledge of the reservoir under study. That uncertainty is data-dependent and most importantly model-dependent. The paper presents a Markov Chain Monte Carlo methodology for sampling from the invariant or stationary probability distribution characterizing the reservoir permeability field. An initial random field is perturbed successively following a Gibbs sampling procedure. The updated values at the perturbed node are obtained by sampling from the corresponding kriged distribution. This iterative updating procedure is continued until a prescribed large number of iterations are performed. Hard data, histogram and variogram model are honored as expected. The multiple point histogram (MPH) and entropy of MPH are used to assess the joint spatial uncertainty exhibited by this MCMC model and the impact of data quantity and configuration on that uncertainty are explored. Convergence characteristics of the proposed methodology are investigated and detailed comparisons of the proposed methodology to other established algorithms are presented.