Abstract The growing human population and demand for food have significantly impacted soil resources. Understanding the spatiotemporal change of soil conditions is important to support food production, environmental sustainability, and climate change adaptation. Nevertheless, spatiotemporal prediction of soil properties could be seriously influenced by the uncertainties of the data and model. Integrated Nested Laplace Approximation (INLA) with the Stochastic Partial Differential Equation (SPDE) was proposed as a general model that can account for the uncertainties in spatiotemporal soil modelling and prediction. INLA-SPDE has significant advantages in computation efficiency over commonly-used geostatistical methods with Markov Chain Monte Carlo. However, until now, only few pedometrics studies used it for soil spatial modelling. This study demonstrates an application of INLA-SPDE within a hierarchical spatiotemporal model for soil organic matter based on soil survey data collected in Jiangsu, China, during three periods, i.e., 1979–1982, 2000 and 2006–2007. Compared with updating digital soil maps using the Bayesian Maximum Entropy approach, the prediction generated using INLA-SPDE is more accurate. For example, the root mean square error using INLA-SPDE (i.e., 6.57 g kg−1) was reduced by 20% compared to the updating approach (i.e., 8.39 g kg−1). Moreover, accounting for sources of uncertainties made the prediction using INLA-SPDE more certain. Nevertheless, the uncertainty in the temporal prediction of soil change is still large due to the scarcity of data across the sampling periods. The INLA-SPDE model predicts much detailed spatiotemporal changes along the sampling periods. Therefore, this study recommends the use of INLA-SPDE within a hierarchical model as an effective method for studying spatiotemporal soil change.