Digital soil mapping is a rapidly growing area of soil research that has great potential for advancing soil survey activities, and knowledge of soil-landscape relationships. To date many successful studies have shown that geospatial data can be successfully used as covariates to quantitatively model soil spatial variation. Most of these successful studies though have only focused on a select number of soil properties over relatively small areas. In addition to this, current methods of deriving terrain parameters from digital elevation models, the most heavily used soil covariate, are dependent on the spatial resolution of the terrain parameters, which can not account for the varying scales over which pedogenesis occurs. To address these issues our objectives were to develop predictive models of a variety of soil properties at the watershed scale, and assess the impact of varying the neighborhood size used to derive terrain parameters on the performance of these predictive models. Throughout the watershed of study, which is approximately 82,500 acres, 97 sites were sampled using a stratified-random design. At each site a soil pit was excavated and described to 140 cm or bedrock, with grab samples collected from each horizon. To generate the spatial models, linear and tree-based models were used, with terrain parameters, Landsat TM+ band ratios, National Land Cover Dataset, and geology as the soil covariates. Early results of this study have show that the distribution of surface and subsurface rock fragments can be successfully modeled, while other physical and morphological properties were unsuccessful. The modeled spatial distribution of rock fragments indicates a greater range of variation across the mountain hillslopes than the published soil survey, and may be used as a knowledge base to infer other soil-landscape relationships.