An accurate characterization of the spatial distribution of the near-surface geological properties is often challenging because of their heterogeneity, yet essential for different activities. Non-invasive geophysical surveys have been proven powerful tools for the collection of virtually spatially continuous high-resolution datasets that can be translated in detailed images of subsurface physical properties. Electrical and electromagnetic near-surface techniques, particularly frequency-domain electromagnetic (FDEM) and electrical resistivity tomography (ERT) methods have demonstrated their efficiency to characterize heterogeneous subsurface systems. Both methods are sensitivity to electrical conductivity (EC), and FDEM data can additionally be linked to the magnetic susceptibility (MS). FDEM and ERT data can be translated into numerical subsurface models by geophysical inverse problem, that is an ill-posed nonlinear problem with a nonunique solution, resulting in uncertainties in the inverse models. Inverting both data sets from these different geophysical methods in a joint inversion methodology is generally a more preferable approach due to the difference in the spatial resolution of both methods and the sensitivity to different physical properties. In this work, we present an iterative geostatistical joint inversion technique for FDEM and ERT, that is capable of modelling the spatial continuity and assess the uncertainty of the modelling procedure.