Abstract As any model for digital soil mapping suffers from different types of errors, including interpolation errors, it is important to quantify the uncertainty associated with the maps produced. The most common approach is some form of regression kriging or variation involving geostatistical simulation. Another way of assessing the spatial uncertainty lies in the Bayesian approach where the uncertainty is described by the posterior density. Typically Markov Chain Monte Carlo is used to compute the posterior density; however, this process is computationally intensive. The aim of this paper is to present an example of Bayesian uncertainty evaluation using (Bayesian) latent Gaussian models fitted using INLA (Integrated Nested Laplace Approximation) and with the SPDE (Stochastic Partial Differential Equation) approach for modelling the spatial correlation. For illustration, soil organic matter content in the Grampian region of Scotland (UK, about 12,100 km2) was modelled for topsoil (2D) and whole-profile data (3D). Results were assessed using in-sample and out-of-sample measures and compared for distribution similarity, variogram and spatial structure reproduction, computational load and uncertainty ranges. The results were also compared with outputs from an extension of scorpan-kriging. The Bayesian framework using INLA offers a viable alternative to existing methods for digital soil mapping, with comparable validation results, important computational gains, good assessment of uncertainty and potential for integrated modelling uncertainty propagation.