Visible near-infrared reflectance spectroscopy (VNIRS) has proven to be a rapid and inexpensive method for soil assessment and could be useful for mapping soil properties. Our objective was to use VNIRS to predict soil organic matter (OM) content and assess the use of VNIRS in combination with geostatistical methods to characterize spatial variability and map OM distribution. Soil samples were collected from up to 4 depth increments: L1 (0-30 cm), L2 (30-60 cm), L3 (60-120 cm), and L4 (120-180 cm) at 152 sites and analyzed for OM content by weight loss on ignition. Spectral calibration models were developed using linear (partial least squares) and nonlinear (regression trees) methods using calibration data. Model evaluation using independent validation data showed that regression tree prediction models were superior and therefore were used to characterize spatial variability and map OM in L1 and L2. Semivariograms derived from true and spectra-predicted OM had similar spatial structure, with nugget-to-sill ratio of 11% and 19% in L1, and 47% and 53% in L2, respectively. Ordinary kriging prediction maps derived from true and spectra-predicted data showed similar spatial patterns of high and low OM values across the study area. Spectra-predicted data mapped 70.9% and 94.5% of study area within ±0.50% of the predictions derived from true data in L1 and L2, respectively. Our results demonstrate the potential of VNIRS for mapping soil OM and indicate prospects to enhance the scope of soil spectroscopy to digital soil mapping.