Understanding the spatial distribution of soil organic carbon (SOC) content over different climatic regions will enhance our knowledge of carbon gains and losses due to climatic change. However, little is known about the SOC content in the contrasting arid and sub-humid regions of Iran, whose complex SOC–landscape relationships pose a challenge to spatial analysis. Machine learning (ML) models with a digital soil mapping framework can solve such complex relationships. Current research focusses on ensemble ML models to increase the accuracy of prediction. The usual ensemble method is boosting or weighted averaging. This study proposes a novel ensemble technique: the stacking of multiple ML models through a meta-learning model. In addition, we tested the ensemble through rescanning the covariate space to maximize the prediction accuracy. We first applied six state-of-the-art ML models (i.e., Cubist, random forests (RF), extreme gradient boosting (XGBoost), classical artificial neural network models (ANN), neural network ensemble based on model averaging (AvNNet), and deep learning neural networks (DNN)) to predict and map the spatial distribution of SOC content at six soil depth intervals for both regions. In addition, the stacking of multiple ML models through a meta-learning model with/without rescanning the covariate space were tested and applied to maximize the prediction accuracy. Out of six ML models, the DNN resulted in the best modeling accuracies, followed by RF, XGBoost, AvNNet, ANN, and Cubist. Importantly, the stacking of models indicated a significant improvement in the prediction of SOC content, especially when combined with rescanning the covariate space. For instance, the RMSE values for SOC content prediction of the upper 0–5 cm of the soil profiles of the arid site and the sub-humid site by the proposed stacking approaches were 17% and 9% respectively, less than that obtained by the DNN models—the best individual model. This indicates that rescanning the original covariate space by a meta-learning model can extract more information and improve the SOC content prediction accuracy. Overall, our results suggest that the stacking of diverse sets of models could be used to more accurately estimate the spatial distribution of SOC content in different climatic regions.

Soil organic C (SOC) and soil moisture (SM) affect the agricultural productivity of soils. For sustainable food production, knowledge of the horizontal as well as vertical variability of SOC and SM at field scale is crucial. Machine learning models using depth‐related data from multiple electromagnetic induction (EMI) sensors and a gamma‐ray spectrometer can provide insights into this variability of SOC and SM. In this work, we applied weighted conditioned Latin hypercube sampling to calculate 25 representative soil profile locations based on geophysical measurements on the surveyed agricultural field, for sampling and modeling. Ten additional random profiles were used for independent model validation. Soil samples were taken from four equal depth increments of 15 cm each. These were used to approximate polynomial and exponential functions to reproduce the vertical trends of SOC and SM as soil depth functions. We modeled the function coefficients of the soil depth functions spatially with Cubist and random forests with the geophysical measurements as environmental covariates. The spatial prediction of the depth functions provides three‐dimensional (3D) maps of the field scale. The main findings are (a) the 3D models of SOC and SM had low errors; (b) the polynomial function provided better results than the exponential function, as the vertical trends of SOC and SM did not decrease uniformly; and (c) the spatial prediction of SOC and SM with Cubist provided slightly lower error than with random forests. Hence, we recommend modeling the second‐degree polynomial with Cubist for 3D prediction of SOC and SM at field scale.

As limited resources, soils are the largest terrestrial sinks of organic carbon. In this respect, 3D modelling of soil organic carbon (SOC) offers substantial improvements in the understanding and assessment of the spatial distribution of SOC stocks. Previous three-dimensional SOC modelling approaches usually averaged each depth increment for multi-layer two-dimensional predictions. Therefore, these models are limited in their vertical resolution and thus in the interpretability of the soil as a volume as well as in the accuracy of the SOC stock predictions. So far, only few approaches used spatially modelled depth functions for SOC predictions. This study implemented and evaluated an approach that compared polynomial, logarithmic and exponential depth functions using non-linear machine learning techniques, i.e. multivariate adaptive regression splines, random forests and support vector machines to quantify SOC stocks spatially and depth-related in the context of biodiversity and ecosystem functioning research. The legacy datasets used for modelling include profile data for SOC and bulk density (BD), sampled at five depth increments (0-5, 5-10, 10-20, 20-30, 30-50 cm). The samples were taken in an experimental forest in the Chinese subtropics as part of the biodiversity and ecosystem functioning (BEF) China experiment. Here we compared the depth functions by means of the results of the different machine learning approaches obtained based on multi-layer 2D models as well as 3D models. The main findings were (i) that 3^rd degree polynomials provided the best results for SOC and BD (R² = 0.99 and R² = 0.98; RMSE = 0.36% and 0.07 g cm^-3). However, they did not adequately describe the general asymptotic trend of SOC and BD. In this respect the exponential (SOC: R² = 0.94; RMSE = 0.56%) and logarithmic (BD: R² = 84; RMSE = 0.21 g cm^-3) functions provided more reliable estimates. (ii) random forests with the exponential function for SOC correlated better with the corresponding 2.5D predictions (R²: 0.96 to 0.75), compared to the 3^rd degree polynomials (R²: 0.89 to 0.15) which support vector machines fitted best. We recommend not to use polynomial functions with sparsely sampled profiles, as they have many turning points and tend to overfit the data on a given profile. This may limit the spatial prediction capacities. Instead, less adaptive functions with a higher degree of generalisation such as exponential and logarithmic functions should be used to spatially map sparse vertical soil profile datasets. We conclude that spatial prediction of SOC using exponential depth functions, in conjunction with random forests is well suited for 3D SOC stock modelling, and provides much finer vertical resolutions compared to 2.5D approaches.

