In many earth sciences applications, the geological structures to be reproduced are curvilinear. Although locally accurate in a minimum error variance sense, traditional methods based on some form of kriging are not sensitive to the presence of strings of connectivity of low or high values. In such cases, techniques allowing reproduction of multiple-point statistics involving jointly three or more points at a time are required (Krishnan and Journel, 2003). However, estimation problems must be solved, taking into account that data at sparse sample localizations are not enough to infer such multiple statistics. As regards this aspect, information about geological structures prevailing in the subsurface is crucial, since multiple-point statistics can be inferred from training images reflecting the expected patterns of geological heterogeneities. The aim of this paper is to show the flexibility of integrating such multiple-statistics in a model in order to allow shape reproduction, as well as heterogeneity patterns, in the subsurface. The single normal equation methodology provides a solution. A comparison between a traditional variogram based simulation algorithm and a multiple-point statistics algorithm is presented through a case study. It is shown that limestone spatial distribution (in Lecce district, Italy) with meandering channels is better reproduced by using the latter algorithm. Multiple-point statistics are scanned directly from an available geological map, used as a training image. A sample data used in this analysis is extracted from an exhaustive limestone data set of 24642 points.