Abstract 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. Introduction The assessment of uncertainty is a crucial step for decisionmakers to reduce the financial risks during the development of a reservoir. The uncertainty in reservoir performance predictions is usually reflected by a number of parameters of interest, such as OIP, production profiles, ultimate recovery, and any forecast figure, etc. Reservoir complexity and the lack of data available are two main reasons that lead to uncertainty in reservoir modeling(1). Several papers(1–4) have addressed the issue of uncertainty assessment in detail. However, all efforts for quantifying uncertainty are model and data specific. The true reservoir permeability field is entirely characterized by a multivariate probability distribution describing the joint uncertainty in permeability at all locations in the reservoir. That multivariate distribution is usually unknown and hence the objective of stochastic simulation is to obtain the requisite multivariate distribution using the available data and the prior model for spatial continuity (or variability). Ultimately, this posterior probability distribution is sampled in order to obtain realizations of the permeability field. In this paper, a Markov Chain Monte Carlo (MCMC) procedure is used to update and sample from the posterior probability density function. The simulation commences from an initial kriged map. The local conditional probability distribution at each location is constructed using the values at neighbouring locations and subsequently, a value is sampled from that distribution. A Gibbs sampler is utilized implying that the sampled value is retained regardless of whether the mismatch against target statistics increases. The procedure is continued at all locations by following a random path. This iterative procedure is continued until a final realization is obtained that identifies target statistics, such as the prior variogram and/or the prior histogram of permeability values. The details of the MCMC algorithm are described next.