Soil sodicity, generally measured through the exchangeable sodium percentage (ESP), is the most prevalent edaphic stress in Australia, particularly in the northern grains-growing region (NGR) of Australia. Farmers, scientists and policy-makers need accurate maps of sodicity for environmental modelling, and to make rational decisions for management. A common approach to mapping soil properties for multiple target depths is to first harmonize soil profile data to the target depth intervals, usually through equal-area spline functions, then undertake a two-dimensional spatial analysis (often regression kriging) for each interval; we refer to this as a ‘spline-then-krige’ (STK) approach. An alternative is to calibrate a single three-dimensional model that describes soil variation both horizontally and vertically, as is done in the increment-averaged kriging (IAK) approach. This study compares IAK and STK for mapping ESP in the NGR. A machine-learning algorithm was used to model trends based on environmental covariates in both approaches. Both IAK and STK captured the ESP variation reasonably well (concordance correlation coefficient values for IAK around 0.5 for all depths), with IAK giving slightly better predictive accuracies for all depths, most evidently in the topsoil. IAK allows the natural use of all available data, explicitly accounts for the sample support of the soil data and the horizontal and vertical auto-correlation between and within soil profiles, and enables the propagation of uncertainty from the original data through to the final predictions and prediction uncertainties. We conclude that IAK can provide a competitive alternative to STK approaches for mapping soil properties in three dimensions.