Reliable prediction of subsurface flow and contaminant transport depends on the accuracy with which the values and spatial distribution of process-relevant model parameters can be identified. Successful characterization methods for complex soil systems are based on (i) an adequate parameterization of the subsurface, capable of capturing both random and structured aspects of the heterogeneous system, and (ii) site-specific data that are sufficiently sensitive to the processes of interest. We present a stochastic approach in which the high-resolution imaging capability of geophysical methods is combined with the process-specific information obtained from the inversion of hydrological data. Geostatistical concepts are employed as a flexible means to describe and characterize subsurface structures. The key features of the proposed approach are (i) the joint inversion of geophysical and hydrological raw data, avoiding the intermediate step of creating a (nonunique and potentially biased) tomogram of geophysical properties and (ii) the concurrent estimation of hydrological and petrophysical parameters, in addition to (iii) the determination of geostatistical parameters from the joint inversion of hydrological and geophysical data. This approach is fundamentally different from inference of geostatistical parameters from an analysis of spatially distributed property data. The approach has been implemented into the iTOUGH2 inversion code and is demonstrated for the joint use of synthetic time-lapse ground penetrating radar (GPR) travel times and hydrological data collected during a simulated ponded infiltration experiment at a highly heterogeneous site.