Three-dimensional digital soil mapping (3D-DSM) quantifies both the horizontal and the vertical variability of soil properties. Most current studies in 3D-DSM were based on either one-dimensional profile depth functions or two-dimensional horizontal interpolation techniques, which did not allow true 3D visualization of spatial soil heterogeneity. Only a few studies have utilized the 3D variograms for mapping. Recent advances in proximal soil sensing technologies allow measurement and prediction of soil properties rapidly at multiple depths which could serve as input data for DSM. Various soil physical and chemical properties have already shown either direct or indirect relationships with the proximal soil sensing data. This study aims to test the methodology of 3D-DSM by incorporating a 3D regression kriging (RK) with multiple proximal soil sensing techniques. In this study, vis-NIR spectra were collected in-situ at 148 locations to about 1-m depth using the Veris® P4000 soil profiler at Field 26 of Macdonald Farm, McGill University. Additionally, 32 soil cores were collected out of the 148 locations to 1-m maximum depth and sectioned at 10-cm depth intervals for laboratory analysis of volumetric water content (VWC), soil organic matter (SOM), and clay content. Cubist spectral models were developed for each soil property at the 32 locations and then predicted to the 148 locations, which were then randomly split into calibration (70%, 103 locations) and validation (30%, 45 locations) datasets for mapping. The 3D-RK method included a trend prediction between calibration dataset and environmental covariates (including apparent soil electrical conductivity, gamma-ray radiation, and elevation) and a residual kriging. The generalized linear model (GLM), regression tree (RT), and random forest (RF) models were compared for trend prediction. The covariates were also simulated 100 times using sequential Gaussian simulations to fit into 3D-RK and calculate model uncertainty. As a result, complete 3D digital soil maps with uncertainty were developed. We found that the RF model outperformed GLM and RT in regard to interpreting non-linear soil-landscape relationships and resulting in marginally higher validation accuracy and smaller prediction uncertainty for VWC and clay. The GLM model resulted in slightly better validation results and smaller model uncertainty for SOM only. SOM and clay showed large horizontal and vertical variability and affected the spatial distribution of VWC. The validation accuracy was higher in the soil surface for most soil properties due to the uniform environment in the plow layer and sufficient environmental covariates collected at the soil surface. The mapping uncertainty increased with depth for VWC and clay content but decreased with depth for SOM because SOM content decreases with depth.