# Mapping and understanding the regional farmland SOC distribution in southern China using a Bayesian spatial model

Bifeng Hu, Yibo Geng, Hanjie Ni, Zhou Shi, Zheng Wang, Nan Wang, Jipeng Luo, Modian Xie, Qian Zou, Thomas Optiz, et al.

► To cite this version:

Bifeng Hu, Yibo Geng, Hanjie Ni, Zhou Shi, Zheng Wang, et al. Mapping and understanding the regional farmland SOC distribution in southern China using a Bayesian spatial model. Geoderma, 2025, 460, pp.117446. 10.1016/j.geoderma.2025.117446. hal-05215815

HAL Id: hal-05215815

https://hal.inrae.fr/hal-05215815v1

Submitted on 20 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Mapping and understanding the regional farmland SOC distribution in southern China using a Bayesian spatial model

Bifeng Hu<sup>a,b,●</sup>, Yibo Geng<sup>a,●</sup>, Hanjie Ni<sup>a</sup>, Zhou Shi<sup>c,●</sup>, Zheng Wang<sup>c</sup>, Nan Wang<sup>d</sup>, Jipeng Luo<sup>c</sup>, Modian Xie<sup>f</sup>, Qian Zou<sup>a</sup>, Thomas Optiz<sup>g</sup>, Hongyi Li<sup>a,b,∗</sup>

## ARTICLE INFO

Handling Editor: Brandon Heung

**Keywords:**

Soil organic carbon

Farmland

INLA-SPDE

Interpretable machine learning model

## ABSTRACT

Information on the spatial distribution of soil organic carbon (SOC) in regional farmland is crucial for improving management and production. Mapping SOC in farmlands is challenging due to the strong variation of SOC caused by the influence of natural and anthropogenic activities. Additionally, currently widely used predictive models usually suffer from a lack of model interpretability. To fill these gaps, here we use a Bayesian spatial model termed Integrated Nested Laplace Approximation with the Stochastic Partial Differential Equation (INLA-SPDE) to produce the fine scale SOC map in the farmland of Jiangxi Province, south China based on an extensive soil survey dataset (n = 16,050). The competitive adaptive reweighted sampling algorithm – partial least square (CARS-PLS) algorithm is adopted to select the most related covariates from the original covariates pool. Then, the performance of Random Forest (RF), Geographically Weighted Regression (GWR), and Ordinary Kriging (OK) was compared with INLA-SPDE. Finally, an interpretable machine learning model, the SHapley Additive exPlanation (SHAP), is used to quantify the environmental covariates' contribution to mapping SOC, as well as mapping spatial varying primary covariates for predicting SOC in the study area. We find that INLA-SPDE was able to handle a large data and performed much better than OK and GWR with an improvement of 38.89 % and 117.39 % in R², respectively. It also outperforms RF. Overall, amount of straw return, mean annual precipitation, mean annual solar radiation are the most important covariates for mapping SOC. Locally, soil management are the most important covariates for mapping SOC in 50.52 % regions of the study area, followed by climate factors (22.06 %), soil properties (17.09 %), terrain (6.38 %), lithology (2.21 %) and biota factors (1.72 %). Our study demonstrates the advantages of INLA-SPDE on mapping SOC compared with geostatistical and RF for SOC mapping and provides valuable implications for interpreting the results of digital soil mapping.

## 1. Introduction

Soil is the largest reservoir of organic carbon, acting as both a major sink and source of carbon, and stores three times more carbon than the atmosphere (Bai et al., 2023; Batjes, 1996; Hu et al., 2023a; Ma et al., 2023; Post and Kwon, 2000; Sun et al., 2021). Even small increases or decreases in soil organic carbon (SOC), which result from interactions

among edaphic, biotic, and abiotic factors, can significantly influence atmospheric CO₂ concentrations and intensify global climate change (Neeraj et al., 2025; Smith, 2007). In addition, SOC is the most important indicator of soil fertility and plays a vital role in soil management, environmental protection, and climate smart agriculture (Dai et al., 2024; Hu et al., 2024; Lal, 2004a; Zhang et al., 2024). Therefore, predicting and understanding the spatial distribution of SOC is crucial for

<footnote><sup>a</sup> Department of Land Resource Management, School of Public Administration, Jiangxi University of Finance and Economics, Nanchang 330013, China</footnote>

<footnote><sup>b</sup> Key Laboratory of Data Science in Finance and Economics of Jiangxi Province, Jiangxi University of Finance and Economics, Nanchang 330013, China</footnote>

<footnote><sup>c</sup> Institute of Agricultural Remote Sensing and Information Technology Application, Zhejiang University, Hangzhou 310058, China</footnote>

<footnote><sup>d</sup> Department of Earth System Science, Ministry of Education Key Laboratory for Earth System Modeling, Institute for Global Change Studies, Tsinghua University, 100084 Beijing, China</footnote>

<footnote><sup>e</sup> Department of Biology, Indiana University, Bloomington, IN, USA</footnote>

<footnote><sup>f</sup> School of Information Management and Mathematics, Jiangxi University of Finance and Economics, Nanchang 330013, China</footnote>

<footnote><sup>g</sup> INRAE, Biostatistics and Spatial Processes, 228 route de l'Aérodrome, Avignon 84914, France</footnote>

<footnote>* Corresponding author at: Department of Land Resource Management, School of Public Administration, Jiangxi University of Finance and Economics, Nanchang 330013, China.</footnote>

<footnote>E-mail address: lihongyi1981@zju.edu.cn (H. Li).</footnote>

<footnote>https://doi.org/10.1016/j.geoderma.2025.117446</footnote>

<footnote>Received 23 September 2024; Received in revised form 8 July 2025; Accepted 14 July 2025</footnote>

<footnote>Available online 17 July 2025</footnote>

<footnote>0016-7061/© 2025 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).</footnote>
broad issues, including optimizing farmland management, improving food production, and regulating climate change (Batjes, 1996; Chen et al., 2022; Georgiou et al., 2022; Hu et al., 2024; Li et al., 2025).

Traditionally, information on SOC is mainly obtained through field surveys and laboratory analysis, which has high accuracy but low efficiency and high cost. This made it is very hard to monitor in large areas (Yang et al., 2022; Yu et al., 2024; Gao et al., 2024). In the past decades, different approaches have been developed to predict spatial distribution of SOC and other soil properties across landscapes at different spatial scales, including linear regression (Wu et al., 2009), GWR (Brunsdon et al., 1996; Hu et al., 2023b; Jiang et al., 2025), OK (Cressie, 1988; Zare et al., 2021), and machine learning methods (Breiman et al., 2017; Hu et al., 2024; Mahjenabadi et al., 2022; Wadoux, 2019). Among these methods, OK is the most widely used geostatistical method, utilizing sample point data through intensive field sampling and making predictions based on a semivariance function (Emadi and Baghernejad, 2014). However, OK requires significant computation resources and its accuracy is heavily influenced by sampling density and spatial coverage (Shen et al., 2019; Ye et al., 2017). Additionally, OK does not incorporate environmental factors closely related to SOC, limiting its application in complex landscapes (McBratney et al., 2003; Hu et al., 2022). In contrast, GWR, an extension of the traditional regression framework, considers both the spatial heterogeneity of soil properties and the spatial dependence of influencing factors (Keser et al., 2012; Kumar et al., 2012). GWR used a weighting matrix where observations near the sample point are weighted more heavily, enhancing the localized regression coefficient's accuracy (Kupfer and Farris, 2006). However, there is a non-linear relationship between the soil property and environmental factors, and machine learning methods could deal with it well (Petersen et al., 2008). Among machine learning models, tree-based algorithms like random forest (RF) have shown the best performance in predicting soil properties (Gomes et al., 2019; Hengl et al., 2015; Hu et al., 2024; Khaledian et al., 2020; Silatsa et al., 2020). Nevertheless, RF could not integrate spatial autocorrelation into the prediction progress, which limits its predictive performance (Hu et al., 2020a; Keskin et al., 2019; Odeh et al., 1995; Szatmári et al., 2019). Based on this, recent studies demonstrate the effectiveness of mixed model approaches in resolving complex soil spatial patterns. Matinfar et al. (2021) demonstrated the applicability of a hybrid RF-OK and QRF-OK framework to field SOC mapping, achieving an acceptable level of accuracy. Molla et al., (2023) who combined empirical Bayesian kriging regression prediction with RF to predict the potentially toxic elements of greenfield soils, their work not only improved prediction accuracy, but also produced more robust estimates of uncertainty. Nevertheless, current researches generally suffer from a lack of ability to explain prediction results (Hu et al., 2025; Wadoux et al., 2021; Wadoux and Molnar, 2022). Hence, a method which overcomes the shortcomings of the above methods is urgently needed.

Integrated nested Laplace approximation and stochastic partial differential equation method (INLA-SPDE) is an advanced Bayesian spatial model based on Bayesian theory (Hu, 2020; Huang et al., 2017; Lom-bardo et al., 2018; Rue et al., 2009; Ryan et al., 2015). The INLA-SPDE has proved to be a flexible, accurate and computation-efficient method for constructing spatial and spatiotemporal models (Huang et al., 2017; Hu, 2020; Hu et al., 2020b; Simpson et al., 2023; Yadav et al., 2022). And it can output the posterior distribution of the model hyper-parameters, which makes the results more explainable (Minasny et al., 2011; Sun et al., 2021). Currently, it has been widely applied in various fields, including ambient air pollution (Wright et al., 2021; Zhong et al., 2025), mapping geochemical patterns (Wang and Zuo, 2021), and the spatiotemporal characteristics of disease transmission (Moraga et al., 2021; Subramani and Dhakshnamoorthy, 2025). Especially, Poggio et al. (2016) first applied the INLA-SPDE model to digital soil mapping (DSM) and verified its feasibility for soil organic matter (SOM) mapping in Scotland's Grampian region, highlighting its uncertainty assessment advantage. Huang et al. (2017) demonstrated that INLA-SPDE is as

robust as the residual maximum likelihood method in mapping SOC, Ti/ Zr ratio, and pH at Nowley Farm, Australia. Sun et al. (2021) further applied INLA-SPDE to model SOM's spatiotemporal variation at the regional scale. Giannini-Kurina et al. (2022) demonstrated the effectiveness of INLA-SPDE in mapping atrazine persistence in soils of central Argentina, where it provided more accurate predictions compared to other traditional geostatistical methods. Yang et al. (2023) found that INLA-SPDE outperformed RF in SOM mapping. The INLA-SPDE's flexibility and robust uncertainty quantification make it particularly valuable for DSM applications that require reliable prediction interval (Minasny et al., 2011; Poggio et al., 2016; Sun et al., 2021). However, the performance and feasibility of INLA-SPDE model for fine scale mapping SOC in regions with complex terrain conditions like hilly regions in southern China still need to be verified. Therefore, the main objectives of this study are: (1) to compare the performance of INLA-SPDE to classic different spatial predictive models on mapping SOC in farmlands of Jiangxi Province in southern China which is a hilly regions featured by very complex terrain condition; (2) to analyze uncertainty of the results produced by INLA-SPDE and RF models; (3) to quantify the overall importance of different covariates and map the spatially varying primary covariates to map SOC at the pixel scale by using INLA-SPDE and an explainable machine learning model termed Shapley Additive explanations (SHAP).

## 2. Materials and methods

### 2.1. Study area

Our study area covers all farmlands across Jiangxi Province, located on the south bank of the middle reaches of the Yangtze River in southern China, with altitudes ranging from 191 m below sea level to 2158 m above sea level (Fig. 1). Jiangxi Province covers 166,900 km², with mountains and hills being the main landforms (Hu et al., 2023). In 2023, Jiangxi Province had 27,216 km² of farmland and produced 21.98 million tons of grain, making it one of China's most important grain-producing areas. The annual average temperature in Jiangxi Province ranges from 16.3 to 19.5°C, and annual precipitation ranges from 1241 to 1943 mm (Hu et al., 2024). According to the Genetic Soil Classification of China, the main soil types in Jiangxi Province are red soil (classified as Acrisol in the WRB system) and paddy soil (classified as Fluvisol in the WRB system).

### 2.2. Data collection

From 2012 to 2013, as part of the Ministry of Agriculture and Rural Affairs' soil testing and formulated fertilization program, 16,050 soil samples were collected from the surface layer (0–20 cm) of farmland in Jiangxi Province. Sampling points were selected based on soil types, topographic characteristics, and transportation convenience. Additionally, land use and spatial distribution were considered. Each soil sample was a composite of five subsamples collected from adjacent areas within a 5-meter radius, resulting in approximately 1 kg of soil after removing surface debris. A portable GPS device was used to record precise sampling coordinates. Stones and tree roots were removed. All soil samples were packed into plastic bags, transported to the laboratory, air-dried, crushed with a mortar and pestle, and sieved through a 2 mm mesh. SOM was evaluated using the potassium dichromate oxidation method (Dong et al., 2019). The organic matter was then converted into the weight percentage of SOC using the Van Bemmelen conversion ratio (0.58) (Lal, 2004b; Ma et al., 2017; Minasny et al., 2020).

### 2.3. Environmental covariates

#### 2.3.1. Covariates collection

