Previous studies have demonstrated that Visible Near InfraRed (Vis-NIR) hyperspectral imagery is a cost-efficient way to map soil properties at fine resolutions (similar to 5 m) over large areas. However, such mapping is only feasible for the soil surface because the effective penetration depths of optical sensors do not exceed several millimeters. This study aims to determine how Vis-NIR hyperspectral imagery can serve to map the subsurface properties at four depth intervals (15-30 cm, 30-60 cm, 60-100 cm and 30-100 cm) when used with legacy soil profiles and images of parameters derived from digital elevation model (DEM). Two types of surface-subsurface functions, namely linear models and random forests, that estimate subsurface property values from surface values and landscape covariates were first calibrated over the set of legacy measured profiles. These functions were then applied to map the soil properties using the hyperspectral-derived digital surface soil property maps and the images of landscape covariates as input. Error propagation was addressed using a Monte Carlo approach to estimate the mapping uncertainties. The study was conducted in a pedologically contrasted 300 km(2)-cultivated area located in the Cap Bon region (Northern Tunisia) and tested on three soil surface properties (clay and sand contents and cation exchange capacity). The main results were as follows: i) fairly satisfactory (cross-validation R-2 between 0.55 and 0.81) surface-subsurface functions were obtained for predicting the soil properties at 15-30 cm and 30-60 cm, whereas predictions at 60-100 cm were less accurate (R-2 between 0.38 and 0.43); ii) linear models outperformed random-forest models in developing surface-subsurface functions; iii) due to the error propagations, the final predicted maps of the subsurface soil properties captured from 1/3 to 2/3 of the total variance with a significantly decreasing performance with depth; and iv) these maps brought significant improvements over the existing soil maps of the region and showed soil patterns that largely agreed with the local pedological knowledge. This paper demonstrates the added value of combining modern remote sensing techniques with old legacy soil databases. (C) 2013 Elsevier B.V. All rights reserved.