A new geostatistical grid-free simulation (GFS) method has been developed recently that represents simulated continuous attributes of the natural phenomena, such as stratigraphic surface boundaries or petrophysical properties, as an analytical function of the coordinates of the simulation locations. Thus, GFS resolves challenges related to model regridding, increasing resolution around already simulated locations and integration of newly available data in a consistent manner. The present paper contains further developments in simulation of categorical variables, such as facies, in a grid-free fashion based on the truncated pluri-Gaussian simulation (TPG) paradigm. The resultant simulation engine allows the entire reservoir system to be represented as an analytical stochastic function: that is, values of any reservoir properties are simulated on demand at requested locations in space. The selection of proper variograms of Gaussian continuous variables for simulation of categorical variables in a TPG framework is proposed through a methodology based on Monte Carlo simulation. The variogram models of the underlying Gaussian continuous variables are obtained by minimizing the difference between numerically computed and target indicator variograms. A local optimization approach is suggested for a fast precise derivation of the variograms. The stable variogram model leads to the closest fit to the experimental variograms of the continuous variables. The automatic establishment of a truncation mask based on multidimensional scaling to convert Gaussian continuous variables to categories is also explained. Finally, the proposed GFS algorithm for petroleum reservoir characterization is demonstrated in its full-scale applicability to the Hekla offshore petroleum reservoir located in the North Sea. The results look promising, and should be beneficial to petroleum reservoir modelling in practice.