Geostatistical simulation methods are able to generate numerical models or relations of the spatial distribution of a continuous geologic variable (grade, thickness, density, etc.) or a categorical variable (geological units and lithofacies or rock types). In this work, a review of traditional simulation techniques, as the Sequential Indicator Simulation (SIS), reveals a major pitfall that comes from theoretical difficulties in the development of a valid cross covariance model. On the contrary, a valid indicator cross covariance model is automatically defined in the framework of the Truncated Gaussian Simulation Method (TGS). This method is based on the concept that the categorical variables are obtained by truncating one standard multigaussian random variable at different thresholds. Plurigaussian Simulation Method (PGS) is an extension of the TGS Method but based on the simultaneous truncation of several multigaussian variables. An application of Plurigaussian method to simulate the lithofacies in the alluvial formations of the West Thessaly Basin is finally presented. This method was shown to be effective in reproducing the spatial characteristics of the different lithofacies and their distribution across the studied area.