Reconstructing the variation of contaminant concentration with a limited number of soil samples is more or less the norm, even though it fails more often than not for problems of even moderate complexity. To overcome the limits inherent to discrete measurements, we propose to integrate soil sampling with continuous surface geophysical measurements in a geostatistical framework. We present this integrated analysis for a PAH contaminated site in France. For the study site, two 3D surveys were acquired: an electrical resistivity tomography survey and a seismic travel time tomography survey. Those two surveys permitted us to infer two spatially continuous physical properties on the whole volume, namely the electrical resistivity and P-wave velocity. The probability density function relating the velocity-resistivity pairs with each of the 75 lab measurements of PAH concentration was modeled using a Gaussian kernel. This probability density function combined with the 3D volumes of resistivity and P-wave velocity provided a means to translate the latter into a 3D map of PAH concentration. This 3D map of concentration was then used as a secondary variable in a cokriging simulation of the 75 lab samples, thus reintroducing the spatial correlation of the initial dataset. Comparing this final 3D PAH concentration model with the simple kriging of the PAH samples, the geophysical integrated model reproduce much better the distribution of measured concentration, shows a much more realistic spatial pattern of the contamination, and lowers the estimated contaminated volume.