The assessment of the spatial variability and uncertainty of subsoil properties is crucial for civil and environmental engineering and mining purposes. This implies that reliable subsoil models describing the 3-dimensional distribution of lithotypes and hydro-mechanical variables of interest, also known as the Engineering Geological Model (EGM), need to be defined properly. In the last forty years, a number of geostatistical methods have been implemented in Geostatistics to take into account non-Gaussian distributions, co-regionalized variables, geographic trends, and external drift that come from complex 3D spatial datasets (Heuvelink and Webster, 2022), all based on the assumption of spatial cross-correlation among measurements. Some attempts to take advantage of these techniques have been successfully defined, e.g., the 3D lithological model from Cone Penetration Data (CPTu), using stationary and non-stationary Kriging methods (Vessia et al., 2020). However, the variables of engineering geological and geotechnical interest in EGMs are rarely measured directly. Especially when using CPTu profiles as initial datasets, transformation equations are needed to obtain the target hydro-mechanical properties of EGMs. For this reason, in this work, we used a stochastic simulation approach to define reliable 3D models of some of the most used geotechnical designing variables (undrained shear resistance–su, friction angle–’, and the Darcy permeability coefficient–k) from tip resistance (qc), sleeve friction (fs), and pore pressure (u2) profiles. The selected method – the Sequential Gaussian Co-Simulation (SGCS) –, in addition to providing more reliable optimized 3D models of the spatial distribution of the variables of interest, allowed quantifying the propagation of the estimation uncertainty associated with the raw measurement models through the transformation equations. The resulting 3D subsoil models can be queried to get estimates and uncertainty, at each point in the space and updated by new CPTs investigations.