Knowledge of the stratigraphic architecture and geotechnical properties of surficial soil sediments is essential for geotechnical risk assessment. In the Saguenay study area, the Quaternary deposits consist of a basal till layer and heterogeneous post-glacial deposits. Considering the stratigraphic setting and soil type heterogeneity, a multistep stochastic methodology is developed for 3D geological modelling and quantification of the associated uncertainties. This methodology is adopted for regional studies and involves geostatistical interpolation and simulation methods. Empirical Bayesian kriging (EBK) is applied to generate the bedrock topography map and determine the thickness of the till sediments and their uncertainties. The locally varying mean and variance of the EBK method enable accounting for data complexity and moderate nonstationarity. Sequential indicator simulation is then performed to determine the occurrence probability of the discontinuous post-glacial sediments (clay, sand and gravel) on top of the basal till layer. The individual thickness maps of the discontinuous soil layers and uncertainties are generated in a probabilistic manner. The proposed stochastic framework is suitable for heterogeneous soil deposits characterised with complex surface and subsurface datasets.