Quantifying uncertainties in subsurface properties and stratigraphy can lead to better understanding of the ground conditions and enhance the design and assessment of geotechnical structures. Several studies have utilized Cone Penetration Test (CPT) data and employed Bayesian and Machine Learning methods to quantify the geological uncertainty, based on the Robertson’s soil classification charts and the Soil Behaviour Type Index (Ic). The incorporation of borehole data can reduce the stratigraphic uncertainty. Significant challenges can arise, however, mainly due to the intrinsic differences between field and laboratory-based soil classification systems, which can potentially lead to inconsistent soil classification. To this end, this study proposes a multivariate Gaussian Process model that utilizes site-specific data and: i) jointly models multiple categorical (USCS labels) and continuous (Ic) variables, ii) learns a (shared) spatial correlation structure and the between-outputs covariance, and iii) produces two types of dependent classification outputs. The results indicate that the integration of geotechnical and geological information into a unified model can provide more reliable predictions of the subsurface stratification, by allowing simultaneous interpretation of USCS and Ic profiles. Importantly, the model demonstrates the potential to integrate multiple variables of different types, aiming to contribute to the development of a methodology for joint modeling of geotechnical, geological and geophysical data.