The assessment of well head protection zones is strongly related to spatial variability of hydrogeological properties. To set up a numerical model for probabilistic wellhead protection zone delineation, field and laboratory scale experiments for subsurface investigation were performed at the “Lauswiesen” site in southern Germany. At this site, spatial variations of hydraulic conductivity mainly dominate flow and transport processes. By means of a multivariate cluster analysis, three different types of aquifer materials (facies) were identified to characterize the heterogeneity of the aquifer lithology. Spatial variability of each material type was separately analysed by a geostatistical procedure based on indicator variography. Three-dimensional recombined distributions of the identified clusters were used to reconstruct the internal aquifer architecture. Then, three-dimensional time related mean well head protection zones and their associated uncertainty were obtained by a Monte Carlo stochastic simulation procedure