The interpretation of in-situ tests often relies on empirical methods or analytical solutions of simplified geometrical models (e.g. cavity expansion for cone testing). Those approaches are performant for the vast majority of textbook soils, such as sands or insensitive clays, but the interpretation of the test is still challenging for non-standard materials, such as sensitive clays or mine tailings. Advances in numerical techniques now allow for the realistic simulation of in situ tests using the actual geometry and incorporating soil behavior through appropriate governing equations and advanced constitutive models, thus providing valuable insights into the mechanisms that take place during testing and, more importantly, the numerical results can be employed to propose new interpretation techniques. In this work, we propose a numerical-aided data-driven approach for the interpretation of Cone Penetration Tests (CPTu). All available data from a site-investigation campaign is analyzed and used to define a range of plausible constitutive parameters of the site. Numerical simulations of cone penetration testing with this range of constitutive parameters are then performed and the outputs of these simulations are used to define site-specific interpretation techniques. We employ this methodology to interpret a CPTu from a geotechnical research site, composed by slightly over-consolidated, sensitive clays. In this site, abundant in-situ tests and laboratory (undrained triaxial tests) data are available. Numerical simulations are performed by means of PFEM (Particle Finite Element method). A fully coupled hydromechanical formulation is employed and the problem domain is discretized using low-order finite elements. To overcome the unstable behavior in the undrained limit and potential volumetric locking, we employ a three-field mixed formulation, encompassing two stabilization terms. The constitutive response of the soil is modelled using a critical state elasto-plastic model (CASM), extended to incorporate soil structure following the approach of Gens and Nova. This model is integrated with a suitable explicit stress integration technique for large strain elasto-plastic models, based on the multiplicative decomposition of the deformation gradient into an elastic and plastic part. Using classical finite elements, the numerical solution of strain-softening dominated simulations is mesh-dependent. This is alleviated here by employing an integral-type non-local reg