Plurigaussian simulation, an improvement of truncated Gaussian method, is a new and effective procedure of reproducing complex geometrical attributes of a reservoir by means of simulating several facies with different spatial structures (anisotropies, for example) and taking into account their global proportions. The method is based on the relationship between indicator variables from facies distribution and the Gaussian random functions chosen to represent them, with the compromise of having equivalent geostatistical behavior. Geological information is incorporated to the model by the definition of the indicator variables, their truncation strategy and the facies proportions. Two critical points arise in this methodology: the relationship between the indicator and the Gaussian variograms during the structural modeling and the correspondence between indicator and continuous variable values during the conditioning. To face them, we use iterative fitting and Gibbs sampler algorithms. As illustration, the method was applied to a training image from a 2D geological interpretative map of the Aux Vases Zone in the Rural Hill Field - Illinois – USA, from which two sets of samples with different grid densities were extracted. Previous knowledge of the facies geometry permitted us to adopt a truncation strategy performed in both stationary and non-stationary approaches.Plurigaussian simulation, an improvement of truncated Gaussian method, is a new and effective procedure of reproducing complex geometrical attributes of a reservoir by means of simulating several facies with different spatial structures (anisotropies, for example) and taking into account their global proportions. The method is based on the relationship between indicator variables from facies distribution and the Gaussian random functions chosen to represent them, with the compromise of having equivalent geostatistical behavior. Geological information is incorporated to the model by the definition of the indicator variables, their truncation strategy and the facies proportions. Two critical points arise in this methodology: the relationship between the indicator and the Gaussian variograms during the structural modeling and the correspondence between indicator and continuous variable values during the conditioning. To face them, we use iterative fitting and Gibbs sampler algorithms. As illustration, the method was applied to a training image from a 2D geological interpretative map of the Aux Vases Zone in the Rural Hill Field - Illinois – USA, from which two sets of samples with different grid densities were extracted. Previous knowledge of the facies geometry permitted us to adopt a truncation strategy performed in both stationary and non-stationary approaches.