The use of machine learning as a method for knowledge discovery is often critically discussed in soil science and related environmental disciplines. Reviews of the use of machine learning in digital soil mapping identified few studies that incorporated existing soil knowledge of transformation and translocation processes in soils and mechanistic relationships between covariates in the modelling process. Even models trained with predictors that are meaningless from a soil science perspective can have high accuracies. To test and widen this perspective, we expanded the setup of a previous study by Wadoux, Samuel‐Rosa, et al. (2020) from one hypothetical case study to a larger set of 668 hypothetical case studies in 334 study areas. We found that the high accuracy of one single model for a specific area was part of a wide range of possible accuracy metrics (concordance correlation coefficient: 0.16–0.91) when applying the same set of meaningless predictors to all study areas. We discuss these spurious correlations in the context of explainable machine learning and highlight how the important elements of model explainability, model input and model output largely depend on discipline‐specific domain knowledge. As soil science knowledge is often incorporated implicitly, we argue that the motivation behind covariate selection should be discussed more explicitly to achieve soil science knowledge beyond spatial prediction.

Multi-scale contextual modelling is an important toolset for environmental mapping. It accounts for spatial dependence by using covariates on multiple spatial scales and incorporates spatial context and structural dependence to environmental properties into machine learning models. For spatial soil modelling, three relevant scales or ranges of scale exist: quasi-local soil formation processes that are independent of the spatial context, short-range catenary processes, and long-range processes related to climate and large-scale terrain settings. Recent studies investigated the spatial dependence of topsoil properties only. We hypothesize that soil properties within a soil profile were formed due to specific interactions between different features and scales of the spatial context, and that there are depth gradients in spatial and structural dependencies. The results showed that for topsoil, features at small to intermediate scales do not increase model accuracy, whereas large scales increase model accuracy. In contrast, subsoil models benefit from all scales—small, intermediate, and large. Based on the differences in relevance, we conclude that the relevant ranges of scales do not only differ in the horizontal domain, but also in the vertical domain across the soil profile. This clearly demonstrates the impact of contextual spatial modelling on 3D soil mapping.

We present the first detailed soil property maps at multiple depths for the northwestern autonomous Kurdistan region of Iraq (Dohuk). A total of 532 soil samples from 122 sites were collected at five depth increments (0–10, 10–30, 30–50, 50–70, and 70–100 cm), and their mid-infrared (MIR) spectra were measured. A subset of 108 samples, selected via Kennard–Stone sampling, was analysed in a laboratory on ten soil properties. A Cubist model was trained and used from these measured values to predict all samples' soil properties from their MIR spectra. Digital soil mapping was conducted using a machine learning regression techniques based on a quantile random forest model, trained on the predicted soil properties and using a total of 85 covariates at 30 m pixel resolution, resulting in 50 prediction maps in total. Results were compared with the SoilGrids 2.0 product and a regional texture model. Soil depth was also mapped using a similar model with 26 covariates. Our model outperformed global SoilGrids 2.0 predictions in resolution and accuracy, with texture RMSEs (sand: $\stackrel{\mathrm{‾}}{x}$ RMSE 11.03; silt: $\stackrel{\mathrm{‾}}{x}$ RMSE 8.82; clay: $\stackrel{\mathrm{‾}}{x}$ RMSE 7.39) comparable to local models. Key predictors included Landsat 8 SWIR, EVI, SAVI, Sentinel 2 SWIR, PET and solar radiation. Spatial patterns reflected the contrast between the flat areas of the Selevani and Zakho plains, as opposed to the shallower and steeper Little Khabur Valley and anticline formations. Furthermore, the soil depth prediction model (R² 0.39; RMSE 30.76 cm) showed strong correlation with slope and a similar pattern distribution with deeper soils in the flat areas of the Selevani and Zakho plains, while shallow soils were predicted in the anticline and strongly erodible areas. Our comprehensive dataset (https://doi.org/10.1594/PANGAEA.973700, Bellat et al., 2024 a ; https://doi.org/10.1594/PANGAEA.973701, Bellat et al., 2024 b ; https://doi.org/10.1594/PANGAEA.973714, Bellat et al., 2024 c ; https://doi.org/10.6084/m9.figshare.31320958.v2, Bellat et al., 2026 a ; https://doi.org/10.57754/FDAT.d5h1h-4x027, Bellat et al., 2026 b ) offers substantial insights for soil knowledge in the region, as well as for aridic and semi-aridic areas.