Quantifying the proportions and distribution of rock types is a key step in modeling sedimentary formations, since it greatly determines the quality of petrophysical and dynamic flow models. While forward stratigraphic models well reproduce deposition processes, they rely on inversion for data conditioning, which critically raises the empty space problem. Therefore, several stochastic methods have been suggested to simulate proportions and distribution of rock types in reservoirs respecting the observation data while reproducing realistic spatial patterns. However, these methods are always limited by the stationarity decision. We propose a method to generate 3D facies probability cubes which account for well and seismic data, stratigraphic interpretation, sedimentological rules describing the spatial distribution of rock types, and sequence stratigraphy principles. Different probability cubes are computed by integrating various information controlling the lithofacies occurrence, for instance the stratigraphic control of the shoreline migration, by defining linear or inequality constraints at the spatial estimation stage. The generated P-field cubes can then be combined considering the redundancy of the data they express such as the tau model. This methodology is demonstrated on typical synthetic data sets and on a North Sea reservoir.