Summary An accurate characterisation of the spatial distribution of near-surface geological properties is often challenging due to their heterogeneity, but is essential for different activities. Non-invasive geophysical surveys have proven to be powerful tools for collecting high-resolution, virtually space-continuous datasets that can be translated into numerical models of subsurface physical properties through an inverse geophysical problem, which is a poorly posed nonlinear problem with a non-unique solution, resulting in uncertainties in the inverse models. Frequency domain electromagnetic (FDEM) and electrical resistivity tomography (ERT) methods, have demonstrated their efficiency in characterising heterogeneous subsurface systems. Both methods are sensitive to electrical conductivity (EC), and FDEM data can be additionally coupled with magnetic susceptibility (MS). Inversion of both datasets from these different geophysical methods in a joint inversion methodology is generally a more preferable approach due to the difference in spatial resolution of both methods and sensitivity to different physical properties. In this work, we present the application of an iterative geostatistical joint inversion technique for FDEM and ERT in a real data set, which is able to model the spatial continuity and assess the uncertainty of the modelling procedure.