Delineation of facies in the subsurface and quantification of uncertainty in their boundaries are significant steps in mineral resource evaluation and reservoir modeling, which impact downstream analyses of a mining or petroleum project. This paper investigates the ability of nonparametric geostatistical simulation algorithms (sequential indicator, single normal equation and filter-based simulation) to construct realizations that reproduce some expected statistical and spatial features, namely facies proportions, boundary regularity, contact relationships and spatial correlation structure, as well as the expected fluctuations of these features across the realizations. The investigation is held through a synthetic case study and a real case study, in which a pluri-Gaussian model is considered as the reference for comparing the simulation results. Sequential indicator simulation and single normal equation simulation based on over-restricted neighborhood implementations yield the poorest results, followed by filter-based simulation, whereas single normal equation simulation with a large neighborhood implementation provides results that are closest to the reference pluri-Gaussian model. However, some biases and inaccurate fluctuations in the realization statistics (facies proportions, indicator direct and cross-variograms) still arise, which can be explained by the use of a single finite-size training image to construct the realizations.