Various factors related to soil properties, climate, terrain, biota, and soil management were collected to model SOC in Jiangxi Province's
<figure><img src="image_3.png" /><figcaption>Fig. 1. The location of study area (a), Digital elevation model (DEM) (b) and the distribution of soil sampling sites (c).</figcaption></figure>

farmland soil (Table 1). Digital elevation data with a 30 m resolution was obtained from the European Space Agency (https://eop-cfi.esa.int/index.php/docs-and-mission-data/dem). Subsequently, SAGA GIS (https://saga-gis.org/) software was used to calculate the Topographic Wetness Index (TWI), Topographic Position Index (TPI), Multi-resolution Valley Bottom Flatness (MRVBF), Slope, Aspect, Valley Depth (VD), Slope Length (SL), Catchment Slope (CAS), Convergence Index (CI), Length of Slope Factor (LSF), and Catchment Area (CA).

To represent the climate factor in soil formation, data on mean annual temperature (MAT), mean annual precipitation (MAP), daytime land surface temperature (LSTD), nighttime land surface temperature (LSTN), evapotranspiration (ET), and mean annual solar radiation (MASR) were obtained from the Resources and Environmental Sciences Data Center, Chinese Academy of Sciences (RESDC) (https://www.redc.cn) at a spatial resolution of 1 km. The mean annual solar radiation (MASR) was computed based on the yearly averages for 2012 and 2013 (https://www.geodata.cn).

Soil class maps correspond to the Genetic Soil Classification of China (GSCC), and parent material maps were downloaded from the National Earth System Science Data Center (https://www.geodata.cn). Maps depicting coarse fragments (>2 mm) (CF), clay, sand, and silt were obtained from SoilGrids250 2.0 (https://soilgrids.org/). Additionally, a cropping intensity map by Zhang et al. (2021) with a spatial resolution of 30 m was included in the study.

### 2.3.2. Covariates preprocessing

Using ArcGIS 10.8, we employed the OK technique to generate maps of pH, bulk density (BD), cation exchange capacity (CEC), total nitrogen (TN), total phosphorus (TP), total potassium (TK), exchangeable magnesium (EMg), available silicon (ASi), available sulphur (AS), and plough layer depth (PLD) for the entire survey area. Soil management data were linked to corresponding polygons using the spatial join tool in ArcGIS 10.8. These polygons were then converted into raster layers

using ArcGIS 10.8. All original environmental covariate raster layers were re-projected to the Universal Transverse Mercator reference system (WGS 1984 UTM 49 N) and standardized to a grid resolution of 30 m using the nearest neighbor resampling method. Finally, all raster data were clipped to the cropland boundaries of Jiangxi Province.

### 2.3.3. Covariates selection

To select the optimal covariate for predicting SOC, we used CARS-PLS. The CARS-PLS algorithm can show advantages in feature selection(Li et al., 2009). This method realizes dynamic variable screening through Monte Carlo sampling and exponential decay strategy, and can effectively complete data dimensionality reduction and optimize computational efficiency (Wu et al., 2023; Zeng et al., 2023). Furthermore, its adaptive reweighting strategy does not require a preset threshold. Compared with the variance inflation factor, it can effectively reduce subjective intervention. Meanwhile, it relies on cross-validation mechanisms (such as RMSECV index optimization) to ensure the predictive efficacy of variable combinations and the robustness of the model. It is worth noting that this method is applicable to large-scale datasets and can significantly alleviate the risk of overfitting (Hu et al., 2024). CARS-PLS involves four steps: (1) the 80 % of soil samples from the calibration set for the modeling set by the Monte Carlo method, while the remaining samples (20 % of calibration) are used as the prediction set for modeling by PLS. (2) The exponentially decreasing function (EDF) removes wavelengths with small absolute regression coefficient weights. (3) At each sampling, adaptive weighted sampling (ARS) selects wavelength variables from the previous sampling, PLS modeling is performed, and RMSECV is calculated. (4) After N samplings, the wavelength variable subset with the minimum RMSECV value is selected as the characteristic wavelength (Li et al., 2009).
<table><caption>**Table 1**
List of covariates collected (n = 44).</caption><thead><tr><th>Environmental factor</th><th>Variable</th><th>Acronym</th><th>Spatial resolution</th><th>Unit</th><th>Source</th></tr></thead><tbody><tr><td rowspan="5">Climate</td><td>Mean annual temperature</td><td>MAT</td><td>1 km</td><td>°C</td><td>https://www.resdc.cn/</td></tr><tr><td>Mean annual precipitation</td><td>MAP</td><td>1 km</td><td>mm</td><td>https://www.resdc.cn/</td></tr><tr><td>Evapotranspiration</td><td>ET</td><td>1 km</td><td>–</td><td>https://lpdaac.usgs.gov/</td></tr><tr><td>Mean annual solar radiation</td><td>MASR</td><td>1 km</td><td>MJ/m²</td><td>https://www.geodata.cn</td></tr><tr><td>Land surface temperature, day time</td><td>LSTD</td><td>1 km</td><td>°C</td><td>https://lpdaac.usgs.gov/</td></tr><tr><td rowspan="10">Terrain</td><td>Land surface temperature, night time</td><td>LSTN</td><td>1 km</td><td>°C</td><td>https://lpdaac.usgs.gov/</td></tr><tr><td>Elevation</td><td>Ele</td><td>30 m</td><td>m</td><td>USGS ASTGTM</td></tr><tr><td>Topographic wetness index</td><td>TWI</td><td>30 m</td><td>–</td><td>Calculated from Elevation</td></tr><tr><td>Topographic position index</td><td>TPI</td><td>30 m</td><td>–</td><td>Calculated from Elevation</td></tr><tr><td>Multi-resolution valley bottom flatness</td><td>MRVBF</td><td>30 m</td><td>–</td><td>Calculated from Elevation</td></tr><tr><td>Slope</td><td>Slope</td><td>30 m</td><td>°</td><td>Calculated from Elevation</td></tr><tr><td>Aspect</td><td>Aspect</td><td>30 m</td><td>°</td><td>Calculated from Elevation</td></tr><tr><td>Valley depth</td><td>VD</td><td>30 m</td><td>m</td><td>Calculated from Elevation</td></tr><tr><td>Slope length</td><td>SL</td><td>30 m</td><td>m</td><td>Calculated from Elevation</td></tr><tr><td>Catchment Slope</td><td>CAS</td><td>30 m</td><td>°</td><td>Calculated from Elevation</td></tr><tr><td rowspan="10">Soil properties</td><td>Convergence Index</td><td>CI</td><td>30 m</td><td>–</td><td>Calculated from Elevation</td></tr><tr><td>Length-Slope Factor</td><td>LSF</td><td>30 m</td><td>–</td><td>Calculated from Elevation</td></tr><tr><td>Catchment Area</td><td>CA</td><td>30 m</td><td>–</td><td>Calculated from Elevation</td></tr><tr><td>Bulk density</td><td>BD</td><td>30 m</td><td>–</td><td>This study</td></tr><tr><td>Cation exchange capacity</td><td>CEC</td><td>30 m</td><td>cmol/kg(+)</td><td>This study</td></tr><tr><td>Soil pH</td><td>pH</td><td>30 m</td><td>–</td><td>This study</td></tr><tr><td>Total nitrogen</td><td>TN</td><td>30 m</td><td>g/kg</td><td>This study</td></tr><tr><td>Total phosphorus</td><td>TP</td><td>30 m</td><td>g/kg</td><td>This study</td></tr><tr><td>Total potassium</td><td>TK</td><td>30 m</td><td>g/kg</td><td>This study</td></tr><tr><td>Exchangeable Magnesium</td><td>EMg</td><td>30 m</td><td>cmol/kg(+)</td><td>This study</td></tr><tr><td>Available silicon</td><td>ASi</td><td>30 m</td><td>mg/kg</td><td>This study</td></tr><tr><td>Available sulphur</td><td>AS</td><td>30 m</td><td>mg/kg</td><td>This study</td></tr><tr><td>Plough layer depth</td><td>PLD</td><td>30 m</td><td>cm</td><td>This study</td></tr><tr><td rowspan="4">Biota</td><td>Clay</td><td>Clay</td><td>250 m</td><td>–</td><td>SoilGrids250m 2.0</td></tr><tr><td>Sand</td><td>Sand</td><td>250 m</td><td>–</td><td>SoilGrids250m 2.0</td></tr><tr><td>Silt</td><td>Silt</td><td>250 m</td><td>–</td><td>SoilGrids250m 2.0</td></tr><tr><td>Percentage of soil coarse fragment (&gt; 2 mm) (%)</td><td>CF</td><td>250 m</td><td>–</td><td>SoilGrids250m 2.0</td></tr><tr><td rowspan="4">Lithology</td><td>Normalized Difference Vegetation Index</td><td>NDVI</td><td>1 km</td><td>–</td><td>https://www.resdc.cn/</td></tr><tr><td>Net primary production</td><td>NPP</td><td>1 km</td><td>–</td><td>https://www.resdc.cn/</td></tr><tr><td>Rotation system</td><td>RS</td><td>–</td><td>–</td><td>This study</td></tr><tr><td>Soil class</td><td>SC</td><td>1 km</td><td>–</td><td>The Second National Soil Survey</td></tr><tr><td rowspan="7">Soil management</td><td>Parental material</td><td>PM</td><td>1 km</td><td>–</td><td>The Second National Soil Survey</td></tr><tr><td>Input of nitrogen fertiliser</td><td>NF</td><td>–</td><td>kg/ha</td><td>This study</td></tr><tr><td>Input of phosphate fertiliser</td><td>PF</td><td>–</td><td>kg/ha</td><td>This study</td></tr><tr><td>Input of potash fertiliser</td><td>KF</td><td>–</td><td>kg/ha</td><td>This study</td></tr><tr><td>Amount of straw return</td><td>SRA</td><td>–</td><td>kg/ha</td><td>This study</td></tr><tr><td>Irrigation capacity</td><td>IC</td><td>–</td><td>–</td><td>This study</td></tr><tr><td>Drainage capacity</td><td>DC</td><td>–</td><td>–</td><td>This study</td></tr><tr><td></td><td>Soil erosion degree</td><td>SED</td><td>–</td><td>–</td><td>This study</td></tr></tbody></table>

## 2.4. INLA-SPDE model

The INLA-SPDE model was developed by Rue et al. (2009). In INLA-SPDE, the INLA allows joint estimation of all model components, including smoothness parameters for any sampling design and with large data size, such that estimation uncertainty is propagated through all model components jointly. Through the hierarchical formulation of INLA-based models, we can accommodate measurement errors at the observation layer. The SPDE approach allows using sparse precision matrices (inverse covariance matrices) in spatio-temporal Gaussian models, tractable in very high dimension and useful for defining flexible prior structures for space-varying components. In addition, the INLA-SPDE method uses to approximate the Matérn covariance structure in spatially continuous processes (Gacutan et al., 2024; He et al., 2024; Hu, 2020; Paige et al., 2022). This covariance function can describe the law of the autocorrelation intensity attenuation with distance (such as range parameters and variance parameters), so as to quantify the spatial dependence. For example, a short-range parameter indicates a stronger autocorrelation between closely adjacent sample points, while a long-range parameter reflects a wider spatial correlation.

In this study, using SOC as an example, if the indices of the SOC data model change continuously in a two-dimensional space $\mathbb{R}$, and the study area $D$ is a subset of this space, $D \subset \mathbb{R}$. Let points vary continuously

within $D$, and let $X(s)$ represent the SOC content at points, a random variable reflecting its spatial randomness. Suppose there are $n$ points within region $D$, denoted as $s_1, s_2, s_3...s_n$. Then, $X(s_1), X(s_2), X(s_3)...X(s_n)$ form a random vector of $n$ random variables. These random variables are spatially autocorrelated. Therefore, the set of random variables $X(s)$ at multiple locations $s$ within the study area $D$ can be called a spatial random field or random function. If the distribution of the random field is Gaussian, it is called a Gaussian random field (GRF) (Wu, 2021). If SOC satisfy or approximate a Gaussian distribution, a Gaussian random field can be defined by the mean and covariance functions. The covariance function used in the INLA-SPDE method is the Matérn model (Hu, 2020; Huang et al., 2017; Lindgren et al., 2011; Li et al., 2018):

$$ \operatorname{cov}(X(s_i), X(s_j)) = \sigma_e^2 \delta_{ij} + \frac{\sigma_x^2}{2^{v-1}\Gamma(v)} (\kappa s_i - s_j)^v K_v(\kappa s_i - s_j) \quad (1) $$

here, $s_i - s_j$ represents the Euclidean distance between two locations $s_i$ and $s_j$; $\sigma_x^2$ is the variance; $\kappa$ is the scale parameter; $K_v$ is the modified Bessel function of the second kind; $\Gamma(v)$ is the Gamma function (Huang et al., 2017; Minasny and McBratney, 2005); $\sigma_e^2$ is the noise variance (nugget effect); $\delta_{ij}$ is the Kronecker delta (Minasny and McBratney, 2016).

When a Gaussian field satisfies the Matérn covariance function
model, it can be regarded as a stationary solution of a SPDE. The stochastic partial differential equation is expressed as (Lindgren et al., 2011):

$$ (\kappa^2 - \Delta)^{\frac{\alpha}{2}} X(s) = W(s), s \in \mathbb{R} \qquad (2) $$

here, $W(s)$ represents white noise, and $\Delta$ denotes the Laplace operator: $\Delta = \frac{\partial^2}{ds_1^2} + \frac{\partial^2}{ds_2^2}$, where $s_1$ and $s_2$ are the coordinates in the two-dimensional space (Wu, 2021). The parameter $\alpha$ is a positive integer related to the smoothness parameter $v$, defined as $\alpha = v+d/2$ (Lindgren and Rue, 2015; Poggio et al., 2016). Based on previous studies, $\alpha$ is set to 2 (with $v=1$ and $d=2$) (Sun et al., 2021; Whittle, 1954). The SPDE can approximate the continuous spatial random variable $X(s)$ as a piecewise linear function, significantly reducing the computational burden (Blangiardo and Cameletti, 2015; Huang et al., 2017). For more details on the theory and implementation of INLA-SPDE, refer to Blangiardo et al. (2013), Martins et al. (2012), Poggio et al. (2016) and Rue et al. (2009).

## 2.5. Ordinary kriging

Based on regionalized variables and theoretical principles, OK is a geostatistical method that uses spatial correlation and variograms as primary tools (Matheron, 1971; Takoutsing and Martín, 2025; Wu et al., 2025). The conventional OK model is represented as follows:

$$ Z_{(x_0)} = \sum_{i=1}^{n} \lambda_i(x_0)Z(x_i) \qquad (3) $$

$Z_{(x_0)}$ denotes the predicted value of the variable at the location $x_0$, $Z(x_i)$ represents the observed data points, $\lambda_i(x_0)$ are the weights corresponding to the observed values, determined through the stationary OK system, and $n$ is the number of observed data points within a specified neighborhood.

## 2.6. Geographically weighted regression

The GWR is a local estimation technique that allows for location-specific variations in regression coefficients (Brunsdon et al., 1996; Fotheringham et al., 2002). Given a set of target soil properties and predictors, a traditional linear regression model fitted using ordinary least squares (OLS) can be represented as:

$$ Y(i) = \beta_0(u_i, v_i) + \sum_k \beta_k(u_i, v_i)x_{ik} + \varepsilon_i \qquad (4) $$

here, $(u_i, v_i)$ represent the coordinates of location $i$, $\beta_0(u_i, v_i)$ denotes the intercept, $\beta_k$ is the regression coefficient, and $x_{ik}$ is the environmental variable at location $i$, with $k$ indicating the number of environmental variables. The regression parameters are estimated at each specific location $(u_i, v_i)$. The parameters near point $i$ are determined using weighted least squares.

## 2.7. Random forest

The RF is a robust machine-learning algorithm based on an ensemble of decision trees and is widely used for predicting the spatial distribution of soil properties (Hu et al., 2024; Kuang et al., 2023). The RF model offers several advantages, including resilience to noise and reduced tendency to overfit (Breiman et al., 2017). The previous studies has been demonstrated that the RF model consistently outperforms alternative methods in predictive performance (Camera et al., 2017; Dai et al., 2024). In this study, the primary parameters of RF algorithm (the number of trees, ntree; the minimum count of terminal nodes, mtry) were optimized to 500 and 5 by 10-fold cross-validation (Breiman et al., 2017; Hu et al., 2020a). In addition, the CARS-PLS was used to select the

suitable covariates to obtain the best prediction results.

To assess the prediction uncertainty of RF, the 100 bootstrap-based realizations of RF predictions were achieved and the results were averaged outputs on each grid node. Then, the uncertainty was assessed through quantifying the total variance, which was calculated by summing the auxiliary program variance and the validation mean squared errors (Chen et al., 2022). For each pixel, the 95 % confidence interval (this study adopted the analysis of the median prediction) is the square root of the overall prediction variance multiplied by the z-score of 1.96, corresponding to the 95 % confidence interval (Zhou et al., 2022).

## 2.8. Shapley additive exPlanation

SHAP Values were originally developed in coalitional game theory (Wadoux et al., 2023). From a game theory perspective, a modeling exercise can be seen as the combination of multiple collaborative games. In each game, agents (explanatory variables) strategically interact to achieve the goal of predicting a single observation (Nash, 1950; Rasmussen, 1989). The RF model is interpreted using SHAP Values (Shapley, 1997; Hu et al., 2025a), as described and applied in soil science by Padarian et al. (2020).

To gain insight into how each feature contributes to the final prediction outcome, methods in this category are explanatory models that represent the relationship between features and predictions in the form of a linear function with simplified binary variables:

$$ f(x) = g(z) = \theta_0 + \sum_i \theta_i z_i^m \qquad (5) $$

where $f(x)$ is the original model; $g(x)$ is the explanation model; $z$ is the simplified input; $x = h_x(z)$ is a mapping function to the original method; and $\theta_0 = f(h_x(0))$ is the model output without all of the simplified inputs.

## 2.9. Evaluation of the predicted soil maps

The data set was randomly divided into a calibration group (11,238 samples) and a validation group (4812 samples) in a ratio of 70 % to 30 % (Barkhordari and Qi, 2025; Chen et al., 2021; de Mello et al., 2024; Hu et al., 2025a; Wang et al., 2025), with 50 replicates. $R^2$, RMSE and Lin's Concordance Correlation Coefficient (CCC) were used to evaluate model accuracy on the validation set. A higher $R^2$ and CCC, as well as a lower RMSE close to 0 means better model accuracy. These indices were calculated as follows:

$$ R^2 = 1 - \frac{\sqrt{\sum_{i}^{n} (\hat{y}_i - y_i)^2}}{\sqrt{\sum_{i}^{n} (y_i - \bar{y})^2}} \qquad (6) $$

$$ RMSE = \sqrt{\frac{\sum_{i}^{n} (\hat{y}_i - y_i)^2}{n}} \qquad (7) $$

$$ CCC = \frac{2\rho\sigma_y\sigma_y}{\sigma_y^2 + \sigma_y^2 + (\hat{y} - \bar{y})^2} \qquad (8) $$

where $y_i$ and $\hat{y}_i$ are observation and prediction for sample $i$, $n$ is the number of samples. $\rho$ is the Pearson correlation coefficient, $\sigma_y$ and $\sigma_y$ are the average and variance of $y$ (observed SOC) and $\hat{y}$ (predicted SOC).

## 2.10. Statistical calculation and analysis

Statistical analysis, correlation tests, covariate importance quantification, estimation of the Shapley values, and SOC prediction were conducted using R Studio Software (R Development Core Team, 2022). Experimental semivariograms were fitted using GS + 9.0 software (Gamma Design Software, Plainwell, MI, USA). The INLA-SPDE model was implemented using the “INLA” package in R language (http://www.rinla.org/). The RF was conducted by the “randomForest”
package in R 4.2.1 (Liaw and Wiener, 2002; R core Team, 2022) and the "raster" package (Hijmans, 2023; R core Team, 2022) to process and manipulate raster data. The GWR model and OK model was implemented using ArcGIS 10.8 (ESRI Inc., Version 10.8). Maps were generated using ArcGIS 10.8. Raster dataset processing was performed using SAGA GIS and ArcGIS 10.8. The workflow of this study is presented in Fig. 2.

# 3. Results

## 3.1. Descriptive statistics of SOC

The descriptive statistics of SOC are shown in Table 2. The SOC in all datasets ranges from 5.22 to 40.08 $g kg^{-1}$ with an average value of 17.91 $g kg^{-1}$. The calibration and validation data sets have an average

SOC content of 17.88 $g kg^{-1}$ and 17.98 $g kg^{-1}$, respectively. According to the Second National Soil Survey of China, the mean SOC content in our survey was higher than the mean SOC content in farmlands of China in 1980s (14.96 $g kg^{-1}$) (Shen, 1998). The standard deviation (SD) of all data sets ranged from 5.54 to 5.55 $g kg^{-1}$ while the skewness varied between 0.27 and 0.30, respectively. Additionally, the data structure of the three datasets was similar, indicating that the calibration set and validation set were representative and consistent.

## 3.2. The correlation between SOC and environmental covariates

Fig. 3 illustrates the correlations between SOC and various environmental covariates. Among them, the amount of straw return (SRA) exhibits the strongest and most significant positive correlations with SOC, playing a critical role in shaping its spatial distribution. Mantel's r

<figure><img src="image_4.png" /><figcaption>Fig. 2. The workflow of this study.</figcaption></figure>
<table><caption>**Table 2**
Descriptive statistics of soil SOC at sample locations.</caption><thead><tr><th></th><th>Database size</th><th>Minimum</th><th>Maximum</th><th>Mean</th><th>Median</th><th>SD</th><th>Skewness</th></tr></thead><tbody><tr><td>All dataset</td><td>16,050</td><td>5.22</td><td>40.08</td><td>17.91</td><td>17.57</td><td>5.54</td><td>0.29</td></tr><tr><td>Calibration</td><td>11,238</td><td>5.22</td><td>40.08</td><td>17.88</td><td>17.52</td><td>5.54</td><td>0.30</td></tr><tr><td>Validation</td><td>4812</td><td>6.73</td><td>37.76</td><td>17.98</td><td>17.75</td><td>5.55</td><td>0.27</td></tr></tbody></table>

Note: SD means standard deviation; All dataset means the dataset including all soil samples; Calibration means calibration dataset; Validation means validation dataset; N means the number of samples.

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Correlation analysis between SOC and environmental covariates.</figcaption></figure>

values for SOC with this factor exceeded 0.04, with *p* values below 0.001. In contrast, the relationships between SOC and other covariates, including Ele, VD, AS, SED, and soil erosion degree (SED), were not significant, with *r* values less than 0.02 and *p* values exceeding 0.05.

### 3.3. Spatial prediction of SOC by INLA-SPDE

#### 3.3.1. Environmental covariates selection

The CARS-PLS method was utilized to select the suitable environmental covariates for driving the spatial prediction models. As shown in Fig. 4, with an increase of the number of covariates used for predicting SOC, the RMSE was firstly sharply decreased and then kept stable, finally it showed steadily increasing trend. When the number of covariates was 19, the model had the best performance. Thus, the 19 covariates listed in Table 3 were used as covariates for mapping SOC.

#### 3.3.2. Summary of hyperparameters of INLA-SPDE model

Table 4 presents descriptive statistics of the fixed effect parameters of INLA-SPDE approach, including the mean value, standard deviation, and quantiles, which could indicate the influence of environmental covariates on SOC contents. When the absolute value of the parameter mean is the largest, it indicates the most significant impact of that environmental covariate on SOC. The positive or negative sign of the value reflects the direction of the effect, meaning SOC increases or decreases as the environmental covariate changes. In sight of the fixed

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Change of the RMSE of the CARS-PLS model with the increase of covariates number.</figcaption></figure>
<caption>**Table 3**</caption>

Selected optimal variables for spatial prediction modeling by CARS-PLS algorithm.

<table><caption>**Table 4**</caption><thead><tr><th>Environmental covariates</th><th>CARS-PLS(N = 19)</th></tr></thead><tbody><tr><td>Climate</td><td>MAT, MAP, MASR, LSTN</td></tr><tr><td>Terrain</td><td>Ele, VD</td></tr><tr><td>Soil properties</td><td>BD, ASi, PLD, Sand, EMg, AS</td></tr><tr><td>Biota</td><td>RS</td></tr><tr><td>Lithology</td><td>SC</td></tr><tr><td>Soil management</td><td>NF, SRA, IC, DC, SED</td></tr></tbody></table>

Summary statistics of the fixed effect parameters the INLA-SPDE model.

<table><thead><tr><th>Variable</th><th>Mean</th><th>2.5 % percentile</th><th>Median</th><th>97.5 % percentile</th><th>Mode</th><th>SD</th></tr></thead><tbody><tr><td>MAT</td><td>0.000</td><td>-0.001</td><td>0.000</td><td>0.001</td><td>0.000</td><td>0.001</td></tr><tr><td>MAP</td><td>-0.001</td><td>-0.003</td><td>-0.001</td><td>0.001</td><td>0.001</td><td>0.001</td></tr><tr><td>MASR</td><td>0.000</td><td>-0.001</td><td>0.000</td><td>0.001</td><td>0.000</td><td>0.001</td></tr><tr><td>LSTN</td><td>-0.001</td><td>-0.003</td><td>-0.001</td><td>0.001</td><td>0.001</td><td>0.001</td></tr><tr><td>Ele</td><td>0.000</td><td>-0.002</td><td>0.000</td><td>0.002</td><td>0.000</td><td>0.001</td></tr><tr><td>VD</td><td>-0.001</td><td>-0.005</td><td>-0.001</td><td>0.002</td><td>0.001</td><td>0.002</td></tr><tr><td>BD</td><td>-0.308</td><td>-0.581</td><td>-0.318</td><td>-0.198</td><td>-0.308</td><td>0.054</td></tr><tr><td>EMg</td><td>1.356</td><td>0.836</td><td>1.356</td><td>1.877</td><td>1.356</td><td>0.266</td></tr><tr><td>ASi</td><td>0.000</td><td>-0.008</td><td>0.000</td><td>0.008</td><td>0.000</td><td>0.004</td></tr><tr><td>AS</td><td>-0.001</td><td>-0.008</td><td>-0.001</td><td>0.007</td><td>0.001</td><td>0.004</td></tr><tr><td>PLD</td><td>0.022</td><td>-0.102</td><td>0.022</td><td>0.147</td><td>0.022</td><td>0.064</td></tr><tr><td>Sand</td><td>0.000</td><td>-0.002</td><td>0.000</td><td>0.003</td><td>0.000</td><td>0.001</td></tr><tr><td>NF</td><td>-0.018</td><td>-0.050</td><td>-0.018</td><td>0.013</td><td>0.018</td><td>0.016</td></tr><tr><td>SRA</td><td>3.238</td><td>3.068</td><td>3.207</td><td>3.437</td><td>3.237</td><td>0.067</td></tr><tr><td>IC</td><td>-0.139</td><td>-0.435</td><td>-0.139</td><td>0.157</td><td>0.139</td><td>0.151</td></tr><tr><td>DC</td><td>-0.106</td><td>-0.304</td><td>-0.106</td><td>0.092</td><td>0.106</td><td>0.101</td></tr><tr><td>SED</td><td>0.103</td><td>-0.056</td><td>0.103</td><td>0.263</td><td>0.103</td><td>0.081</td></tr><tr><td>SC</td><td>0.140</td><td>-0.041</td><td>0.140</td><td>0.321</td><td>0.140</td><td>0.092</td></tr><tr><td>RS</td><td>-0.142</td><td>-0.284</td><td>-0.142</td><td>-0.001</td><td>0.142</td><td>0.072</td></tr></tbody></table>

effect parameter, SRA had the most significant positive influence on the spatial pattern of SOC with a mean value of 3.238. There was a positive relationship between the SOC and EMg with a mean value of 1.356, representing the SOC contents increased with the SRA and EMg. In addition, the BD, Rotation system (RS), SC (soil class) and IC (irrigation capacity) also showed the significant impact on the spatial distribution of the SOC.

### 3.3.3. Comparison prediction accuracies of the INLA-SPDE with other models

To verify the feasibility and performance of INLA-SPDE, three other models including RF, GWR, and OK were used to predict the SOC. As presented in Fig. 5, INLA-SPDE performs much better than the GWR and OK with R² of 0.50. Compared with GWR and OK, the R² of INLA-SPDE increased 117.39 % and 38.89 %, respectively. INLA-SPDE also slightly outperformed RF with higher R² (0.50 Vs 0.41). Although the RMSE is slightly higher (3.92 Vs 2.73 g kg⁻¹). Furthermore, CCC provided an additional measure of agreement between observed and predicted SOC values, where INLA-SPDE (CCC = 0.65), slightly surpassing RF (CCC = 0.59). These results indicate that INLA-SPDE not only offers better predictive accuracy but also maintains stronger consistency with observed values compared to the other models.

### 3.3.4. The predicted SOC maps

Fig. 6 displays the maps of farmland SOC in Jiangxi Province produced by INLA-SPDE, RF, OK and GWR. Generally, all the maps showed a similar spatial pattern (Fig. 6). The high SOC values were mainly observed in the central part of the study area while the low SOC values were majorly distributed in the north part of the study area. However, there are some differences in the SOC spatial distributions predicted by each model. The predicted values produced by INLA-SPDE model had a broader range than RF and GWR while the predicted value range of OK is

closest to observed values.

## 3.4. Relative importance of covariates

The Fig. 7 ranks the contribution of different covariates to mapping SOC, evaluated by average absolute SHAP values derived from the RF model. It showed that SRA and MAP is the most important variable for mapping SOC, followed by MASR, MAT, and AS. Soil management (50.52 %) and climate variables (22.06 %) made largest contribution to modelling SOC, followed by soil properties (17.09 %), terrain (6.38 %), lithology (2.21 %) and biota factors (1.72 %) (Fig. 7). More specifically, SRA and MAP have the largest SHAP value range than others variables (Fig. 8). For instance, the SHAP values for SRA range from -8 to 4. Overall, soil samples with high SRA content showed positive effects on SOC while low SRA content showed negative effects on SOC. Similar patterns are observed in other covariates like MAP and MASR.

# 4. Discussion

## 4.1. Different methods for spatial prediction

The INLA-SPDE model explicitly describes uncertainty by treating model parameters as probabilistic variables with posterior probability density functions (Yang et al., 2025). This method provides predictions and uncertainty estimates by numerically approximating the marginal densities of hyperparameters and latent variables. Furthermore, it offers posterior distributions for hyperparameters, enhancing the interpretability of predictions (He and Wong, 2024; Yang et al., 2023). In contrast, machine learning methods like RF are data-driven approaches that average predictions from multiple decision trees to improve accuracy (Hu et al., 2024; Khaledian and Miller, 2020). However, RF's "black box" nature necessitates additional interpretability methods to understand its internal prediction mechanisms and variable contributions (Hu et al., 2025a; Wadoux et al., 2021). Despite its high predictive accuracy in this study, RF's interpretability lags behind that of INLA-SPDE.

Different methods also exhibit distinct characteristics in handling tail distributions (extreme values) of spatial data. In regions with uneven or sparse sampling, OK predictions tend to smooth the data, leading to averaged tail distribution predictions that may miss true extremes (McBratney et al., 2003; Hu et al., 2022). GWR may also struggle with tail distributions, as its accuracy in estimating local regression coefficients decreases in sparsely sampled areas, potentially leading to inaccurate extreme value predictions (Kupfer and Farris, 2006). While RF can handle nonlinear relationships and extreme values, its black-box nature makes it difficult to fully assess the credibility of its tail distribution predictions (Hu et al., 2025a; Wadoux et al., 2021). Conversely, INLA-SPDE's incorporation of random effects and spatial autocorrelation allows it to better capture spatial variability and extreme values. Its Bayesian-based uncertainty estimation further ensures more robust tail distribution predictions (Sun et al., 2021; Wu, 2021).

In summary, this study compared the performance of different methods in spatial prediction, exploring their inherent differences and the advantages and disadvantages for spatial data processing. Future research should further investigate hybrid models to combine multiple strengths, aiming to construct more accurate and interpretable spatial prediction models suitable for a broader range of soil mapping and land management applications.

## 4.2. Uncertainty analysis

We reported and compared the uncertainty of our prediction produced by INLA-SPDE and RF. The maps of 95 % upper and lower prediction limit produced by INLA-SPDE model and RF model were presented in Fig. 9. As shown in Fig. 9, INLA-SPDE model reported a clearly wider range of SOC prediction values. The uncertainty quantification of the RF model derives from the ensemble of all regression tree
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Scatter plot of SOC contents of external validation dataset ($N_v$ = 4812) and corresponding histogram (The red dash line means the 1:1 line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

results, based on frequentist methods, producing a smoother mapping result (Fig. 9) (Besic et al., 2025; Meinshausen, 2006). Numerous previous studies have shown that RF tends to overestimate extremely low values while underestimating high values (Han and Suh, 2024; Janitza and Hornung, 2018; Wolfensberger et al., 2021). In contrast, the INLA-SPDE method uses a joint posterior probability density function based on Bayesian theorem, making it more effective for uncertainty prediction (Huang et al., 2017; Poggio et al., 2016). The areas have high uncertainty in the SOC map produced by INLA-SPDE mainly appears at the edges of the study area and under sampled regions (Fig. 9), which is consistent with Yang et al. (2023). This may be due to edge effects from the model calculation method since the INLA-SPDE make predictions based on a mesh network (Huang et al., 2017, Hu, 2020). The uncertainty map can guide additional soil sampling to improve prediction accuracy in further effort. Studies have already used prediction uncertainty to guide supplementary soil sampling and improve mapping accuracy (Blackford et al., 2022; Bohn and Miller, 2024).

However, although our constructed 50 % of variance of the SOC, more than half variance is still not explained. The uncertainty of our predictions is attributed to several sources. Firstly, although our

constructed model considering lots of factors covering climate, terrain, biota, lithology, soil properties, and soil management factors, some other important factors that could affect the variation of SOC were still not included due to data availability and the cost, like microbial diversity and enzyme activity (Bhattacharyya et al., 2021; Han et al., 2024a; Han et al., 2024b; Xiao et al., 2025). Secondly, part of our covariates (e.g., MAP, MAT) has relative coarse original spatial resolution, this could also result in additional uncertainty in our predictions. Finally, limited by the visiting accessibility and cost, our soil sampling did not well spatially coverage the study area, which can lead to bias on the prediction results (Brus and DeGruijter, 1993; Brus, 2015). The map of prediction uncertainty we produced can guide in selecting the positions of supplemental sampling to reduce uncertainty and improve prediction accuracy (Brus and Riele, 2001; Li et al., 2016).

### 4.3. Spatially varying primary controls of SOC across study area

The SHAP method outputs decomposition factors evaluated by SHAP values, aiding in understanding the importance of each feature value for RF prediction outcomes (Baptista et al., 2022; Somala et al., 2024; Wang
<figure><img src="image_8.png" /><figcaption>Fig. 6. Spatial distribution of SOC predicted by model: INLA-SPDE(a); RF(b); GWR (c); OK (d).</figcaption></figure>

et al., 2025). The maps of SHAP values for the 19 covariates used for mapping SOC were produced and presented in Fig. S2. Based on that, we produced the map of the primary variables contributing to mapping SOC across the study area (Fig. 10). As showed in Fig. 10, SRA is the most important variable for mapping SOC in 69.25 % of the study area followed by MAP (17.20 %), MAT (5.18 %), MASR (3.75 %), and Sand

(2.17 %). Others variables had a proportion less than 1 %.

In pedology, SRA underscores the impact of agricultural practices on SOC, with straw return serving as a major source of organic inputs (Han et al., 2025; Xiong et al., 2025). This practice enriches soil organic matter, boosts microbial biomass, and enhances carbon turnover, thereby shaping the spatial distribution of SOC (Wang et al., 2021a).
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Relative importance of the predictor variables in the random forest model. The importance of each variable was determined by the absolute mean Shapley values.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Shapley values for each covariate for the RF model, i.e. the contribution of the covariate to the prediction of the SOC at that location. The colour scale represents the size of covariate value.</figcaption></figure>

Environmental factors like MAP, and MAT influence SOC dynamics by shaping climatic and physical conditions, though their impact is sec- ondary. The lower importance of MASR, Sand, and other variables suggests more localized or indirect effects. Integrating these insights enhances DSM's predictive accuracy and deepens understanding of SOC variability, supporting sustainable soil management.

Additionally, as illustrated in Fig. 11, farmland SOC in Jiangxi Province is primarily influenced by soil management (72.39 %) and climate (26.55 %), followed by soil properties (3.13 %), terrain (1.03 %) and lithology factors (0.01 %).

Combined with the results of SHAP, the effects of different envi- ronmental covariates on SOC changes can be better understood and validated. We used SHAP to analyze the impact of 19 covariates selected by the CARS-PLS method (Table 2). SHAP analysis indicated that the impact of covariates on SOC change followed this order: Soil manage- ment (69.26 %) > Climate (26.58 %) > Soil properties (3.13 %) > Terrain (1.03 %) > Lithology (<0.01 %) (Fig. 7), consistent with Hu et al. (2024). The lower C₀/(C₀ + C) ratio (Table S1) of SOC (0.5)

indicates that human activities, including soil management, signifi- cantly impact SOC changes, showing that SOC change is primarily affected by natural factors and human activities. SOC is not significantly affected by Terrain, Biota, and Lithology, likely because most farmland is in flat areas like the Poyang Lake Plain and river valley plain (Fig. 1), where topographic factors have less impact than soil management fac- tors (Li et al., 2020; Li et al., 2021; Zhao et al., 2022). This confirms the positive role of soil management measures in enhancing SOC content in Jiangxi Province's farmland.

Our study showed that SHAP method identified soil management as most closely related to SOC (Fig. 11). The SHAP model confirmed that the impact of SRA on SOC change was highest, aligning with findings by Hu et al. (2024) (Fig. 7, Fig. 8, Fig. 10). The SHAP values for SRA were generally higher, indicating that increased SRA positively impacted SOC (Fig. 7, Fig. 8). SRA significantly influences SOC within a specific range, stabilizing once it surpasses a threshold (Fig. 12). Straw contains nu- trients like carbon and nitrogen, which, after decomposition, help in- crease SOC content (Lal, 2004a; Li et al., 2024; Zhao et al., 2025). Additionally, SRA improves soil structure, increases porosity and water retention, enhances microbial activity, promotes organic matter decomposition, and increases SOC content (Blanco-Canqui and Lal, 2009; Li et al., 2024a; Li et al., 2024b). Liu et al. (2015) used a meta- analysis of 176 field studies to show that SRA improves SOC accumu- lation, soil quality, and crop yields.

Furthermore, the climate category factor in the environmental covariates plays a very important role in various processes that affect the SOC transformation mechanism. When the MAP reaches 1600 mm, its influence on SOC is minimal. However, once the MAP exceeds 1600 mm, the impact on SOC becomes increasingly pronounced (Fig. 12). Tem- perature and precipitation significantly impact the decomposition and mineralization rate of SOM (Galluzzi et al., 2024; Hu et al., 2022; Zhao et al., 2019). The increase in water replenishment can promote the accumulation of surface vegetation biomass and the improvement of net primary production efficiency, thereby strengthening the retention ef- fect of SOC in the soil. The decrease in temperature will inhibit the decomposition rate of organic matter in soil aggregates and significantly slow down the process of organic carbon loss. (Wang et al., 2021b; Zhuang et al., 2010). Furthermore, for MAP, the spatial distribution pattern of the SHAP values of the MAP in the study area further reveals the complex influence of precipitation on the spatial distribution of SOCD (Fig. S2). Specifically, except for a small part of the northwest and southwest regions where precipitation has a certain inhibitory effect on SOC accumulation, most areas in the study area show a positive zone, indicating that precipitation is one of the main positive factors pro- moting SOC accumulation in the study area. Additionally, MASR was identified as the third most important variable for predicting SOC ac- cording to the SHAP model in this study, while the fixed effect parameter of MASR in the INLA-SPDE model was 0 (Table 4). Specifically, the INLA-SPDE model indicated no significant direct relationship between MASR and SOC. Therefore, further research is needed to clarify MASR's impact on SOC changes.

Although the SHAP model shows the less significant relationship between SOC and covariate categories such as terrain and biology (Fig. 7, Fig. 8, Fig. 12), attention is still needed. The study area is mainly composed of plain landforms, which may result in a relatively low de- gree of spatial differentiation of each sampling point in terms of climatic conditions, soil types and topographic features (Zhang et al., 2022). Taking Ele as an example, the spatial heterogeneity of its SHAP value reveals the spatial differentiation characteristics of the relationship be- tween this factor and the SOC distribution in the study area (Fig. S2). Specifically, only the local areas in the central and southwestern parts of the study area showed positive SHAP values, indicating that the increase in Ele values has a promoting effect on the accumulation of SOC. However, the Ele-SHAP values in the vast majority of areas approach zero, suggesting that the regulatory effect of terrain elevation changes on SOC fixed memory is not significant in these wide areas. This
<figure><img src="image_11.png" /><figcaption>Fig. 9. The maps of 95 % upper prediction limit for INLA-SPDE (a) and RF (c) model; 95 % lower prediction limit for INLA-SPDE (b) and RF (d) model.</figcaption></figure>

phenomenon may be closely related to the limited scale of cultivated land and the relatively high flatness of the terrain in the study area (Chen et al., 2025). Meanwhile, the variability of some environmental factors at the regional scale is relatively weak, resulting in their driving effect on the organic carbon pool of farmland systems being in a secondary position (Wiesmeier et al., 2019). TN plays a pivotal role in SOC dynamics as it is a critical nutrient for microbial activity, directly influencing organic matter decomposition and carbon stabilization (Li et al., 2020). In addition, it is worth noting that higher TN levels promote microbial growth and enzyme activity, enhancing the formation of stable organic compounds through microbial byproducts (Shu et al., 2025). This makes TN not only an indicator of soil fertility but also a key predictor for SOC mapping (Hu et al., 2022; Hu et al., 2025b). However, since the covariate set after variable screening in this study does not contain TN, different covariate sets should be considered in the future to reflect the potential change process of SOC more accurately.

Our results provide important implications for enhancing digital soil

mapping in future work in our study area as well as other regions. More important, our constructed framework provide a better alternative tool on explaining DSM results with a more intuitive and easily understood way, offering an alternative to the traditional variable importance metrics produced by RF, which is commonly used in DSM (Ellili et al., 2019; Heung et al., 2014; Hu et al., 2024; Wiesmeier et al., 2010). More importantly, this method enables the evaluation of variable importance at a pixel level, allows for the generation of covariate contribution maps, and provides insights into the model’s spatial behavior—an uncommon application of SHAP values.

### 4.4. Limitations and perspectives

Firstly, although we have made important progress in this study, this study inevitably has some limitations like all studies. Firstly, although INLA-SPDE’s modeling time is significantly reduced compared to traditionally used Markov chain Monte Carlo (MCMC) (Bondo et al.,
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Map of the primary variables contributing to mapping SOC in the study area. The number followed each variable represents the proportion of the areas which corresponding variable was identified as primary variables contributing to mapping SOC.</figcaption></figure>

2024; Poggio et al., 2016; Wu, 2021), it still takes more time than RF due to fitting more hyperparameters than RF, with slight improvement in prediction accuracy. Beguin et al. (2017) reported that INLA-SPDE took about 134 times longer than the RF model, while Yang et al. (2023) found it took about 30 times longer. In this study, the time difference was approximately 47 times. If model complexity increases, such as denser mesh grids forming more small triangles, more covariates, or larger regions, time demands and computer power requirements will increase exponentially (Hu, 2020; Lindgren and Rue, 2015). In addition, this study was conducted at regional scale, and the applicability of INLA-SPDE for high-resolution mapping at the national scale or even continental scale still needs further exploration.

Additionally, the current study's evaluation framework exhibits methodological constraints related to geographical data partitioning. While the application of a 70%-30% split ratio represents a conventional validation approach across multiple disciplines, this strategy insufficiently addresses spatially correlated phenomena in the dataset. The presence of spatial interdependencies, where geographical proximity between observations correlates with measurement similarities, may substantially influence model outcomes. Conventional data partitioning protocols that neglect such spatial organization potentially compromise the predictive accuracy and external validity of analytical findings. Consequently, more sophisticated validation methodologies, particularly spatial cross-validation, should be implemented when analyzing spatially structured datasets. Subsequent studies addressing geographically dependent variables must adopt these enhanced evaluation protocols to properly account for spatial autocorrelation effects and improve the scientific rigor of model assessments. Future investigations should integrate these specialized validation frameworks to

<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Map of the most important groups contributing to farmland SOC prediction in Jiangxi Province. The size of each pie represents the proportion of the covariate group for which the specific value is labeled using a number.</figcaption></figure>

enhance evaluation accuracy in spatial correlation contexts, particularly when modeling complex geographical data patterns.

Secondly, some covariate maps (e.g., TP) in this study were generated using OK. While the substantial sampling size enhances interpolation accuracy, interpolation errors cannot be entirely eliminated (Hu et al., 2024). This introduces additional uncertainty into our SOC predictions. Notably, several covariates originally exceeded 30 m spatial resolution, necessitating resampling to 30 m using interpolation algorithm by ArcGIS. This resampling process inherently introduces secondary uncertainty into the modeling framework. Higher-resolution covariate inputs, ideally maintaining native spatial fidelity without resampling, would likely reduce propagation of these uncertainties and improve predictive performance.

Finally, our study was the exclusive use of CARS-PLS Recursive Feature Elimination for variable selection, without incorporating comparative assessments against other feature selection techniques or examining the potential contributions of original covariates. To enhance methodological rigor, subsequent research should implement structured comparisons of emerging feature selection strategies, with particular emphasis on cutting-edge algorithms such as Recursive Feature Elimination, Boruta feature ranking, particle swarm optimization, ant colony optimization, and interval-valued random frog methods (Ma et al., 2019). In addition, future studies should focus on incorporating more comprehensive and higher-resolution covariates into the modeling process to improve prediction accuracy and reduce uncertainty in subsequent studies.

To fully realize the application potential of INLA-SPDE in DSM, further experimental verification is needed. More attention can be paid to investigating mesh grid settings to balance mapping accuracy and time spent, and quantitatively evaluate uncertainty results.
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Dependence plot of covariates for the SOC. The x axis shows the covariate value and the y axis shows the SHAP value of that covariate. The red line is a smooth curve fitted to the SHAP value for visualization purposes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

**5. Conclusion**

This study proved the feasibility of INLA-SPDE on mapping soil properties in hilly regions in farmlands of Jiangxi Province, south China. In comparison with existing studies which use geostatistical and machine learning methods like RF, this study makes progress in the following aspects:

(1) INLA-SPDE model shows clearly better performance than classical geostatistical methods (OK and GWR) and slightly outperformed RF.

(2) INLA-SPDE has an advantage on reporting the posterior distribution of the model hyperparameters.

(3) This study demonstrated that INLA-SPDE can deal with large data size (16,050 samples), and could include large and various kinds of covariates (19 covariates) in the modelling process.

(4) By integrating INLA-SPDE and interpretable machine learning model, we can quantify the overall importance of different variables for mapping soil properties and the spatially varying primary variables for mapping soil properties across the study area at pixel level.
<figure><img src="image_15.png" /><figcaption>Fig. 12. (continued).</figcaption></figure>

Based on these findings, we recommend the application of INLA-SPDE for mapping soil properties as an alternative to traditional geostatistical models and RF. In practice, to effectively apply this framework, it is important to reasonably allocate computational resources to handle large datasets and perform spatial modeling efficiently. Additionally, employing well-designed sampling strategies is crucial for capturing spatial variability and ensuring robust predictions. It is also essential to carefully select relevant covariates that accurately represent the environmental factors influencing soil properties. The model's flexibility allows it to be customized for different environmental and agricultural contexts, making it suitable for a wide range of soil mapping applications.

Therefore, we recommend the application of INLA-SPDE to mapping soil properties as an alternative method for geostatistical model and RF. In addition, more powerful computation resources, better sampling designing scheme and more suitable covariates is also required to improve the prediction accuracy of INLA-SPDE model. Moreover, further efforts are also necessary to confirm the feasibility and applicability of INLA-SPDE on mapping different soil properties in large spatial scale, e.g., national or even continental scale. Our study made positive contribution to deal with the *Challenge 3: In what ways can we use data-driven models to learn about pedological processes?* of the ten challenges for the future of pedometrics reported by Wadoux et al. (2021). By applying the INLA-SPDE model, our study demonstrates how data-driven approaches can uncover the spatial patterns and relationships governing SOC dynamics. This model not only accounts for spatial heterogeneity and dependencies but also highlights key pedological processes, such as the influence of nitrogen, straw return, and climatic factors on SOC distribution. Furthermore, integrating SHAP values provides a deeper and more interpretable understanding of the interactions between environmental covariates and soil processes. These advancements showcase the potential of data-driven models to bridge the gap between predictive mapping and the exploration of pedological mechanisms, offering valuable insights for soil science and management.

### CRediT authorship contribution statement

**Bifeng Hu:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Funding acquisition, Conceptualization. **Yibo Geng:** Writing – review & editing, Data curation. **Hanjie Ni:** Writing – original draft, Visualization, Validation, Software, Data curation. **Zhou Shi:** Supervision, Resources, Conceptualization. **Zheng Wang:** Writing – review & editing, Software. **Nan Wang:** Writing – review & editing, Data curation. **Jipeng Luo:** Writing – review & editing. **Modian Xie:** Writing – review & editing, Writing – original draft, Data curation. **Qian Zou:** Software, Data curation. **Thomas Optiz:** Software. **Hongyi Li:** Supervision, Resources, Project administration, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

This study was supported by the National Science Foundation of China (No. 42201073), the Jiangxi “Double Thousand plan” (No. jxsq202301091), the Natural Science Foundation of Jiangxi Province (No. 20232BAB213058), the National Science Foundation of China (No. 42471060), and Open funding of Key Laboratory of Data Science in Finance and Economics, Jiangxi University of Finance and Economics.

## Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.geoderma.2025.117446.

### Data availability

The authors do not have permission to share data.
## References

Bai, Z., Chen, S., Hong, Y., Hu, B., Luo, D., Peng, J., Shi, Z., 2023. Estimation of soil inorganic carbon with visible near-infrared spectroscopy coupling of variable selection and deep learning in arid region of China. Geoderma 437, 116589. https://doi.org/10.1016/j.geoderma.2023.116589.

Baptista, M.L., Goebel, K., Henriques, E.M., 2022. Relation between prognostics predictor evaluation metrics and local interpretability SHAP values. Artif. Intell. 306, 103667. https://doi.org/10.1016/j.artint.2022.103667.

Barkhordari, M.S., Qi, C., 2025. Prediction of soil organic carbon concentration in European soils: a dimensionality reduction and ensemble learning approach. J. Hazard. Mater. Adv. 100604. https://doi.org/10.1016/j.hazadv.2025.100604.

Batjes, N., 1996. Total carbon and nitrogen in the soils of the world. Eur. J. Soil Sci. 47 (2), 151–163. https://doi.org/10.1111/j.1365-2389.1996.tb01386.x.

Beguin, J., Fuglstad, G., Mansuy, N., Paré, D., 2017. Predicting soil properties in the Canadian boreal forest with limited data: Comparison of spatial and non-spatial statistical approaches. Geoderma 306, 195–205. https://doi.org/10.1016/j.geoderma.2017.06.016.

Besic, N., Picard, N., Vega, C., Bontemps, J., Hertzog, L., Renaud, J., Fogel, F., Schwartz, M., Pelissier-Tanon, A., Destouet, G., Mortier, F., Planells-Rodriguez, M., Ciais, P., 2025. Remote-sensing-based forest canopy height mapping: some models are useful, but might they provide us with even more insights when combined? Geosci. Model Dev. 18 (2), 337–359. https://doi.org/10.5194/gmd-18-337-2025.

Bhattacharyya, R., Rabbi, S.M., Zhang, Y., Young, I.M., Jones, A.R., Dennis, P.G., Menzies, N.W., Kopitte, P.M., Dalal, R.C., 2021. Soil organic carbon is significantly associated with the pore geometry, microbial diversity and enzyme activity of the macro-aggregates under different land uses. Sci. Total Environ. 778, 146286. https://doi.org/10.1016/j.scitotenv.2021.146286.

Blackford, C., Heung, B., Webster, K.L., 2022. Incorporating spatial uncertainty maps into soil sampling improves digital soil mapping classification accuracy in Ontario, Canada. Geoderma Regional 29, e00495. https://doi.org/10.1016/j.geodrs.2022.e00495.

Blanco-Canqui, H., Lal, R., 2009. Crop residue removal impacts on soil productivity and environmental quality. Crit. Rev. Plant Sci. 28 (3), 139–163. https://doi.org/10.1080/07352680902776507.

Blangiardo, M., Cameletti, M., Baio, G., Rue, H., 2013. Spatial and spatio-temporal models with R-INLA. Spatial Spatio-temporal Epidemiol. 7, 39–55. https://doi.org/10.1016/j.sste.2013.07.003.

Blangiardo, M., & Cameletti, M. (2015). Spatial and Spatio-temporal Bayesian Models with R-INLA. doi: 10.1002/9781118950203.

Bohn, M.P., Miller, B.A., 2024. Locally enhanced digital soil mapping in support of a bottom-up approach is more accurate than conventional soil mapping and top-down digital soil mapping. Geoderma 442, 116781. https://doi.org/10.1016/j.geoderma.2024.116781.

Bondo, K.J., Rosenberry, C.S., Stainbrook, D., Walter, W.D., 2024. Comparing risk of chronic wasting disease occurrence using Bayesian hierarchical spatial models and different surveillance types. Ecol. Model. 493, 110756. https://doi.org/10.1016/j.ecolmodel.2024.110756.

Breiman, L., Friedman, J. H., Olshen, R. A., & Stone, C. J. (2017). Classification and regression trees. In Routledge Books. doi: 10.1201/9781315139470.

Brunsdon, C., Fotheringham, A.S., Charlton, M.E., 1996. Geographically weighted regression: a method for exploring spatial nonstationarity. Geogr. Anal. 28 (4), 281–298. https://doi.org/10.1111/j.1538-4632.1996.tb00936.x.

Brus, D.J., DeGruijter, J.J., 1993. Design-based versus model-based estimates of spatial means: Theory and application in environmental soil science. Environmetrics 4 (2), 123–152. https://doi.org/10.1002/env.3170040202.

Brus, D., Riele, W.T., 2001. Design-based regression estimators for spatial means of soil properties: the use of two-phase sampling when the means of the auxiliary variables are unknown. Geoderma 104 (3–4), 257–279. https://doi.org/10.1016/s0016-7061 (01)00084-2.

Brus, D., 2015. Balanced sampling: a versatile sampling approach for statistical soil surveys. Geoderma 253–254, 111–121. https://doi.org/10.1016/j.geoderma.2015.04.009.

Camera, C., Zomeni, Z., Noller, J.S., Zissimos, A.M., Christoforou, I.C., Bruggeman, A., 2017. A high resolution map of soil types and physical properties for Cyprus: a digital soil mapping optimization. Geoderma 285, 35–49. https://doi.org/10.1016/j.geoderma.2016.09.019.

Chen, S., Arrouays, D., Mulder, V.L., Poggio, L., Minasny, B., Roudier, P., Libohova, Z., Lagacherie, P., Shi, Z., Hannam, J., Meersmans, J., Richer-De-Forges, A.C., Walter, C., 2022. Digital mapping of GlobalSoilMap soil properties at a broad scale: a review. Geoderma 409, 115567. https://doi.org/10.1016/j.geoderma.2021.115567.

Chen, S., Xu, H., Xu, D., Ji, W., Li, S., Yang, M., Shi, Z., 2021. Evaluating validation strategies on the performance of soil property prediction from regional to continental spectral data. Geoderma 400, 115159. https://doi.org/10.1016/j.geoderma.2021.115159.

Cressie, N., 1988. Spatial prediction and ordinary kriging. Math. Geol. 20 (4), 405–421. https://doi.org/10.1007/bf00892986.

Dai, L., Xue, J., Lu, R., Wang, Z., Chen, Z., Yu, Q., Shi, Z., Chen, S., 2024. In-situ prediction of soil organic carbon contents in wheat-rice rotation fields via visible near-infrared spectroscopy. Soil Environ. Health. 2 (4), 100113. https://doi.org/10.1016/j.seh.2024.100113.

de Mello, D.C., Silvero, N.E., Miller, B.A., Rosin, N.A., Rosas, J.T.F., dos Anjos Bartsch, B., Dematté, J.A.M., 2024. Mapping soil drainage classes: Comparing expert knowledge and machine learning strategies. Soil Adv. 100028. https://doi.org/10.1016/j.soilad.2024.100028.

Dong, S., Zhang, J., Li, Y., Liu, S., Dong, Q., Zhou, H., Yeomans, J., Li, Y., Li, S., Gao, X., 2019. Effect of grassland degradation on aggregate-associated soil organic carbon of alpine grassland ecosystems in the Qinghai-Tibetan Plateau. Eur. J. Soil Sci. 71 (1), 69–79. https://doi.org/10.1111/ejss.12835.

Ellili, Y., Walter, C., Michot, D., Pichelin, P., Lemercier, B., 2019. Mapping soil organic carbon stock change by soil monitoring and digital soil mapping at the landscape scale. Geoderma 351, 1–8. https://doi.org/10.1016/j.geoderma.2019.03.005.

Emadi, M., Baghernejad, M., 2014. Comparison of spatial interpolation techniques for mapping soil pH and salinity in agricultural coastal areas, northern Iran. Arch. Agronomy Soil Sci./archiv Für Acker- Und Pflanzenbau Und Bodenkunde 60 (9), 1315–1327. https://doi.org/10.1080/03650340.2014.880837.

Fotheringham, A. S., Brunsdon, C., & Charlton, M. (2002). Geographically Weighted Regression: The analysis of spatially varying relationships. doi: 10.1111/j.1538-4632.2003.tb01114.x.

Gacutan, J., Tait, H., Clark, G.F., 2024. Assessing drivers of estuarine debris using a Bayesian spatial modelling approach (INLA-SPDE). Estuar. Coast. Shelf Sci. 296, 108592. https://doi.org/10.1016/j.ecss.2023.108592.

Galluzzi, G., Plaza, C., Priori, S., Giannetta, B., Zaccone, C., 2024. Soil organic matter dynamics and stability: climate vs. In: Time. Science of the Total Environment. https://doi.org/10.1016/j.scitotenv.2024.172441.

Gao, H., Wang, S., Shen, D., Li, W., Zhang, Y., Deng, J., Jiao, J., Huang, P., Ma, Q., Huang, Y., Yuan, F., Wang, X., Zhang, J., Meersman, J., Li, Y., Li, Q., Cheng, Y., 2024. Protected cultivation can mitigate carbon emissions for tomato production. Agric Commun. 2 (4), 100065. https://doi.org/10.1016/j.agrcom.2024.100065.

Georgiou, K., Jackson, R.B., Vindusková, O., Abramoff, R.Z., Ahlström, A., Feng, W., Harden, J.W., Pellegrini, A.F.A., Polley, H.W., Soong, J.L., Riley, W.J., Torn, M.S., 2022. Global stocks and capacity of mineral-associated soil organic carbon. Nat. Commun. 13 (1). https://doi.org/10.1038/s41467-022-31540-9.

Giannini-Kurina, F., Borello, J., Canías, I., Hang, S., Balzarini, M., 2022. Mapping atrazine persistence in soils of central Argentina using INLA. Soil Tillage Res. 219, 105320. https://doi.org/10.1016/j.still.2022.105320.

Gomes, L.C., Faria, R.M., De Souza, E., Veloso, G.V., Schaefer, C.E.G., Filho, E.I.F., 2019. Modelling and mapping soil organic carbon stocks in Brazil. Geoderma 340, 337–350. https://doi.org/10.1016/j.geoderma.2019.01.007.

Han, H., Suh, J., 2024. Spatial prediction of soil contaminants using a hybrid random Forest-Ordinary kriging model. Appl. Sci. 14 (4), 1666. https://doi.org/10.3390/app14041666.

Han, X., Doménech-Pascual, A., Casas-Ruiz, J.P., Donhauser, J., Jordaan, K., Ramond, J., Priemá, A., Romaní, A.M., Frossard, A., 2024a. Soil organic matter properties drive microbial enzyme activities and greenhouse gas fluxes along an elevational gradient. Geoderma 449, 116993. https://doi.org/10.1016/j.geoderma.2024.116993.

Han, Y., Lu, L., Wang, L., Liu, Z., Huang, P., Chen, S., Li, Y., Sun, M., He, C., Wang, J., Yan, Y., 2024b. In-situ straw return, combined with plastic film use, influences soil properties and tomato quality and yield in greenhouse conditions. Agric. Commun. 2 (1), 100028. https://doi.org/10.1016/j.agrcom.2024.100028.

He, S., Wong, S.W., 2024. Spatio-temporal data fusion for the analysis of in situ and remote sensing data using the INLA-SPDE approach. Spatial Stat. 64, 100863. https://doi.org/10.1016/j.spasta.2024.100863.

Hengl, T., Heuvelink, G.B.M., Kempen, B., Leenaars, J.G.B., Walsh, M.G., Shepherd, K.D., Sila, A., MacMillan, R.A., De Jesus, J.M., Tamene, L., Tondoh, J.E., 2015. Mapping Soil Properties of Africa at 250 m Resolution: Random Forests significantly improve current predictions. PLoS One 10 (6), e0125814. https://doi.org/10.1371/journal.pone.0125814.

Heung, B., Bulmer, C.E., Schmidt, M.G., 2014. Predictive soil parent material mapping at a regional-scale: a Random Forest approach. Geoderma 214–215, 141–154. https://doi.org/10.1016/j.geoderma.2013.09.016.

Hijmans R (2023). raster: Geographic Data Analysis and Modeling. R package version 3.6-26. doi: 10.32614/CRAN/package.raster.

Hu, B., Xue, J., Zhou, Y., Shao, S., Fu, Z., Li, Y., Chen, S., Qi, L., Shi, Z., 2020a. Modelling bioaccumulation of heavy metals in soil-crop ecosystems and identifying its controlling factors using machine learning. Environ. Pollut. 262, 114308. https://doi.org/10.1016/j.envpol.2020.114308.

Hu, B. (2020). Joint space-time modelling of soil properties using data from national soil information system. A Bayesian modelling approach on phosphorus across mainland France. Orleans University.

Hu, B., Saby, N., Bourennane, H., Opitz, T., Denoroy, P., Lemercier, B., & Shi, Z. (2020, May). Bayesian spatio-temporal modeling of soil phosphorus in Britany in western France (1995-2014) with INLA-SPDE. In EGU General Assembly Conference Abstracts (p. 19135). doi: 10.5194/egusphere-egu2020-19135.

Hu, B., Xie, M., Li, H., Zhao, W., Hu, J., Jiang, Y., Ji, W., Li, S., Hong, Y., Yang, M., Optiz, T., 2022. Stoichiometry of soil carbon, nitrogen, and phosphorus in farmland soils in southern China: Spatial pattern and related dominates. Catena 217, 106468. https://doi.org/10.1016/j.catena.2022.106468.

Hu, B., Ni, H., Xie, M., Li, H., Wen, Y., Chen, S., Zhou, Y., Teng, H., Bourennane, H., Shi, Z., 2023a. Mapping soil organic matter and identifying potential controls in the farmland of Southern China: Integration of multi-source data, machine learning and geostatistics. Land Degrad. Dev. 34 (17), 5468–5485. https://doi.org/10.1002/ldr.4858.

Hu, B., Xie, M., He, R., Shi, Z., Zhou, Y., Ni, H., Li, H., 2023b. Spatio-temporal evolution of cropland ecosystem services value and its spatially varying dominate over the past two decades: a case study in Jiangxi Province in southern China. J. Clean. Prod. 427, 139228. https://doi.org/10.1016/j.jclepro.2023.139228.

Hu, B., Xie, M., Zhou, Y., Chen, S., Zhou, Y., Ni, H., Peng, J., Ji, W., Hong, Y., Li, H., Shi, Z., 2024. A high-resolution map of soil organic carbon in cropland of Southern China. Catena 237, 107813. https://doi.org/10.1016/j.catena.2024.107813.
Hu, B., Geng, Y., Lin, Y., Ni, H., Xie, M., Wang, N., Hu, J., Zou, Q., Chen, S., Zhou, Y., Li, H., Shi, Z., 2025a. Comparison of machine learning and geostatistical methods on mapping soil organic carbon density in regional croplands and visualizing its location-specific dominators via interpretable model. Land Degrad. Dev. 36 (10), 3405-3427. https://doi.org/10.1002/ldr.5573.

Hu, B., Geng, Y., Shi, K., Xie, M., Ni, H., Zhu, Q., Qiu, Y., Zhang, Y., Bourennane, H., 2025b. Fine-resolution baseline maps of soil nutrients in farmland of Jiangxi Province using digital soil mapping and interpretable machine learning. Catena 249, 108635. https://doi.org/10.1016/j.catena.2024.108635.

Huang, J., Malone, B.P., Minasny, B., McBratney, A.B., Triantafilis, J., 2017. Evaluating a Bayesian modelling approach (INLA-SPDE) for environmental mapping. Sci. Total Environ. 609, 621-632. https://doi.org/10.1016/j.scitotenv.2017.07.201.

Janitza, S., Hornung, R., 2018. On the overestimation of random forest's out-of-bag error. PLoS One 13 (8), e0201904. https://doi.org/10.1371/journal.pone.0201904.

Jiang, Y., Li, F., Gong, Y., Yang, X., Zhang, Z., 2025. Multiple environmental variables as covariates to improve the accuracy of spatial prediction models for SOM on Karst Aera. Land Degrad. Dev. https://doi.org/10.1002/ldr.5454.

Keser, S., Duzgun, S., Aksoy, A., 2012. Application of spatial and non-spatial data analysis in determination of the factors that impact municipal solid waste generation rates in Turkey. Waste Manag. 32 (3), 359-371. https://doi.org/10.1016/j.wasman.2011.10.017.

Keskin, H., Grunwald, S., Harris, W.G., 2019. Digital mapping of soil carbon fractions with machine learning. Geoderma 339, 40-58. https://doi.org/10.1016/j.geoderma.2018.12.037.

Khaledian, Y., Miller, B.A., 2020. Selecting appropriate machine learning methods for digital soil mapping. App. Math. Model. 81, 401-418. https://doi.org/10.1016/j.aplmmd.2019.12.016.

Kuang, Z., Zhao, Y., Yang, X., 2023. Machine learning approaches for plant miRNA prediction: challenges, advancements, and future directions. Agric. Commun. 1 (2), 100014. https://doi.org/10.1016/j.agrcom.2023.100014.

Kumar, S., Lal, R., Liu, D., 2012. A geographically weighted regression kriging approach for mapping soil organic carbon stock. Geoderma 189-190, 627-634. https://doi.org/10.1016/j.geoderma.2012.05.022.

Kupfer, J.A., Farris, C.A., 2006. Incorporating spatial non-stationarity of regression coefficients into predictive vegetation models. Landsc. Ecol. 22 (6), 837-852. https://doi.org/10.1007/s10980-006-9058-2.

Lal, R., 2004a. Soil carbon sequestration impacts on global climate change and food security. Science 304 (5677), 1623-1627. https://doi.org/10.1126/science.1097396.

Lal, R., 2004b. Soil carbon sequestration to mitigate climate change. Geoderma 123 (1-2), 1-22. https://doi.org/10.1016/j.geoderma.2004.01.032.

Li, H., Xie, M., Hu, B., Zhou, Y., Chen, S., He, R., Li, Y., Teng, H., Zhao, W., Jia, X., Peng, J., Bourennane, H., & Shi, Z. (2024). Effects of straw return on soil carbon sequestration, soil nutrients and rice yield of in acidic farmland soil of Southern China. Environment, Development and Sustainability, pp.1-21. doi: 10.1007/s10668-024-04895-w.

Li, H., Yang, J., Hu, B., Geng, Y., Zhu, Q., Hong, Y., Lin, Y., Peng, J., Ji, W., Chen, S., Shi, Z., 2025. Improving in-situ spectral estimation of wetland soil organic carbon by integrating multiple optimization strategies. Catena 255, 109078. https://doi.org/10.1016/j.catena.2025.109078.

Li, H., Liang, Y., Xu, Q., Cao, D., 2009. Key wavelengths screening using adaptive competitive reweighted sampling method for multivariate calibration. Anal. Chim. Acta 648 (1), 77-84. https://doi.org/10.1016/j.aca.2009.06.046.

Li, J., Jian, S., Lane, C.S., Lu, Y., Ye, He, X., Wang, G., Mayes, M.A., Dzantor, K.E., Hui, D., 2020. Effects of nitrogen fertilization and bioenergy crop type on topsoil organic carbon and total nitrogen contents in middle Tennessee USA. PLoS One 15 (3), e0230688. https://doi.org/10.1371/journal.pone.0230688.

Li, J., Zhao, J., Liao, X., Hu, P., Wang, W., Ling, Q., Xie, L., Xiao, J., Zhang, W., Wang, K., 2024b. Pathways of soil organic carbon accumulation are related to microbial life history strategies in fertilized agroecosystems. Sci. Total Environ. 927, 172191. https://doi.org/10.1016/j.scitotenv.2024.172191.

Li, N., Zare, E., Huang, J., Triantafilis, J., 2018. Mapping soil cation-exchange capacity using bayesian modeling and proximal sensors at the field scale. Soil Sci. Soc. Am. J. 82 (5), 1203-1216. https://doi.org/10.2136/sssaj2017.10.0356.

Li, R., Gao, Y., Chen, Q., Li, Z., Gao, F., Meng, Q., Li, T., Liu, A., Wang, Q., Wu, L., Wang, Y., Liu, Z., Zhang, M., 2021. Blended controlled-release nitrogen fertilizer with straw returning improved soil nitrogen availability, soil microbial community, and root morphology of wheat. Soil Tillage Res. 212, 105045. https://doi.org/10.1016/j.still.2021.105045.

Li, Y., Zhu, A., Shi, Z., Liu, J., Du, F., 2016. Supplemental sampling for digital soil mapping based on prediction uncertainty from both the feature domain and the spatial domain. Geoderma 284, 73-84. https://doi.org/10.1016/j.geoderma.2016.08.013.

Liaw, A., Wiener, M., 2002. Classification and Regression by randomForest. R News 2 (3), 18-22.

Lindgren, F., Rue, H., Lindström, J., 2011. An explicit link between gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. J. R. Stat. Soc. Ser. B Stat Methodol. 73 (4), 423-498. https://doi.org/10.1111/j.1467-9868.2011.00777.x.

Lindgren, F., Rue, H., 2015. Bayesian spatial modelling withR-INLA. J. Stat. Softw. 63 (19). https://doi.org/10.18637/jss.v063i19.

Liu, Y., Guo, L., Jiang, Q., Zhang, H., Chen, Y., 2015. Comparing geostatistical techniques to predict SOC stocks. Soil Tillage Res. 148, 46-58. https://doi.org/10.1016/j.still.2014.12.002.

Lombardo, L., Optiz, T., Huser, R., 2018. Point process-based modeling of multiple debris flow landslides using INLA: an application to the 2009 Messina disaster. Stoch. Env. Res. Risk A. 32, 2179-2198. https://doi.org/10.1007/s00477-018-1518-0.

Ma, Y., Minasny, B., Malone, B.P., McBratney, A.B., 2019. Pedology and digital soil mapping (DSM). Eur. J. Soil. Sci. 70 (2), 216-235. https://doi.org/10.1111/ejs.12790.

Ma, Y., Minasny, B., Viaud, V., Walter, C., Malone, B., McBratney, A., 2023. Modelling the whole profile soil organic carbon dynamics considering soil redistribution under future climate change and landscape projections over the Lower Hunter Valley, Australia. Land-Base 12 (1), 255. https://doi.org/10.3390/land12010255.

Ma, Y., Minasny, B., Wu, C., 2017. Mapping key soil properties to support agricultural production in Eastern China. Geoderma Reg. 10, 144-153. https://doi.org/10.1016/j.geodrs.2017.06.002.

Mahjenabadi, V., Mousavi, S.R., Rahmani, A., Karami, A., Rahmani, H.A., Khavazi, K., Rezaei, M., 2022. Digital mapping of soil biological and wheat yield using remotely sensed, soil chemical data and machine learning approaches. Comput. Electron. Agric. 197, 106978. https://doi.org/10.1016/j.compag.2022.106978.

Martins, I.F., Teixeira, A.L., Pinheiro, L., Falcao, A.O., 2012. A Bayesian approach to in silico Blood-Brain barrier penetration modeling. J. Chem. Inf. Model. 52 (6), 1686-1697. https://doi.org/10.1021/ci300124c.

Matheron, G. (1971). The Theory of Regionalized Variables and its Applications.

Matinfar, H.R., Maghsodi, Z., Mousavi, S.R., Rahmani, A., 2021. Evaluation and prediction of topsoil organic carbon using machine learning and hybrid models at a field-scale. Catena 202, 105258. https://doi.org/10.1016/j.catena.2021.105258.

McBratney, A., Santos, M.M., Minasny, B., 2003. On digital soil mapping. Geoderma 117 (1-2), 3-52. https://doi.org/10.1016/s0016-7061(03)00223-4.

Minasny, B., McBratney, A.B., 2005. The Matèrn function as a general model for soil variograms. Geoderma 128 (3-4), 192-207. https://doi.org/10.1016/j.geoderma.2005.04.003.

Minasny, B., Vrugt, J.A., McBratney, A.B., 2011. Confronting uncertainty in model-based geostatistics using Markov Chain Monte Carlo simulation. Geoderma 163 (3-4), 150-162. https://doi.org/10.1016/j.geoderma.2011.03.011.

Minasny, B., McBratney, A., 2016. Digital soil mapping: a brief history and some lessons. Geoderma 264, 301-311. https://doi.org/10.1016/j.geoderma.2015.07.017.

Minasny, B., McBratney, A.B., Wadoux, A.M., Akoeb, E.N., Sabrina, T., 2020. Precocious 19th century soil carbon science. Geoderma Reg. 22, e00306. https://doi.org/10.1016/j.geodrs.2020.e00306.

Molla, A., Zhang, W., Zuo, S., Ren, Y., Han, J., 2023. A machine learning and geostatistical hybrid method to improve spatial prediction accuracy of soil potentially toxic elements. Stoch. Env. Res. Risk A. 37 (2), 681-696. https://doi.org/10.1007/s00477-022-02284-1.

Moraga, P., Dean, C., Inoue, J., Morawiecki, P., Noureen, S.R., Wang, F., 2021. Bayesian spatial modelling of geostatistical data using INLA and SPDE methods: a case study predicting malaria risk in Mozambique. Spatial Spatio-temporal Epidemiol. 39, 100440. https://doi.org/10.1016/j.sste.2021.100440.

Nash, J.F., 1950. Equilibrium points in n-person games. Proc. Natl. Acad. Sci. 36 (1), 48-49. https://doi.org/10.1073/pnas.36.1.48.

Neeraj, M.S., Kn, P., Deori, C., Kumar, K., Ishaq, T., Gouthami, D., Gourkhede, P.H., Kashyap, M., Das, S., 2025. The impact of conservation agriculture on soil carbon sequestration and climate mitigation: a comprehensive review. Int. J. Res. Agronomy 8 (1), 151-158. https://doi.org/10.3354/2618060x.2025.v8.i1c.2362.

Odeh, J., McBratney, A., Chittleborough, D., 1995. Further results on prediction of soil properties from terrain attributes: heterotopic cokriging and regression-kriging. Geoderma 67 (3-4), 215-226. https://doi.org/10.1016/0016-7061(95)00007-b.

Padarian, J., McBratney, A.B., Minasny, B., 2020. Game theory interpretation of digital soil mapping convolutional neural networks. Soil 6 (2), 389-397. https://doi.org/0.5194/soil-6-389-2020.

Paige, J., Fuglstad, G.A., Riebler, A., Wakefield, J., 2022. Bayesian multiresolution modeling of georeferenced data: an extension of 'LatticeKrig'. Comput. Stat. Data Anal. 173, 107503. https://doi.org/10.1016/j.csda.2022.107503.

Petersen, C., Trautner, A., Hansen, S., 2008. Spatio-temporal variation of anisotropy of saturated hydraulic conductivity in a tilled sandy loam soil. Soil Tillage Res. 100 (1-2), 108-113. https://doi.org/10.1016/j.still.2008.05.004.

Poggio, L., Gimona, A., Spezia, L., Brewer, M.J., 2016. Bayesian spatial modelling of soil properties and their uncertainty: the example of soil organic matter in Scotland using R-INLA. Geoderma 277, 69-82. https://doi.org/10.1016/j.geoderma.2016.04.026.

Post, W.M., Kwon, K.C., 2000. Soil carbon sequestration and land-use change: processes and potential. Glob. Chang. Biol. 6 (3), 317-327. https://doi.org/10.1046/j.1365-2486.2000.00308.x.

R Core Team, 2022. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Lanzhou, China. (URL) https://cran.rproject.org/.

Rasmusen, E., 1989. Games and Information: An Introduction to Game Theory. Wiley-Blackwell.

Rue, H., Martino, S., Chopin, N., 2009. Approximate Bayesian Inference for Latent Gaussian models by using Integrated Nested Laplace Approximations. J. R. Stat. Soc. Ser. B Stat Methodol. 71 (2), 319-392. https://doi.org/10.1111/j.1467-9868.2008.00700.x.

Ryan, E.G., Drovandi, C., McGree, J., Pettitt, A., 2015. A review of modern computational algorithms for bayesian optimal design. Int. Stat. Rev. 84 (1), 128-154. https://doi.org/10.1111/insr.12107.

Shapley, L., 1997. 7. A Value for N-Person Games. Contributions to the Theory of Games II (1953) 307-317. In Princeton University Press eBooks (pp. 69-79). doi: 10.1515/9781400829156-012.

Shen, S.M., 1998. Soil fertility in China (In Chinese). China Agriculture Press, Beijing.
Shen, Q., Wang, Y., Wang, X., Liu, X., Zhang, X., Zhang, S., 2019. Comparing interpolation methods to predict soil total phosphorus in the Mollisol area of Northeast China. Catena 174, 59–72. https://doi.org/10.1016/j.catena.2018.10.052.

Shu, Q., Gao, S., Liu, X., Yao, Z., Wu, H., Qi, L., Zhang, X., 2025. Soil enzyme activities and microbial carbon pump promote carbon storage by influencing bacterial communities under Nitrogen-Rich conditions in tea plantation. Agriculture 15 (3), 238. https://doi.org/10.3390/agriculture15030238.

Silatsa, F.B., Yemefack, M., Tabi, F.O., Heuvelink, G.B., Leenaars, J.G., 2020. Assessing countrywide soil organic carbon stock using hybrid machine learning modelling and legacy soil data in Cameroon. Geoderma 367, 114260. https://doi.org/10.1016/j.geoderma.2020.114260.

Simpson, E.S., Opitz, T., Wadsworth, J.L., 2023. High-dimensional modeling of spatial and spatio-temporal conditional extremes using INLA and Gaussian Markov random fields. Extremes 26 (4), 669–713. https://doi.org/10.1007/s10687-023-00468-8.

Smith, P., 2007. Land use change and soil organic carbon dynamics. Nutr. Cycl. Agrocosyst. 81 (2), 169–178. https://doi.org/10.1007/s10705-007-9138-y.

Somala, S.N., Chanda, S., AlHamaydeh, M., Mangalathu, S., 2024. Explainable XGBOOST-SHAP machine-learning model for prediction of ground motion duration in New Zealand. nat. Hazard. Rev. 25 (2). https://doi.org/10.1061/nhrefo.nhngre-1837.

Subramani, P., Dhakshnamoorthy, K., 2025. Spatio-temporal analysis and prediction by logistic regression of respiratory diseases in India. Contemp. Math. 346–366. https://doi.org/10.37256/cm.6120255923.

Sun, X., Minasny, B., Wang, H., Zhao, Y., Zhang, G., Wu, Y., 2021. Spatiotemporal modelling of soil organic matter changes in Jiangsu, China between 1980 and 2006 using INLA-SPDE. Geoderma 384, 114808. https://doi.org/10.1016/j.geoderma.2020.114808.

Szatmári, G., Pirkó, B., Koós, S., Laborczi, A., Bakaci, Z., Szabó, J., Pásztor, L., 2019. Spatio-temporal assessment of topsoil organic carbon stock change in Hungary. Soil Tillage Res. 195, 104410. https://doi.org/10.1016/j.still.2019.104410.

Takoutsing, B., Martín, J.A.R., 2025. Effects of landscape attributes on the spatial distribution of soil organic carbon stocks in smallholder farming systems in Cameroon. Catena 249, 108640. https://doi.org/10.1016/j.catena.2024.108640.

Wadoux, A.M., 2019. Using deep learning for multivariate mapping of soil with quantified uncertainty. Geoderma 351, 59–70. https://doi.org/10.1016/j.geoderma.2019.05.012.

Wadoux, A.M., Heuvelink, G.B., Lark, R.M., Lagacherie, P., Bouma, J., Mulder, V.L., Libohova, Z., Yang, L., McBratney, A.B., 2021. Ten challenges for the future of pedometrics. Geoderma 401, 115155. https://doi.org/10.1016/j.geoderma.2021.115155.

Wadoux, A.M., Molnar, C., 2022. Beyond prediction: methods for interpreting complex models of soil variation. Geoderma 422, 115953. https://doi.org/10.1016/j.geoderma.2022.115953.

Wadoux, A.M.J., Saby, N.P.A., Martin, M.P., 2023. Shapley values reveal the drivers of soil organic carbon stock prediction. Soil 9 (1), 21–38. https://doi.org/10.5194/soil-9-21-2023.

Wang, J., Deng, Y., Huang, Z., Li, D.A., Zhang, X., 2025. Identification of driving factors for heavy metals and polycyclic aromatic hydrocarbons pollution in agricultural soils using interpretable machine learning. Sci. Total Environ. 960, 178384. https://doi.org/10.1016/j.scitotenv.2025.178384.

Wang, J., Zuo, R., 2021. Spatial modelling of hydrothermal mineralization-related geochemical patterns using INLA+SPDE and local singularity analysis. Comput. Geosci. 154, 104822. https://doi.org/10.1016/j.cageo.2021.104822.

Wang, Y., Wu, P., Mei, F., Ling, Y., Qiaoy, Y., Liu, C., Leghari, S.J., Guan, X., Wang, T., 2021a. Does continuous straw returning keep China farmland soil organic carbon continued increase? A Meta-Analysis. Journal of Environmental Management 288, 112391. https://doi.org/10.1016/j.jenvman.2021.112391.

Wang, Y., Shen, X., Jiang, M., Tong, S., Lu, X., 2021b. Spatiotemporal change of aboveground biomass and its response to climate change in marshes of the Tibetan Plateau. Int. J. Appl. Earth Obs. Geoinf. 102, 102385.

Whittle, P., 1954. ON stationary processes in the plane. Biometrika 41 (3–4), 434–449. https://doi.org/10.1093/biomet/41.3.4.434.

Wiesmeier, M., Barthold, F., Blank, B., Kögel-Knabner, I., 2010. Digital mapping of soil organic matter stocks using Random Forest modeling in a semi-arid steppe ecosystem. Plant and Soil 340 (1–2), 7–24. https://doi.org/10.1007/s11104-010-0425-z.

Wolfensberger, D., Gabella, M., Boscacci, M., Germann, U., Berne, A., 2021. RainForest: a random forest algorithm for quantitative precipitation estimation over Switzerland. Atmos. Meas. Tech. 14 (4), 3169–3193. https://doi.org/10.5194/amt-14-3169-2021.

Wright, N., Newell, K., Lam, K.B.H., Kurmi, O., Chen, Z., Kartsonaki, C., 2021. Estimating ambient air pollutant levels in Suzhou through the SPDE approach with R-INLA. Int. J. Hyg. Environ. Health 235, 113766. https://doi.org/10.1016/j.ijheh.2021.113766.

Wu, C., Wu, J., Luo, Y., Zhang, L., DeGloria, S.D., 2009. Spatial prediction of soil organic matter content using cokriging with remotely sensed data. Soil Sci. Soc. Am. J. 73 (4), 1202–1208. https://doi.org/10.2136/ssaj008.0045.

Wu, W., 2021. Multiscale Bayesian Spatial Modelling and Mapping of Soil Organic Carbon with the INLA-SPDE Approach (in chinese). Zhejiang University. D Thesis.

Wu, X., Jiang, N., Li, A., Yang, Y., Cheng, H., 2025. Spatial distribution pattern of soil organic matter in the wind erosion region of northeastern China based on the

cokriging method. Catena 248, 108575. https://doi.org/10.1016/j.catena.2024.108575.

Xiao, Y., Cheng, C., Cheng, A., Kang, W., Shen, T., Yang, Q., Peng, C., Yu, L., Li, W., 2025. Dynamics of microbial communities and organic carbon pools in karst wetland soils. Catena 249, 108691. https://doi.org/10.1016/j.catena.2024.108691.

Xiong, Z., Gao, Z., Lu, J., Zhang, Y., Li, X., 2025. Straw return combined with potassium fertilization improves potassium stocks in large-macroaggregates by increasing complex iron oxide under rice-oilseed rape rotation system. Soil Tillage Res. 248, 106404. https://doi.org/10.1016/j.still.2024.106404.

Yadav, R., Huser, R., Opitz, T., 2022. A flexible Bayesian hierarchical modeling framework for spatially dependent peaks-over-threshold data. Spatial Stat. 51, 100672. https://doi.org/10.1016/j.spasta.2022.100672.

Yang, C., Yang, L., Zhang, L., Zhou, C., 2023. Soil organic matter mapping using INLA-SPDE with remote sensing based soil moisture indices and Fourier transforms decomposed variables. Geoderma 437, 116571. https://doi.org/10.1016/j.geoderma.2023.116571.

Yang, P., Hu, J., Hu, B., Luo, D., Peng, J., 2022. Estimating soil organic matter content in desert areas using in situ hyperspectral data and feature variable selection algorithms in southern Xinjiang, China. Remote Sens. 14 (20), 5221. https://doi.org/10.3390/rs14205221.

Ye, H., Huang, W., Huang, S., Huang, Y., Zhang, S., Dong, Y., Chen, P., 2017. Effects of different sampling densities on geographically weighted regression kriging for predicting soil organic carbon. Spatial Stat. 20, 76–91. https://doi.org/10.1016/j.spasta.2017.02.001.

Yu, M., Wang, Y., Deng, Q., Jiang, J., Cao, N., Tang, X., Zhang, D., Yan, J., 2024. Soil acidification enhanced soil carbon sequestration through increased mineral protection. Plant and Soil. https://doi.org/10.1007/s11104-024-06608-8.

Zare, S., Abtahi, A., Shamsi, S.R.F., Lagacherie, P., 2021. Combining laboratory measurements and proximal soil sensing data in digital soil mapping approaches. Catena 207, 105702. https://doi.org/10.1016/j.catena.2021.105702.

Zeng, S., Wu, X., Wu, B., Zhou, H., Wang, M., 2023. Rapid determination of cadmium residues in tomato leaves by Vis-NIR hyperspectral and Synergy interval PLS coupled Monte Carlo method. Food Sci. Technol. (43, e113422). https://doi.org/10.1590/fst.113422.

Zhang, M., Wu, B., Zeng, H., He, G., Liu, C., Tao, S., Zhang, Q., Nabil, M., Tian, F., Bofana, J., Beyene, A.N., Elnashar, A., Yan, N., Wang, Z., Liu, Y., 2021. GCI30: a global dataset of 30 m cropping intensity using multisource remote sensing imagery. Earth Syst. Sci. Data 13 (10), 4799–4817. https://doi.org/10.5194/essd-13-4799-2021.

Zhang, Y., Dong, X., Liu, H., Gao, T., Ren, Y., Meng, Y., Teng, C., Zhang, J., 2024. A novel approach for predicting the carbon content of birch leaf litter using Fourier transform infrared (FTIR) spectroscopy. Microchem. J. 200, 110262. https://doi.org/10.1016/j.microc.2024.110262.

Zhao, H., Ning, P., Chen, Y., Liu, J., Ghaffar, S.A., Xiaohong, T., Shi, J., 2019. Effect of straw amendment modes on soil organic carbon, nitrogen sequestration and crop yield on the North-Central Plain of China. Soil Use Manag. 35 (3), 511–525. https://doi.org/10.1111/sum.12482.

Zhao, X., Liu, P., Feng, Y., Zhang, W., Njoroge, B., Long, F., Zhou, Q., Qu, C., Gan, X., Liu, X., 2022. Changes in soil physico-chemical and microbiological properties during natural succession: a case study in lower subtropical China. Front. Plant Sci. 13. https://doi.org/10.3389/fpls.2022.878908.

Zhao, X., Ji, J., Wu, J., 2025. Effects of adding livestock and poultry manure to field strips of straw on soil organic carbon components. Appl. Sci. 15 (2), 577. https://doi.org/10.3390/app15020577.

Zhong, R., Amaral, A.V.R., Moraga, P., 2025. Spatial data fusion adjusting for preferential sampling using integrated nested Laplace approximation and stochastic partial differential equation. J. R. Stat. Soc. A. Stat. Soc. https://doi.org/10.1093/jrsss/qnae058.

Zhou, Y., Zhao, X., Guo, X., Li, Y., 2022. Mapping of soil organic carbon using machine learning models: Combination of optical and radar remote sensing data. Soil Sci. Soc. Am. J. 86 (2), 293–310. https://doi.org/10.1002/saj2.20371.

Zhuang, Q., He, J., Lu, Y., Ji, L., Xiao, J., Luo, T., 2010. Carbon dynamics of terrestrial ecosystems on the Tibetan Plateau during the 20th century: an analysis with a process-based biogeochemical model. Glob. Ecol. Biogeogr. 19 (5), 649–662. https://doi.org/10.1111/j.1466-8238.2010.00559.x.

Zhang, H., Ouyang, Z., Jiang, P., Li, M., Zhao, X., 2022. Spatial distribution patterns and influencing factors of soil carbon, phosphorus, and C: P ratio on farmlands in southeastern China. Catena 216, 106409. https://doi.org/10.1016/j.catena.2022.106409.

Chen, Y., Xi, H., Cheng, W., 2025. Spatiotemporal dynamics of soil organic carbon in desert region and its response to climate change: a case study of Heihe river. Ecol. Ind. 170, 113005. https://doi.org/10.1016/j.ecolind.2024.113005.

Wu, X., Zeng, S., Fu, H., Wu, B., Zhou, H., Dai, C., 2023. Determination of corn protein content using near-infrared spectroscopy combined with A-CARS-PLS. Food Chem.: X 18, 100666. https://doi.org/10.1016/j.fochx.2023.100666.

Wiesmeier, M., Urbanski, L., Hobley, E., Lang, B., von Lützow, M., Marin-Spiotta, E., Kögel-Knabner, I., 2019. Soil organic carbon storage as a key function of soils-a review of drivers and indicators at various scales. Geoderma 333, 149–162. https://doi.org/10.1016/j.geoderma.2018.07.026.