Accurate representation of reservoir heterogeneity using stochastic modelling 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 identifying reservoir properties arising from our imperfect knowledge of the reservoir under study. That uncertainty is data-dependent and most importantly, model-dependent. The paper presents a Markov Chain Monte Carlo (MCMC) 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 local probability distributions obtained by kriging. This iterative updating procedure is continued until a prescribed large number of iterations are performed. Hard data, histogram, and variogram models are honoured, as expected. The multiple point histogram (MPH) and entropy of MPH are used to assess the joint spatial uncertainty exhibited by the MCMC model. The paper also discusses and demonstrates a new approach for integrating secondary information within the MCMC framework. The reduction in joint spatial uncertainty subsequent to integrating the secondary information is demonstrated.