A geostatistical model was developed and applied to predict six soil properties and soil horizon thickness for mineral A, B and C soil horizons at the European scale and quantify the associated prediction uncertainties. The soil properties are pH, organic carbon content, organic nitrogen content, clay and sand contents and bulk density. The geostatistical model takes a regression cokriging approach, in which correlations between soil properties and across soil horizons are taken into account. Non-stationarities in the means and variances are represented by mapping units of the generalised European soil and land cover maps. The model was calibrated using the combined WISE, SPADE 1 and EFSDB databases, which jointly contain approximately 3600 soil profiles, irregularly distributed over Europe. The resulting model showed for most soil properties strong dependencies on soil type and land cover, moderate correlations between soil property residuals, strong correlations across horizons, and moderate spatial correlation of regression residuals. Kriging predictions and simulations were made on a 5 km by 5 km grid. Uncertainties in the resulting maps are large, particularly in under-sampled parts of Europe and in strata with large spatial variation. We conclude that geostatistical prediction and simulation are useful techniques to quantify uncertainties in soil property maps at the European scale, but that many more observations are required to fully exploit the relationship with explanatory variables and improve mapping accuracy. One important advantage of the techniques used is that they yield a full probabilistic model, as required by Monte Carlo uncertainty propagation analyses of spatially distributed dynamic models that use soil properties as uncertain input. In particular, the results of this study have been used to analyse how uncertainty in soil properties propagate through terrestrial greenhouse gas emission models.