Shi et al. (2026). Diagnosing global soil organic carbon mapping discrepancies with interpretable machine learning: A comparison of soil type linkage and digital soil mapping

McBratney et al., 2003; Ma et al., 2019

Global soil property datasets are ultimately derived from soil survey

data and produced using two primary approaches (

Dai et al., 2019

): (i)

digital soil mapping (DSM), which predicts continuous soil properties

based on the relationship between environmental covariates and ob-

servations (

), and (ii) the soil type

linkage (STL) method, which assigns soil attribute values to map units

through polygon linkage or taxotransfer rules (

Shangguan et al., 2014;

Batjes, 2016

). In this study, we compare SoilGrids 2.0 (SG;

Poggio, et al.,

2021

), a DSM-based product, and the Global Soil Dataset for Earth

System Modelling (GSDE,

Shangguan et al., 2014

), a STL dataset, as

representative exemplars of these contrasting approaches. We selected

SG and GSDE because they provide global coverage, are widely applied

in Earth system and land-management research (

Zhang et al., 2023,

Huang et al., 2022; Xu et al., 2023

), and can be harmonised to consistent

depths and spatial resolutions for a robust, paradigm-level comparison.

While SG and GSDE provide multiple soil properties (e.g., sand, clay,

bulk density), this study focuses on SOC as the primary variable due to

its ecological significance. Using both topsoil (0

–

5 cm) and subsoil

(30

–

60 cm) SOC in Southwest China as focal depths, this study explores

the practical significance of differences between SG and GSDE and

clarifies their implications for cautious dataset interpretation and

regional-scale applications. Although the analysis is limited to South-

west China, the area's complex terrain and strong climate gradients

provide a rigorous test case, with the resulting patterns offering insights

potentially relevant to other environmentally heterogeneous regions

globally. The purpose of focusing on a regional case study is not to

develop a region-specific product, but to provide a diagnostic window

for understanding how global soil datasets behave under complex

environmental conditions.

In addition to the main comparison between SG and GSDE, we

included CSDLv2 as an auxiliary national DSM-based comparison for

Southwest China. This analysis was designed as a sensitivity test to

examine whether similar SOC

–

environment response patterns can be

observed in a China-specific gridded soil product developed under a

regional data context. For national-scale SOC mapping across China,

both CSDLv2 and GSDE derive from the same underlying SNSSC soil

profile dataset. Accordingly, CSDLv2 was not adopted as a fully inde-

pendent validation dataset. Instead, it was employed for auxiliary

sensitivity

analysis

to

verify

whether

DSM-consistent

SOC

–

environmental response patterns can be reproduced in China-

specific gridded soil products, while eliminating potential biases intro-

duced by disparate soil profile sources.

Most existing studies infer dominant SOC controls from mechanistic

reasoning or empirical regression models. However, traditional regres-

sion struggles to reliably capture the direction and magnitude of co-

variate effects in high-dimensional settings with strong multicollinearity

and nonlinear interactions. In contrast, machine learning models offer

greater flexibility for capturing nonlinear responses and higher-order

interactions among environmental drivers (

Zeng et al., 2024

). Among

them, Extreme Gradient Boosting (XGBoost), a gradient boosting tree

algorithm, is particularly well-suited for soil mapping applications due

to its robustness to overfitting and ability to handle high-dimensional

data (

Chi et al., 2022; Ma et al., 2021

). Despite these advantages, the

complexity of the ensemble tree model hinders its interpretability,

posing challenges for scientific validation and for integrating domain

knowledge. We applied Shapley Additive Explanations (SHAP) together

with partial dependence plots (PDPs) (

El Bilali et al., 2023

;

Kakhani

et al., 2025; Carbajal-Llosa et al., 2025

). In this study, interpretability

tools are not used merely for feature importance visualization, but as a

diagnostic framework to systematically analyse how different soil

datasets encode environmental controls on SOC (

Schweng et al., 2026

).

Specifically, it enables (i) quantitative comparison of driver importance,

(ii) identification of nonlinear response patterns and thresholds, and (iii)

diagnosis of how environmental gradients and interactions are repre-

sented under different mapping paradigms (

Reddy and Gopinath, 2025

).

Accordingly, this study aims to:

(1) compare SOC patterns between SG and GSDE at both topsoil (0

–

5

cm) and subsoil (30

–

60 cm) depths;

(2) evaluate broad product

–

profile agreement using legacy soil pro-

file observations; and.

(3) diagnose differences in how DSM- and STL-based products

represent SOC

–

environment relationships using XAI (explainable

AI) techniques.

To achieve these objectives, XGBoost models were trained separately

using SOC values from SG, GSDE, and soil profile observations as

response variables. Model behaviour was interpreted using SHAP and

partial dependence plots to compare driver importance, nonlinear re-

sponses, thresholds, and interaction structures. Soil profile observations

from China

’

s Second National Soil Survey were used as a legacy refer-

ence baseline to assess broad product

–

profile agreement and

SOC

–

environment consistency.

The framework (

Fig. 1

) provides a basis for cautious dataset inter-

pretation and regional-scale diagnostic assessment. Rather than seeking

to identify the single most suitable dataset for regional applications or to

produce a regionally optimized SOC map, this study's primary objective

is to diagnose the practical performance of global soil datasets at

regional scales. We aim to elucidate structural differences and associ-

ated uncertainties among them. Given the complexity of global soil-

environment relationships, which preclude detailed analyses across all

regions worldwide in a single study, we have selected a representative

region as a case study to enable a more focused and process-relevant

comparison.

2. Materials and methods

2.1. Overall methodological framework

The overall workflow of this study consisted of five major steps

(

Fig. 1

). The following subsections describe each step in detail.

1. Data preprocessing, including masking, cleaning, and quality control

2. Selection of environmental covariates relevant to SOC modelling

3. Model training using XGBoost with cross-validation

4. Model prediction and performance evaluation

5. Model interpretation using SHAP and partial dependence plots

(PDPs)

2.2. Study area

Southwest China is one of the country's seven major physical-

geographical regions, comprising the provinces of Sichuan, Yunnan,

Guizhou, Guangxi, and the municipality of Chongqing. The region spans

approximately 1.376 million km

2

, accounting for about 14.3% of China's

total land area (

Song et al., 2024

). It extends from 21

◦

06

′

to 34

◦

18

′

N and

from 97

◦

24

′

to 112

◦

06

′

E, stretching from the Qinling-Huaihe line in the

north to the vicinity of the Tropic of Cancer in the south. The region is

defined by a complex mosaic of mountains and basins, generating pro-

nounced topographic relief, with elevations ranging from several hun-

dred meters in the Sichuan Basin to 1

–

2 km across the Yunnan-Guizhou

Plateau (

Wu et al., 2023

). This diverse terrain, characterized by moun-

tain ranges, deep gorges, valleys, hills, and intermontane basins, drives

intense lateral and vertical pedological differentiation, resulting in

highly heterogeneous SOC patterns (

Zhang et al., 2008

). The combina-

tion of complex relief, monsoonal climate regimes, and heterogeneous

land use makes Southwest China a suitable region for evaluating how

different soil mapping paradigms represent coupled pedogenic and

environmental controls. However, the region was not selected to sta-

tistically represent global SOC conditions. Instead, it was used as a

diagnostic testbed because its strong topographic relief, climate gradi-

ents, and mountain

–

basin transitions provide a stringent setting for

examining how global SOC products encode terrain

–

climate controls.

G. Shi et al.

Geoderma

473

(2026)

117998

2

2.3. Data sources

This study used four primary data sources: (i) SoilGrids 2.0 (

Poggio

et al., 2021

), (ii) the Global Soil Dataset for Earth System Modelling

(

Shangguan et al., 2014

), (iii) soil profile observations from China's

Second National Soil Survey, and (iv) CSDLv2 (

Shi et al., 2025

).

SoilGrids 2.0 (SG)

. Developed by ISRIC-World Soil Information, SG

is a global digital soil mapping product generated using machine

learning (

Hengl et al., 2017

). The models were trained on more than

230,000 harmonized soil profile observations curated in the World Soil

Information Service (WoSIS) together with over 400 environmental

covariates, including climate, topography, vegetation, and lithology. SG

provides global predictions of key soil properties, such as SOC, bulk

density, particle-size fractions, and coarse fragments, at six standard

depth intervals (0

–

5, 5

–

15, 15

–

30, 30

–

60, 60

–

100, and 100

–

200 cm) at

a native resolution of 250 m (

Poggio et al., 2021

The Global Soil Dataset for Earth System Modelling (GSDE)

.

Developed by the Land-Atmosphere Interaction Research Group at Sun

Yat-sen University, GSDE integrates the FAO world soil map with mul-

tiple regional and national soil databases under a harmonized data

framework. For China, GSDE uses the 1:1,000,000 Chinese soil type map

as the spatial basis for soil property linkage. In this map, each polygon

represents a single soil type, with a minimum polygon support approx-

imately comparable to the 1-km grid scale. Soil properties are estimated

by applying taxon-based type-to-property transfer functions and linking

profile information to polygonal soil mapping units (

Shangguan et al.,

2014

). The dataset provides global soil property grids at 1-km resolu-

tion for eight depth intervals (0

–

4.5, 4.5

–

9.1, 9.1

–

16.6, 16.6

–

28.9,

28.9

–

49.3, 49.3

–

82.9, 82.9

–

138.3, and 138.3

–

229.6c m). To enable a

consistent comparison, SG was resampled to 1 km to match GSDE, and

GSDE layer-specific SOC values were harmonized to match the six

standard depth intervals used in SG using a depth-weighted averaging

approach. For each target depth interval, SOC values were computed

based on the proportional overlap between GSDE layers and the target

depth range (

Malone et al., 2009

Shangguan et al., 2013, Shi

Soil profile observations

. Soil profile reference data were taken

from the Second National Soil Survey of China (SNSSC), which collected

8,979 profiles between 1979 and 1985 (

et al., 2025

). The dataset includes site descriptions such as longitude,

latitude, elevation, slope, land cover, mean annual temperature and

laboratory measurements of SOC, total nitrogen, cation-exchange ca-

pacity, and clay content (

Shi et al., 2024

). Within Southwest China,

1,367 SNSSC soil profiles were initially identified. After removing pro-

files with missing or non-physical SOC values, invalid coordinates, or

incomplete environmental covariate information, 1,320 profiles

remained for modelling and evaluation. The final 1,320 SNSSC profiles

were randomly partitioned into training and test subsets using an 80/20

split. The reserved test subset was used as a common profile-based

reference for assessing product

–

profile agreement: SG and GSDE SOC

values were extracted at the locations of these test profiles and

compared with field-observed SOC. Separately, the SG- and GSDE-target

XGBoost surrogate models were trained and evaluated using their own

gridded samples to assess within-product surrogate fidelity. Therefore,

product

–

profile agreement and surrogate-model fidelity were inter-

preted as distinct analyses. The SNSSC profiles were used as a legacy soil

profile reference rather than as strict ground truth for present-day SOC

Fig. 1.

Research framework of the study. SG: SoilGrids 2.0; GSDE: Global Soil Dataset for Earth System Modeling; XGBoost: Extreme Gradient Boosting. The XGBoost

models for SG and GSDE are surrogate models used to approximate mapped SOC

–

environment relationships, rather than independent product validation models.

G. Shi et al.

Geoderma

473

(2026)

117998

3

conditions. Although these profiles were collected during 1979

–

1985,

they represent one of the most spatially extensive and systematically

collected soil profile datasets available for Southwest China. Therefore,

they provide a useful historical reference baseline for assessing broad

product

–

profile agreement and SOC

–

environment relationships. How-

ever, because substantial land-use and vegetation changes have

occurred since the 1980 s, comparisons with the gridded SOC products

analysed in this study should be interpreted cautiously, particularly

when assessing absolute SOC differences.

Because SoilGrids 2.0 was developed using WoSIS and other global

soil profile compilations, potential overlap between the SNSSC profiles

used in this study and the source data underlying SoilGrids cannot be

fully excluded. Therefore, the SNSSC profiles were not treated as a

strictly independent validation dataset for SG. Instead, they were used as

a legacy soil profile reference to assess broad product

–

profile agreement

and SOC

–

environment consistency.

CSDLv2 dataset.

CSDLv2 is a China-specific DSM-based gridded soil

dataset developed using Chinese soil profile observations and environ-

mental covariates. Further details on its role in this study are provided in

Section 2.4

.

Environmental covariates

. Following the SCORPAN framework

(

Minasny and McBratney, 2016

), approximately 150 environmental

covariates (

Table 1

) were compiled to represent relief, climate, vegeta-

tion, and surface properties. Topographic variables were derived from

the MERIT digital elevation model at approximately 90-meter resolution

(

Yamazaki et al., 2019

), which reduces stripe noise and bias compared to

earlier global SRTM products. Terrain metrics such as elevation standard

deviation (Elvstd) and the topographic position index (TPI) were

computed using SAGA GIS (

Conrad et al., 2015

). Biotic variables

encompassed vegetation indices from Landsat 8 Collection 2, MODIS

products such as net primary production (NPP), leaf area index (LAI),

GlobeLand30 (

Chen and Chen, 2018

), the Global Accessibility Map, and

ESA GlobCover. Climate variables were obtained from MODIS, World-

Clim v2.1, and CHELSA, including land surface temperature (LST),

precipitation, solar radiation, wind speed, and cloud-regime descriptors.

Additional surface and soil-related covariates included GlobeLand30

land cover classes and Landsat-derived albedo indices (

Baig et al., 2014

Predictor selection and data processing

. To minimize multi-

collinearity and redundancy, candidate predictors were first screened

using Pearson correlation thresholds (|r|

>

0.8) and variance inflation

factor criteria (VIF

>

5). Recursive feature elimination (RFE) was then

applied, yielding approximately 20 predictors per dataset. For cross-

product comparison, a harmonized subset of 15 key covariates was

used, with the full list provided in

Table S2

. All spatial layers were

projected into the interrupted Goode Homolosine equal-area projection

(WGS 84 datum) and resampled to a 1-km resolution. Continuous var-

iables were resampled using bilinear interpolation, while categorical

layers were resampled using nearest-neighbour assignment (

Moreira De

Sousa et al., 2019

Padarian et al., 2020; Reddy and

2.4. CSDLv2-based auxiliary sensitivity analysis

CSDLv2 was included as an auxiliary national DSM-based compari-

son for sensitivity analysis. It is a China-specific gridded soil dataset

developed using Chinese soil profile observations, including profiles

from the Second National Soil Survey collected during 1979

–

1985,

together with environmental covariates and machine-learning-based

digital soil mapping. The original CSDLv2 SOC product was produced

at 90 m spatial resolution and was harmonized to the 1-km spatial

support and corresponding depth interval used in this study before

comparison with SG and GSDE. For SOC mapping in China, CSDLv2 and

GSDE used the same underlying SNSSC profile source. Therefore,

CSDLv2 was not treated as a fully independent validation dataset.

Instead, it was used as an auxiliary sensitivity analysis to examine

whether DSM-like SOC

–

environment response patterns can also be

observed in a China-specific gridded soil product.

2.5. Model development and interpretation

Extreme Gradient Boosting (XGBoost) was adopted as the primary

modelling framework for SOC prediction. Prior to model training, re-

cords containing missing values, non-physical values or numeric

sentinel placeholders were removed to ensure consistency across data-

sets. The final dataset comprised 1,380,919 pixels from SG, 1,386,831

pixels from GSDE, and 1,320 soil profiles within Southwest China.

Model training and evaluation were conducted using the predefined

training and test subsets described in

Section 2.3

. Ten-fold cross-vali-

dation was employed for hyperparameter tuning and performance

assessment to mitigate overfitting and produce robust out-of-sample

performance estimates. It is important to note that the XGBoost

models trained on SG- and GSDE-derived SOC are not intended to serve

as independent predictive models. Instead, they are used as surrogate

models to approximate how each dataset encodes relationships between

SOC and environmental covariates. In this context, surrogate-model

performance should be interpreted as surrogate fidelity, that is, the

degree to which mapped SOC values from each product can be

approximated by the selected environmental covariates. It should not be

interpreted as an independent measure of product

–

profile agreement. In

particular, high surrogate fidelity may partly reflect spatial smoothness,

large sample size, and the degree to which a gridded product was

originally generated from environmental covariates.

2.5.1. Extreme gradient boosting (XGBoost)

XGBoost is a tree-ensemble algorithm within the gradient-boosting

paradigm, in which predictors are constructed additively, with each

Table 1

List and description of environmental covariates used in this study.

Category

Covariates

Description

Source data set

Climate factors

LSTdm

mean daytime

land surface

temperature (LST,

◦

C)

https://modis.gsfc.nasa.

gov

LSTds

standard

deviation of

daytime LST (

◦

C)

LSTd0405

mean daytime LST

in Apr-May (

◦

C)

Srad

solar radiation (kJ

m

−

2

day

−

1

)

https://worldclim.org/

data/worldclim21.html

Climatic-

moisture and

atmospheric

humidity

Wind

wind speed (m

s

−

1

)

https://worldclim.org/

data/worldclim21.html

Prec

precipitation

(mm)

MODCFma

Mean annual

cloud cover

https://www.earthenv.

org/cloud

MODCFiasd

Annual standard

deviation of cloud

cover

Vegetation and

land cover

NPP

Net Primary

Production

https://lpdaac.usgs.gov/

products/

mod17a3hgfv061/

LAI

Leaf Area Index

https://modis.gsfc.nasa.

gov

Landtypes

Land cover types

https://www.resdc.cn/

DOI/DOI.aspx?DOIID

=

54

Topographic

factors

DEM

Land surface

elevation

https://hydro.iis.u-

tokyo.ac.jp/~yamadai/

MERIT_DEM/

Elvstd

Standard

deviation of

elevation

TPI

Topographic

Position Index

Surface radiative

property

Albedo

Surface

reflectance

Calculated from Landsat

8 Collection 2 Level-2 on

the GEE platform

G. Shi et al.

Geoderma

473

(2026)

117998

4

tree correcting the residuals of the preceding ensemble. The final model

output is obtained as the sum of the contributions from all individual

trees (

Friedman, 2001

). Compared to conventional gradient-boosting

decision trees, XGBoost introduces several enhancements, including

the use of a second-order Taylor expansion of the loss for more precise

optimization, an explicit structural regularization term that penalises

model complexity (e.g., number of leaves, L2 norm of leaf weights), and

the incorporation of shrinkage, subsampling, and a sparsity-aware split

finding algorithm (

Chen and Guestrin, 2016

). Together with parallel

computation, these features improve efficiency, scalability, and gener-

alization performance in high-dimensional, multi-source modelling

contexts. In this study, SOC was modelled separately for SG, GSDE, and

soil profile datasets. Data cleaning procedures were applied uniformly

across all datasets to remove anomalous or incomplete records. The

learning objective was squared-error regression. The model's perfor-

mance was first evaluated internally via 10-fold cross-validation during

hyperparameter optimisation, and subsequently assessed on the test set.

Although Pearson correlation and VIF screening were used to reduce

multicollinearity among predictors, collinearity among DEM-derived

topographic variables cannot be fully eliminated. DEM, Elvstd, and

TPI are geometrically related because they are derived from the same

elevation surface. Therefore, their SHAP contributions should not be

interpreted as fully independent causal effects of individual terrain

metrics. Instead, they were interpreted collectively as a grouped topo-

graphic effect, representing the broader influence of terrain structure on

SOC

–

environment relationships.

2.5.2. Shapley additive explanations (SHAP) and partial dependence plots

(PDPs)

To enhance model interpretability, Shapley Additive Explanations

(SHAP) were applied to attribute machine learning model predictions to

individual environmental covariates (

Kakhani et al., 2025

). SHAP is

grounded in cooperative game theory and decomposes each prediction

into additive contributions for each feature relative to a baseline, typi-

cally the global mean prediction (

Gopinath, 2025

). This approach ensures desirable properties such as

local accuracy, symmetry, and consistency, enabling both global and

instance-level interpretation of feature influence. In the context of SOC

prediction, where environmental covariates often exhibit nonlinear in-

teractions and multicollinearity, SHAP provides a robust and theoreti-

cally grounded method for quantifying feature importance and

interaction effects. Compared to the default gain-based importance

metrics in XGBoost, SHAP offers improved comparability across datasets

and more reliable interpretation when predictors are correlated (

Meng

et al., 2020

).

Formally, the SHAP value for feature

i

at instance

x

is expressed as:

∅

i

(

x

) =

∑

S

⊆

F

\{

i

}

|

S

|

!

(

M

− |

S

| −

1

)

!

M

!

[

f

S

∪{

i

}

(

x

S

∪{

i

}

)

−

f

S

(

x

S

)

]

(1)

where

F

is the full set of

M

features,

S

is a coalition (subset) excluding

feature

i

,

f

S

is the model restricted to

S

, and

x

S

the corresponding input

vector.

PDPs were employed as a complementary diagnostic to visualize the

marginal response of SOC predictions to selected covariates by aver-

aging model predictions over the empirical distribution of the remaining

predictors. PDPs provide intuitive insights into overall response shapes,

including monotonic trends, saturation behaviour, thresholds, and

nonlinear relationships, thereby supporting interpretation of systematic

model behaviour across the covariate domain (

Kakhani et al., 2025;

Carbajal-Llosa et al., 2025

). In this study, PDPs also served to cross-

validate SHAP-based attributions, offering a model-agnostic perspec-

tive on feature importance and response patterns. Because PDPs

implicitly assume functional independence among covariates, an

assumption that is often violated in geospatial and environmental con-

texts due to strong collinearity or spatial autocorrelation, their

interpretations were evaluated in conjunction with SHAP-based attri-

butions. Agreement or divergence between the two methods was used to

assess the robustness of inferred feature effects and to identify contexts

in which marginal and conditional interpretations may differ. In this

study, SHAP and PDPs are jointly used not only to interpret model

behaviour, but to quantitatively compare how different datasets repre-

sent SOC

–

environment relationships. Specifically, these analyses are

used to (i) identify differences in the relative importance of environ-

mental drivers, (ii) characterize differences in response patterns such as

nonlinearity and thresholds, and (iii) diagnose how environmental

gradients and interactions are encoded across datasets.

To explicitly quantify pairwise interactions among environmental

predictors, SHAP interaction values were additionally calculated for the

trained XGBoost models. While standard SHAP values were used to

summarize the main contribution of individual predictors, SHAP inter-

action values were used as the primary basis for interpreting interactions

such as DEM

×

LSTdm and DEM

×

precipitation.

2.6. Evaluation criteria

Agreement was quantified using RMSE, ME, MEC, and Lin

’

s

concordance correlation coefficient (CCC). RMSE quantifies the average

magnitude of prediction errors (

Chai and Draxler, 2014

), ME captures

systematic bias (

Morley et al., 2018

), MEC provides an out-of-sample

goodness-of-fit measure, and CCC evaluates agreement with the 1:1

line by accounting for both correlation and systematic bias (

Lawrence I-

Kuei Lin, 1989

). For product

–

profile comparisons, these metrics were

used to describe broad agreement between gridded products and legacy

profile observations. For SG-, GSDE-, and CSDLv2-target surrogate

models, the same metrics were interpreted as within-target surrogate

fidelity rather than product accuracy.

MEC

=

1

−

∑

n

i

=

1

(

y

i

−

̂

y

i

)

2

∑

n

i

=

1

(

y

i

−

y

)

2

(2)

CCC

=

2

r

σ

o

σ

p

σ

2

o

+

σ

2

p

+ (

O

i

−

P

)

2

(3)

RMSE

=

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅

1

n

∑

n
i

=

1

(

̂

y

i

−

y

i

)

2

√

(4)

ME

=

1

n

∑

n
i

=

1

(

̂

y

i

−

y

i

)

(5)

where

y

i

is observed SOC,

̂

y

i

is predicted SOC,

y

the mean of observa-

tions, and

n

the number of validation samples.

P

represent the averages

of the predicted values; the standard deviations are denoted as

σ

o

and

σ

p

;

the correlation coefficient value between the predicted and observed

values is represented by

r

.

3. Results

3.1. Spatial comparison of SOC in SoilGrids 2.0 and GSDE

This study focused on two representative depth intervals, topsoil

(0

–

5 cm) and subsoil (30

–

60 cm), to maintain a focused and interpret-

able comparison.

Fig. 2

a

–

b compares the spatial distributions of SOC at

the 0

–

5 cm depth interval predicted by SG and GSDE across Southwest

China. Both datasets show a pronounced northwest-southeast gradient,

with higher SOC in the humid, high-relief mountains of the west and

northwest and lower values toward the Sichuan Basin and adjacent

southeastern lowlands. SG displays smoother gradients that are visually

aligned with topographic features, whereas GSDE shows patchier, block-

like patterns aligned with soil map-unit boundaries. A similar spatial

pattern is observed for the subsoil layer (30

–

60 cm), with SG exhibiting

smoother gradients and GSDE showing block-like structures (

Fig. S1

G. Shi et al.

Geoderma

473

(2026)

117998

5

The difference map (GSDE

−

SG;

Fig. 2

c) is dominated by negative

values, indicating that GSDE generally predicts lower SOC than SG

across most of the region. Positive anomalies are mainly confined to

northern Sichuan, forming narrow belts along the Minshan-Qionglai-

Daxue ranges. These differences weaken toward the Sichuan Basin and

remain predominantly negative across Yunnan, Guizhou, Guangxi, and

Chongqing. The largest discrepancies occur in rugged, high-elevation

terrain, whereas differences are smaller in lowlands and basins.

To illustrate local-scale morphology, a 1

◦

×

1

◦

subregion in northern

Sichuan (28

◦

45

′

–

29

◦

45

′

N, 102

◦

30

′

–

103

◦

30

′

E; red box in

Fig. 2

c) was

examined at higher resolution (

). The DEM (

a) reveals deeply

dissected ridge

–

valley topography, and the precipitation layer (

b)

shows strong orographic gradients. Within this window, SG-SOC

(

c) follows terrain-parallel gradients and shows smooth valley-

to-hillslope transitions, whereas GSDE-SOC (

d) is dominated by

coarse polygonal patches with abrupt boundaries. These local differ-

ences are consistent with the regional difference map and suggest that

SG shows more continuous environmental gradients, whereas the grid-

ded GSDE product reflects the 1:1,000,000 source soil type map, single

soil-type polygon representation, and taxon-to-property linkage used in

the Chinese part of GSDE.

3.2. Evaluation against soil profile observations

To provide a consistent legacy profile-based reference, SOC values

from SG and GSDE were compared against observed SOC from soil

profiles using a common test dataset. For the topsoil layer (0

–

5 cm), SG

showed relatively stronger agreement with the legacy profile observa-

tions than GSDE (

Fig. 4

). Data points from SG were more closely

distributed around the 1:1 line, indicating relatively stronger agreement

with the legacy profile reference in terms of both central tendency and

variability. In contrast, GSDE exhibits larger scatter and a systematic

tendency to underestimate SOC, particularly at higher values. This is

supported by evaluation metrics, with SG achieving lower RMSE and

higher model efficiency compared to GSDE. The fitted regression lines

further highlight structural differences between the two products. SG

captures a broader dynamic range of SOC values, whereas GSDE pro-

duces a more compressed distribution. However, this comparison should

be interpreted cautiously because resampling SG from 250 m to 1 km

may also affect variance and local spatial contrasts. A similar but weaker

tendency was observed for the subsoil layer (30

–

60 cm). SG showed

relatively stronger product

–

profile agreement under some metrics,

whereas GSDE tended to show lower variability and a different bias

direction (

Fig. S4

). However, the differences between the two products

were not large for all metrics, and the opposite signs of ME should be

interpreted as different bias directions rather than a substantial differ-

ence in bias magnitude. Therefore, the cross-depth comparison is

interpreted cautiously as broadly consistent tendencies rather than

definitive evidence of large performance differences.

3.3. Surrogate fidelity of XGBoost models

Because the SG-, GSDE-, and profile-target models were trained using

different response variables and sample structures, the metrics reported

in this section were interpreted as within-target surrogate fidelity rather

than cross-dataset accuracy rankings. Product-level agreement between

SG and GSDE was evaluated separately using the common soil-profile

test dataset in

Section 3.2

.

For the SG-target surrogate model, predictions on the test set showed

strong agreement with withheld SG mapped SOC values (RMSE

=

6.54 g

kg

⁻

1

; ME

= −

0.01 g kg

⁻

1

; MEC

=

0.90;

a). This high surrogate

fidelity indicates that SG-derived SOC values can be readily approxi-

mated from the selected environmental covariates, but it should not be

interpreted as independent evidence of higher dataset quality. Obser-

vations cluster closely around the 1:1 line, indicating minimal system-

atic bias across most of the SOC range. Heteroscedasticity is limited and

occurs mainly at higher SOC values. The kernel density estimates

(

d) show a slight right skew, consistent with an extended high-SOC

tail.

For the GSDE-target surrogate model, agreement between XGBoost

predictions and withheld GSDE mapped SOC values was lower (RMSE

=

12.36 g kg

⁻

1

; ME

= −

0.02 g kg

⁻

1

; MEC

=

0.604;

b). This lower

surrogate fidelity may reflect the more discrete, polygon-linked struc-

ture of the gridded GSDE product, which is constrained in China by the

1:1,000,000 source soil type map, single soil-type polygon representa-

tion, and taxon-to-property linkage, rather than necessarily indicating

lower product quality. The regression slope is lower than unity, indi-

cating underestimation at higher SOC values and overestimation at

lower SOC values. Scatter is greater than in the SG-based model and

shows distinct vertical banding, reflecting discretized SOC values

inherited from polygon-based soil mapping units.

For the profile-target model, the XGBoost predictions were compared

with field-observed SOC values from the soil profile test subset. Unlike

the SG- and GSDE-target models, this comparison reflects the learn-

ability of field-observed SOC from the selected environmental cova-

riates. The profile-target model showed considerable residual variance

relative to field-observed SOC values (RMSE

=

17.02 g kg

−

1; ME

=

0.10 g kg

−

1; MEC

=

0.389; CCC

=

0.582;

c), likely reflecting the

Fig. 2.

Spatial distribution of soil organic carbon content (SOC; g kg

⁻

1

) at the 0

–

5 cm depth interval across Southwest China: (a) SOC predictions from SoilGrids 2.0

(SG-SOC). (b) SOC predictions from the Global Soil Dataset for Earth System Modelling (GSDE-SOC). (c) Difference map (GSDE minus SG), with blue indicating areas

where GSDE is lower than SG and red indicating areas where GSDE is higher. The red box delineates the subregion selected for detailed visualization of local

SOC patterns.

G. Shi et al.

Geoderma

473

(2026)

117998

6

smaller sample size and strong site-level heterogeneity of field-observed

SOC.

Kernel density estimates (

d) further highlight distributional

differences among the three data sources. SG shows a bimodal distri-

bution with a pronounced high-value tail, GSDE is unimodal with a

lower central tendency and few values above 60 g kg

⁻

1

, and soil profiles

form a unimodal distribution centered at 20

–

25 g kg

⁻

1

. These patterns

are consistent with the spatial comparison in

Section 3.1

: SG predicts

broader and generally higher SOC distributions, whereas GSDE shows

lower central tendencies and a more truncated range. Therefore, these

metrics should be interpreted as within-target surrogate fidelity, indi-

cating how well each SOC data source can be approximated by the

selected environmental covariates, rather than as direct evidence for

ranking dataset accuracy across models. A similar pattern is observed for

the subsoil layer (30

–

60 cm), with consistent model performance across

datasets (

Fig. S2

3.4. Relationships with environmental covariates

To investigate the spatial variability of SOC and its relationships with

environmental covariates across datasets, we applied SHAP-based

attribution analysis to three SOC targets: SG predictions, GSDE pre-

dictions, and observed SOC from soil profiles.

For the soil-profile model (

a), the leading predictors were DEM,

TPI, and precipitation, followed by MODCFma, MODCFiasd, and

LSTdm. DEM and TPI capture complementary aspects of terrain, while

precipitation and cloud-cover variables showed non-negligible contri-

butions. However, because DEM, Elvstd, and TPI are derived from the

same elevation surface, their SHAP values should be interpreted as

indicating the overall importance of topographic structure rather than

Fig. 3.

Landscape-scale spatial patterns within a 1

◦

×

1

◦

subregion in northern Sichuan: (a) digital elevation model with hillshade (m); (b) mean annual precipitation

(mm); (c) SG-SOC (g kg

−

1); and (d) GSDE-SOC (g kg

−

1). A semi-transparent hillshade layer was overlaid to improve visualization of terrain-controlled

spatial patterns.

G. Shi et al.

Geoderma

473

(2026)

117998

7

independent causal effects of each terrain metric. The cloud-cover effect

should be interpreted cautiously because MODCFma may act as an in-

direct proxy for regional moisture

–

radiation conditions rather than a

direct pedogenic driver. Vegetation-related predictors played a

comparatively minor role, indicating that climatic and topographic

factors dominate SOC variability at this scale.

For the SG-target model (

b), DEM and LSTdm were the most

influential predictors, followed by precipitation. High DEM values had

mainly positive SHAP contributions, suggesting more favourable SOC

accumulation and preservation in higher terrain (

Chen et al., 2023

),

whereas high LSTdm values had negative contributions, consistent with

enhanced decomposition under warmer conditions (

Li et al., 2025; Wu

et al., 2025

). Precipitation showed a positive effect at higher levels,

consistent with SOC accumulation under humid conditions (

Chen et al.,

2023; Wang et al., 2025

For the GSDE-target model (

c), radiative and atmospheric

variables were more prominent than topographic predictors (

Wu et al.,

2025

). LSTdm, LSTd0405, srad, and cloud metrics ranked highest,

whereas DEM, Elvstd, and land cover contributed much less than in the

SG model. Although LSTdm retained a negative effect, the precipitation

response was weak and inconsistent, indicating that hydrological re-

sponses were less evident in the GSDE-target surrogate model.

PDPs were used to examine the marginal effects of DEM, LSTdm, and

precipitation across the three models (

Fig. 7

). All models showed posi-

tive DEM effects, but with different functional forms: SG displayed a

near-linear increase, GSDE showed a comparatively coarser altitudinal

response, and the profile-based model showed the strongest nonline-

arity. For LSTdm, all three models showed a nonlinear negative

response, strongest in the profile-based model. Precipitation, by

contrast, showed weak and inconsistent responses in SG and GSDE,

while the profile-based model showed a modest positive marginal effect

(

Wang et al., 2025

).

Taken together, the joint SHAP

–

PDP analysis suggests that topo-

graphic and climatic variables played important but dataset-dependent

roles in the surrogate models. Topographic variables showed strong

contributions in the profile- and SG-target models, while temperature-

related variables were consistently important across datasets and pre-

cipitation showed weaker but still informative responses (

Chen et al.,

2023

). The agreement between SHAP and PDP patterns in the profile-

based model provides a useful legacy reference perspective for inter-

preting SOC

–

environment relationships, but should not be interpreted as

strict validation of present-day SOC conditions. Compared with SG, the

GSDE-target surrogate model showed less evident responses to precipi-

tation and elevation, whereas SG represented environmental gradients

more continuously, although response magnitudes were still model-

dependent. These results suggest that differences in mapping frame-

works and source-data structures are reflected in how environmental

factors and their interactions are represented in gridded SOC products.

This highlights the value of multi-method interpretability frameworks

for identifying potential structural differences among datasets and

supporting cautious dataset interpretation.

3.5. Interaction structure of environmental factors on SOC

SHAP interaction values were used to explicitly quantify pairwise

interactions among environmental predictors, with particular focus on

DEM

×

LSTdm and DEM

×

precipitation. Standard SHAP values and

PDPs were used as supporting diagnostics for main effects and marginal

response patterns, whereas the interpretation of interaction structure

was based primarily on SHAP interaction values. The SHAP interaction

analysis revealed that DEM interacted strongly with climatic variables,

particularly LSTdm and precipitation, with

using the soil-profile

model as a legacy reference perspective. DEM emerged as an impor-

tant topographic predictor in both SG and the soil-profile model, sup-

porting the important role of topographic structure in SOC variation at

regional scales (

Minasny et al., 2024

). In contrast, DEM showed lower

importance in the GSDE-target surrogate model, suggesting weaker

terrain responsiveness in the gridded GSDE product. This pattern should

not be attributed to the STL framework alone, but may also reflect the

1:1,000,000 source-map scale, single soil-type polygon representation,

and taxon-to-property linkage used in the Chinese part of GSDE, which

can limit the representation of continuous terrain

–

SOC relationships.

The SHAP interaction analysis indicated that DEM interacted with

climatic variables, particularly LSTdm and precipitation, although the

strength and shape of these interactions differed among the profile-, SG-,

and GSDE-target models. The interaction pairs shown in

were

selected because they had high mean absolute SHAP interaction values

and represent pedologically interpretable topography

–

climate cou-

plings. DEM

×

LSTdm represents topographic modulation of thermal

conditions, whereas DEM

×

precipitation represents topographic mod-

ulation of moisture availability. The same interaction pairs were shown

for all three target models to enable direct comparison across datasets.

In the soil-profile model, the DEM-LSTdm interaction shows a clear

nonlinear pattern (

a): strongly negative values occur under high-

elevation, low-temperature conditions, whereas lower-elevation,

warmer conditions show weaker or slightly positive effects. SG

showed a broadly similar pattern in the surrogate framework, while

GSDE showed attenuated responses, particularly under low-elevation,

high-temperature conditions.

The DEM

–

precipitation interaction reveals a similarly complex top-

ography

–

hydrology coupling (

b, e, h). In the soil-profile model,

high-elevation, low-precipitation conditions are associated with nega-

tive interaction effects, whereas low-elevation, high-precipitation set-

tings showed mixed positive and negative contributions, potentially

reflecting unmeasured local hydrological, vegetation, or edaphic con-

ditions. SG showed generally weak interaction strength across the pre-

cipitation gradient, which may reflect the generalized nature of global

DSM products and their reduced sensitivity to localized environmental

combinations. GSDE showed simplified or less responsive patterns,

which may reflect the 1:1,000,000 source soil type map, single soil-type

polygon representation, and taxon-to-property linkage used in the

Fig. 4.

Comparison of mapped soil organic carbon (SOC) content at 0

–

5 cm

depth from SoilGrids (SG) and GSDE against observed SOC from soil profiles in

the test subset. Points represent individual samples, with blue and orange dots

corresponding to SG and GSDE, respectively. The dashed black line denotes the

1:1 relationship, and solid lines indicate linear regression fits. Kernel density

estimates show the distribution of mapped SOC values for each dataset.

G. Shi et al.

Geoderma

473

(2026)

117998

8

Chinese part of GSDE. This highlights the non-monotonic influence of

precipitation on SOC (

Carvalhais et al., 2014

Climate variables, particularly LSTdm and precipitation, form the

second most important predictor group (

Fig. 7

). In the soil-profile

model, LSTdm is negatively associated with SOC, whereas precipita-

tion is positively associated with SOC through its role in supporting

biomass production and soil moisture. SG broadly reproduces these

trends, whereas GSDE shows inconsistent responses, likely because

climate effects are only indirectly represented through soil classes. The

LSTdm-precipitation interaction further highlights the complexity of

climatic controls (

c, f, i). In the soil-profile model, interaction ef-

fects vary with environmental conditions, indicating a threshold-

dependent balance between moisture and temperature controls. SG

showed a broadly similar pattern but smooths extreme responses,

whereas GSDE shows much weaker interaction effects, particularly

under high-temperature conditions. Vegetation-related predictors, such

as NPP and LAI, play a comparatively minor role. In the soil-profile

model, both variables show positive contributions at high values

(

a), consistent with their role in providing organic inputs. SG

captures similar patterns (

b), whereas GSDE shows minimal

sensitivity to these variables (

c), reflecting its lack of explicit

vegetation inputs. The categorical Landtypes variable was removed from

the SHAP beeswarm plots because the continuous color scale is not

appropriate for unordered categorical variables.

Overall, the attribution analysis suggests that topographic and cli-

matic covariates were important contributors to SOC variation in the

surrogate models. SG generally showed more continuous main-effect

and interaction patterns, whereas the gridded GSDE product showed

weaker or more simplified responses. For the Chinese part of GSDE,

these differences may reflect the 1:1,000,000 source soil type map,

Fig. 5.

Surrogate-model fidelity for SOC values derived from different data sources. The x-axis represents withheld SOC values from each target data .

Source: (a) SG mapped SOC values, (b) GSDE mapped SOC values, and (c) field-observed SOC values from soil profiles. The y-axis represents the corresponding

XGBoost-predicted SOC values. Metrics describe within-target surrogate fidelity and should not be interpreted as cross-dataset accuracy rankings. (d) SOC density

curves for SG (blue), GSDE (purple), and soil profiles (orange). n denotes the number of samples in the test subset

G. Shi et al.

Geoderma

473

(2026)

117998

9

single soil-type polygon representation, and taxon-to-property linkage,

rather than the STL framework alone (

Aas et al., 2021

Arrouays et al., 2020; Bonfanti et al., 2025

4. Discussion

4.1. Contributions of covariates and their interaction effects

Soil profiles provide direct field-based observations and therefore

offer a useful legacy reference baseline for assessing broad product

–

-

profile agreement (

). How-

ever, in this study, the SNSSC profiles should not be interpreted as strict

ground truth for present-day SOC conditions because they were

collected during 1979

–

1985 and may also overlap partly with the legacy

profile information used in gridded products. Interpretation analyses

based on observed SOC from these profiles therefore provide a field-

observation-based perspective on SOC

–

environment relationships,

rather than a definitive validation of contemporary SOC conditions. This

reference perspective helps identify whether covariate

–

SOC relation-

ships inferred from gridded products are broadly consistent with legacy

profile observations, while reducing the risk of overinterpreting patterns

introduced by target-data construction or model structure. It may also

help diagnose systematic product

–

profile differences in global products

such as SG and GSDE, with potential value for cautious dataset inter-

pretation, regional bias diagnosis, and future multi-source soil mapping

strategies (

Sylvain et al., 2021

).

The combined SHAP and PDPs analyses (

Figs. 6-8

) show that

topography and climate jointly structure SOC patterns across Southwest

China. Elevation consistently emerged as an important predictor in both

the SG-based model and the profile-target model, underscoring the role

of terrain in shaping SOC distributions at regional scales. Elevation co-

varies with temperature, precipitation, and vegetation zonation,

thereby indirectly regulating both organic matter inputs and decompo-

sition rates (

Chen et al., 2023; Bonfanti et al., 2025

). These results also

demonstrate that the XAI framework can reveal how different modelling

paradigms represent environmental controls, providing insights beyond

map comparison alone.

Although MODCFma appeared as an important predictor in some

surrogate models, this result should not be interpreted as evidence that

cloud cover directly controls SOC formation. Cloud cover is more

plausibly interpreted as a proxy for regional moisture

–

radiation condi-

tions, monsoonal cloudiness, vegetation productivity, and temperature

regimes. These factors can indirectly influence SOC through plant car-

bon inputs and decomposition rates. Therefore, MODCFma was treated

as a proxy variable rather than a direct pedogenic signal.

Despite this overall consistency, the magnitude and functional form

of elevation effects differ across mapping approaches. The soil profile-

based model shows strong nonlinear responses, especially above mid-

elevation thresholds, whereas the SG-based model captures a more

Fig. 6.

SHAP beeswarm plots showing the contribution of environmental predictors to SOC predictions: (a) soil profiles (b) SG; (c) GSDE. The descriptions of

predictors are given in Table 1. Predictors are ordered by mean absolute SHAP value. Individual points represent samples, with color representing the feature values

from low (blue) to high (red). Bars overlaid on each feature indicate mean absolute SHAP values.

Fig. 7.

Partial dependence plots (PDPs) showing the marginal effects of elevation (DEM), mean daytime land-surface temperature (LSTdm), and precipitation (prec)

on SOC. Black, blue, and red curves represent PDPs derived from soil profiles, SoilGrids 2.0 (SG-OC), and GSDE (GSDE-OC), respectively. Light-red histograms (right

y-axes) indicate the frequency distributions of covariate values used for model training.

G. Shi et al.

Geoderma

473

(2026)

117998

10

continuous but attenuated response and the GSDE-based model shows

sensitivity mainly at coarse altitudinal scales. Similar contrasts are

evident for temperature and precipitation. In the soil-profile model,

temperature effects are strongly modulated by elevation and moisture

availability, and precipitation shows heterogeneous interaction effects

that reflect hydrological and geomorphic controls. By comparison, SG

reproduces many dominant marginal trends but shows weaker interac-

tion structures, whereas GSDE exhibits the weakest and least consistent

interactions, especially for elevation

–

precipitation and temper-

ature

–

precipitation couplings.

Vegetation-related predictors contributed positively at high levels in

the soil profile-based and SG-based models, but ranked lower than cli-

matic and topographic variables (

Cui et al., 2025

). This suggests that, at

the spatial scale considered here, SOC stocks are more strongly con-

strained by decomposition and stabilization processes than by vari-

ability in organic matter inputs. In contrast, GSDE shows minimal

sensitivity to vegetation-related covariates, reflecting the absence of

explicit vegetation drivers in its soil class-based framework.

Taken together, these findings indicate that both mapping paradigms

capture broad-scale SOC patterns, but they differ in how coupled envi-

ronmental gradients are represented. In this study, the DSM-based SG

product tended to preserve more continuous gradients and context-

dependent interactions, whereas the gridded GSDE product analysed

here emphasized spatial consistency within taxonomic units. For the

Chinese part of GSDE, its reduced fine-scale environmental responsive-

ness may reflect the combined influence of the 1:1,000,000 source soil

type map, single soil-type polygon representation, and taxon-to-

property linkage.

4.2. Methodological differences between GSDE and SoilGrids 2.0

The XAI-based analysis in this study serves as a diagnostic framework

for identifying structural differences between soil datasets. Rather than

focusing solely on predictive performance, the combined SHAP and PDP

results provide quantitative evidence of how different modelling stra-

tegies encode environmental controls on SOC. GSDE exemplifies the STL

paradigm, in which soil property values are assigned to polygonal soil

mapping units using taxon-based transfer rules. For China, GSDE uses

the 1:1,000,000 Chinese soil type map as the spatial basis for soil

property linkage. In this map, each polygon represents a single soil type,

with a minimum polygon support approximately comparable to the 1-

km grid scale (

Shangguan et al., 2014

). Within this framework, soil

profile information is linked to soil types and then assigned to polygonal

mapping units using type-to-property rules .

This approach capitalizes on established pedological classification

systems, ensures broad taxonomic consistency, and is computationally

efficient, making it suitable for applications that require stable and

taxonomically coherent soil inputs, such as parameterizing Earth system

models. However, because each polygon is represented by a single soil

type and assigned soil property values through taxon-to-property link-

age, the gridded GSDE product may have limited ability to represent

continuous within-polygon variation in SOC associated with contem-

porary climate, terrain, vegetation, or land-use gradients. This can lead

to block-like spatial patterns or stepped transitions at polygon bound-

aries. Therefore, for the Chinese part of GSDE analysed here, reduced

fine-scale environmental responsiveness is more appropriately inter-

preted as reflecting the 1:1,000,000 source soil type map, single soil-

type polygon representation, and taxon-to-property linkage, rather

Fig. 8.

SHAP interaction plots for key covariate pairs across three SOC models: (a-c) soil profiles, (d-f) SG, and (g-i) GSDE. Columns represent interactions between

LSTdm

and

DEM

,

prec

and

DEM

, and

LSTdm

and

prec

, respectively. Each point represents an observation. The x-axis represents the first covariate; the y-axis shows the

SHAP interaction value (with the zero line indicated). Positive values indicate synergistic interactions, while negative values indicate antagonistic interactions. Colors

encode the interacting covariate.

G. Shi et al.

Geoderma

473

(2026)

117998

11

than polygon map-unit averaging or additional aggregation from a finer

polygon map to 1 km. However, the GSDE also adopts polygon map-unit

averaging or further aggregation in other regions, including the United

States with a soil map at the 1:250,000 scale.

SG represents the quantitative DSM approach: it is trained using soil

profiles from WoSIS, integrates environmental covariates spanning

climate, topography, vegetation, geology, and hydrology, and applies

feature screening and machine-learning algorithms such as quantile

regression forests to generate continuous predictions at ~ 250 m reso-

lution, with uncertainty quantified through prediction intervals (

Poggio

et al., 2021, Lagacherie, 2025

). This approach captures environmental

gradients, particularly terrain and hydrothermal regimes, and supports

reproducible, globally consistent mapping that can be updated as new

data become available. However, in profile-sparse regions and areas

outside the covariate domain, the model performance is sensitive to

sample representativeness and geographic bias, and may be further

constrained by legacy profile errors or covariate collinearity

(

Lagacherie, 2025

). Consequently, model interpretability, predictive

stability, and uncertainty estimates remain contingent upon the quality,

representativeness, and spatial distribution of the underlying soil data.

Thus, while DSM enhances environmental responsiveness, it does not

fully resolve uncertainties arising from data sparsity and model gener-

alization at regional scales.

GSDE and SG are not competing products but complementary rep-

resentations of soil information. Within the gridded comparison and

surrogate-model framework used here, GSDE prioritizes structural

consistency and taxonomic continuity, whereas SG shows stronger

environmental responsiveness and spatial detail. In the gridded product

analysed here, GSDE showed weaker context-dependent interactions

among terrain, temperature, and moisture than SG. This pattern should

be interpreted in light of the 1:1,000,000 Chinese soil type map, single

soil-type polygon representation, and taxon-to-property linkage pro-

cedure used in the Chinese part of GSDE. These factors can limit the

representation of continuous terrain-, temperature-, and moisture-

related SOC gradients in the gridded GSDE product. SG, in contrast,

showed interaction structures more similar to those derived from the

soil-profile model. These findings suggest that future soil mapping ef-

forts may benefit from hybrid approaches that integrate the strengths of

both paradigms. For example, hybrid approaches could use soil map

units to constrain or inform DSM predictions in areas of sharp taxonomic

transitions, while allowing DSM-based fields, where supported by ob-

servations and covariates, to recover continuous gradients and reduce

artificial boundary effects.

Building on the XAI-based attribution results, several methodological

approaches for advancing soil mapping emerge. First, DSM workflows

should explicitly prioritize covariates and interactions that are strongly

aligned with pedogenic understanding (

Minasny et al., 2024

). XAI tools

such as SHAP- and PDPs-based diagnostics should be used not only for

interpretation but also to actively screen, diagnose and refine model

behaviour. Second, fusion strategies that combine traditional soil maps

with DSM outputs can improve both spatial realism and taxonomic

coherence. Examples include unit-aware smoothing, quantile mapping,

or uncertainty-weighted blending (

Lagacherie, 2025

). Third, the

continued development of harmonized, high-quality global soil profile

databases remains critical, as data representativeness ultimately con-

strains both predictive accuracy and model interpretability. These ap-

proaches may support the development of global soil products that are

more transparent, process-informed, and better characterized in terms

of uncertainty (

Ma et al., 2019

).

Compared to previous evaluation studies, which mainly focus on

spatial agreement and predictive performance, the XAI-based frame-

work adopted here provides additional insights into the mechanisms

underlying dataset differences (

El Bilali et al., 2023; Greenwell, 2017,

Kakhani et al., 2025; Carbajal-Llosa et al., 2025

). These findings are

consistent with recent studies highlighting the growing role of XAI in

soil mapping, where interpretability tools provide additional insights

into model behaviour and environmental controls beyond conventional

evaluation metrics (

Sahoo et al., 2025; Abderrahmane et al., 2025

). In

particular, it reveals how different modelling paradigms encode envi-

ronmental controls on SOC, including differences in nonlinearity,

threshold effects, and interaction structures, which are not readily

captured by conventional evaluation metrics (

Schweng et al., 2026

). It is

important to note that STL-based products remain valuable in contexts

where classification consistency, data harmonization, and stability

across regions are prioritized over fine-scale environmental

responsiveness.

Practically, the framework developed here can support DSM work-

flows in three ways. First, SHAP and PDP-based diagnostic approaches

can serve as a quality control procedure to compare datasets/products

featuring distinct environmental response mechanisms, going beyond

simple evaluation of model predictive accuracy. Second, the inferred

driver importance, response thresholds, and interaction patterns can be

checked against pedological knowledge to identify whether a product

relies on plausible soil-forming relationships or on potentially spurious

correlations. Third, such diagnostics can guide hybrid mapping and

product-selection strategies by identifying where DSM-based products

preserve continuous environmental gradients and where soil-map-based

products retain useful taxonomic consistency. In future DSM develop-

ment, these diagnostics could be incorporated alongside spatial cross-

validation and uncertainty assessment to screen covariates, evaluate

model behaviour across environmental domains, and improve the pro-

cess consistency of gridded soil products.

4.3. Influence of data-source differences on dataset comparison

A key concern in comparing SG and GSDE is that their differences

may arise not only from mapping methodology, but also from differ-

ences in underlying data sources. To explicitly assess the influence of

data provenance as part of the study design, we conducted an additional

comparative analysis using CSDLv2, a Chinese gridded soil dataset

developed primarily from nationally consistent soil profile observations

(

Shi et al., 2025

). CSDLv2 is derived from soil profiles collected within

China, primarily from the Second National Soil Survey. Importantly, for

SOC mapping in China, GSDE also relies on the same SNSSC dataset.

Therefore, for this specific variable and region, the underlying soil

profile data used by CSDLv2 and GSDE are effectively consistent. This

provides useful complementary regional context for interpreting

SG

–

GSDE differences, but it does not constitute a fully controlled com-

parison because the gridded GSDE product also reflects the 1:1,000,000

Chinese soil type map, single soil-type polygon representation, and

taxon-based transfer rules.

The CSDLv2-target surrogate model showed relatively high surro-

gate fidelity and structured responses to environmental covariates

(

Fig. 9

). These patterns provide complementary regional context for

interpreting SG

–

GSDE differences, but they should not be interpreted as

fully independent validation or as definitive evidence. For SOC mapping

in China, CSDLv2 and GSDE used the sae underlying SNSSC profile

source, and the profile-based reference was also derived from SNSSC. In

addition, the gridded GSDE product reflects the 1:1,000,000 Chinese soil

type map, single soil-type polygon representation, and taxon-to-

property linkage procedure. Therefore, the CSDLv2 analysis should be

viewed as an auxiliary sensitivity analysis rather than a controlled

comparison.

Separately from the CSDLv2 auxiliary analysis, we also extended the

SG

–

GSDE comparison to the 30

–

60 cm depth interval to examine

whether the observed differences were specific to the highly dynamic

surface layer. Although the contrast between SG and GSDE was atten-

uated relative to the 0

–

5 cm layer, the overall spatial divergence

remained evident. SG showed relatively more pronounced spatial gra-

dients and stronger associations with selected environmental covariates,

whereas GSDE retained smoother spatial patterns and lower respon-

siveness to environmental variability. The additional SHAP analysis at

G. Shi et al.

Geoderma

473

(2026)

117998

12

30

–

60 cm further showed that the relative importance and response

structure of key environmental drivers were broadly consistent with

those observed for the surface layer, although the magnitude of differ-

ences should be interpreted cautiously.

These results suggest that differences in source data may contribute

to discrepancies between products, but they may not fully account for all

observed differences. Together, these results suggest that mapping

methodology may contribute to the observed SG

–

GSDE differences, but

the effects of data provenance, source-map scale, single soil-type poly-

gon representation, and taxon-to-property linkage cannot be fully

separated in this study. This is consistent with the growing consensus in

the DSM community that model structure, covariate selection, and

machine-learning algorithms strongly influence the spatial continuity

and ecological interpretability of soil predictions (

Padarian et al., 2019;

Wadoux et al., 2020; Poggio et al., 2021

2021; Shangguan et al., 2014

4.4. Limitations and outlook

The profile-based comparison is limited by the temporal mismatch

between the SNSSC observations and the gridded SOC products analysed

in this study. Because the SNSSC profiles were collected during

1979

–

1985, before major post-1985 land-use changes, ecological

restoration, reforestation, and agricultural intensification in parts of

Southwest China, product

–

profile differences may partly reflect real

temporal SOC changes rather than mapping errors alone. Therefore, the

profiles were used as a legacy reference baseline for evaluating broad

product

–

profile agreement and SOC

–

environment relationships, not as

strict ground truth for present-day SOC conditions. Another limitation

concerns the interpretation of the CSDLv2-based comparison. For SOC

mapping in China, CSDLv2 and GSDE used the same underlying SNSSC

profile source, and the profile-based reference was also derived from

SNSSC. Therefore, CSDLv2 was not treated as a fully independent vali-

dation dataset, but was used only as an auxiliary sensitivity analysis to

provide complementary context for interpreting SOC

–

environment

response patterns, rather than as definitive evidence for separating data-

source effects from mapping-method effects.

Another source of uncertainty is the resampling of SG from its native

250 m resolution to 1 km to match GSDE. Because SOC is spatially

heterogeneous, bilinear resampling may smooth local variability,

suppress extremes, and reduce spatial variance before comparison.

Therefore, differences between SG and GSDE should be interpreted in

light of spatial support and resampling effects, in addition to differences

in mapping methodology. This limitation is particularly relevant when

interpreting variance, local gradients, and fine-scale spatial contrasts

between the two products.

This study compared two mainstream global soil products, GSDE and

SG, in their ability to represent the spatial distribution of topsoil (0

–

5

cm) and subsoil (30

–

60 cm) SOC in Southwest China (

Poggio et al.,

). The analysis revealed nonlinear re-

sponses and threshold effects of key environmental drivers, providing

useful information for interpreting SOC map discrepancies and selecting

appropriate soil datasets for regional applications (

Lundberg and Lee,

2017

). Despite these insights, several limitations arising from design

choices and methodological similarities point to important directions for

future research. First, the vertical analysis was intentionally limited to

two representative depth intervals, topsoil (0

–

5 cm) and subsoil (30

–

60

cm), to maintain a focused and interpretable comparison between SG

and GSDE. Although this design choice enabled a depth-based robust-

ness check, it does not characterize the full vertical SOC profile or

profile-integrated carbon stocks. Future work should extend the analysis

to all available standard depth intervals and explicitly evaluate SOC

stock estimates. Second, methodological singularity may introduce

interpretability bias. The XGBoost-SHAP framework efficiently uncovers

nonlinear effects, but SHAP values, which estimate feature contributions

while assuming that other variables remain unchanged, remain sensitive

to feature collinearity and spatial autocorrelation (

Aas et al., 2021;

Lundberg and Lee, 2017

). Although multicollinearity was reduced

through correlation filtering and variance inflation factor (VIF)

screening, residual dependence among environmental covariates cannot

be fully eliminated. Therefore, SHAP results should be interpreted as

indicative of relative importance and interaction patterns rather than

definitive causal relationships. Similarly, partial dependence plots rely

on the unrealistic assumption that all other variables remain fixed,

which rarely holds in complex geographic systems (

Apley and Zhu,

2020

). Employing multi-model ensembles and complementary inter-

pretability methods, such as accumulated local effects, conditional

permutation importance, and interaction screening metrics, would

improve robustness and cross-validation of inferred relationships. Third,

Fig. 9.

XGBoost surrogate modelling and SHAP interpretation for CSDLv2-derived SOC at 0

–

5 cm depth. (a) Comparison between CSDLv2 mapped SOC and XGBoost-

predicted SOC for the test subset. Points are density-shaded; the black dashed line indicates the 1:1 line, and the red line indicates the fitted regression line. Marginal

histograms show the SOC distributions, and inset statistics report RMSE, ME, MEC, and n. (b) SHAP beeswarm plot for the CSDLv2-target surrogate model. Predictors

are ranked by mean absolute Shapley value, and point colors indicate feature values from low to high.

G. Shi et al.

Geoderma

473

(2026)

117998

13

spatial uncertainty was not explicitly analysed in this study. Although

SoilGrids provides spatial uncertainty layers (

Poggio et al., 2021

),

comparable uncertainty estimates are not consistently available for

GSDE and CSDLv2, making a balanced cross-product uncertainty com-

parison difficult. Therefore, the present analysis focused on product

–

-

profile agreement and SOC

–

environment response structures rather than

formal uncertainty propagation. This omission limits our ability to

assess the spatial reliability and confidence of the mapped SOC differ-

ences. Future work should integrate available uncertainty layers,

particularly those from SoilGrids, and develop comparable uncertainty

estimates for STL- and regional DSM-based products.Fourth, although

SHAP effectively highlights dominant drivers and nonlinear responses

(

Aas et al., 2021; Wadoux and Molnar, 2022

), different XAI methods can

provide different perspectives on model behaviour, and no single

explainer should be treated as a definitive causal interpretation.

The validity and reliability of any single explainer are context-

dependent; beyond its intended scope, explanations may become un-

stable or misleading. These shortcomings highlight a critical challenge

in large-scale digital soil mapping: global or semi-global models like SG

and GSDE, while valuable for broad-scale estimation, may lack the

process fidelity needed to resolve landscape- to regional-scale nonlinear

interactions shaped by complex terrain and monsoonal climate regimes.

Guided by these limitations, future work should focus on three di-

rections directly related to the scope of this study. First, the present

analysis should be extended to all available standard depth intervals to

characterize the full vertical SOC profile. Where harmonized bulk den-

sity and coarse-fragment information are available, this could further

support comparisons of SOC stock estimates rather than SOC concen-

tration alone (

Batjes, 2016

). Second, the robustness of SHAP-based in-

terpretations should be tested using complementary interpretability

methods and spatially explicit validation strategies (

Malone et al.,

2014

). Third, spatial uncertainty should be explicitly incorporated,

particularly by using SoilGrids uncertainty layers and by developing

comparable uncertainty estimates for STL- and regional DSM-based

products. These extensions would help assess whether the

SOC

–

environment relationships identified here remain stable across

depth, scale, model structure, and uncertainty conditions (

Lilburne

et al., 2024; Minasny et al., 2024; Poggio et al., 2021

org/10.1016/j.geoderma.2026.117998

Although interpretable machine learning provides valuable insights

into model behaviour, explanations remain context-dependent and

should be interpreted as diagnostic tools rather than definitive repre-

sentations of causal mechanisms. Integrating XAI-based diagnostics with

targeted field sampling and process-based modelling would further

strengthen the physical grounding of SOC mapping efforts. It should be

noted that this study does not aim to develop new mechanistic theories

of SOC formation or to improve predictive models. Instead, its contri-

bution lies in providing a quantitative and interpretable framework for

diagnosing structural differences among existing soil datasets.

5. Conclusions

Using topsoil SOC (0

–

5 cm) as the primary focus, with an additional

30

–

60 cm analysis as a supplementary depth check, this study evaluated

how two global soil products, SG and GSDE, represent SOC spatial

patterns in Southwest China. Machine learning models, interpreted with

XAI tools, were evaluated with reference to legacy soil profile observa-

tions from China

’

s Second National Soil Survey. Both products repro-

duced the broad regional pattern of higher SOC in the western

mountains and lower values toward the Sichuan Basin. However,

profile-referenced diagnostics revealed systematic deviations, particu-

larly in humid high-elevation regions and across sharp orographic

transitions. In the legacy profile-referenced comparison, SG showed

lower errors and better agreement with the observed variability than

GSDE, whereas GSDE tended to show lower SOC values and reduced

variance.

XAI analyses identified topography and climate as important sources

of divergence between the two products. Temperature-related variables

generally showed negative associations with SOC, while precipitation

showed positive or nonlinear associations. Interaction effects were most

evident along coupled gradients of energy, water, and terrain. In these

settings, elevation modulated thermal and moisture-related responses.

SG reproduced more continuous SOC

–

environment response patterns

that were broadly consistent with the legacy profile-target model,

whereas the gridded GSDE product showed weaker representation of

continuous environmental gradients and interaction structures. For the

Chinese part of GSDE, this pattern may reflect the 1:1,000,000 source

soil type map, single soil-type polygon representation, and taxon-to-

property linkage.

From a practical perspective, the choice between SG and GSDE

should depend on the intended application. Within the 1-km gridded

comparison conducted here, SG may be more suitable for studies

requiring continuous representation of landscape- to regional-scale SOC

gradients, environmental covariate responses, and process-oriented

interpretation. In contrast, GSDE may remain useful for applications

that prioritize taxonomic consistency, structural stability, and compat-

ibility with legacy soil maps, such as Earth system modelling or large-

scale harmonized datasets. These differences should be interpreted as

design trade-offs between DSM- and STL-based products, rather than as

a universal ranking of product quality.

Overall, this study demonstrates that anchoring XAI analyses to soil

profiles can help identify where global SOC products differ and diagnose

how structural assumptions influence the representation of

SOC

–

environment relationships. The integrated workflow combining

profile-based reference comparison, feature attribution, marginal re-

sponses, and interaction screening provides a transferable framework

for cautious dataset interpretation, regional bias diagnosis, and the

development of hybrid soil mapping strategies that combine the taxo-

nomic coherence of traditional soil maps with the environmental

responsiveness of digital soil mapping. Extending this approach to

additional depths, soil properties, spatial scales, and comparable un-

certainty information would further strengthen process-oriented un-

derstanding of global soil property products.

CRediT authorship contribution statement

Gaosong Shi:

Writing

–

original draft, Visualization, Validation,

Resources, Methodology, Data curation, Conceptualization.

Zeyu Song:

Investigation, Formal analysis, Data curation.

Wenye Sun:

Methodol-

ogy, Formal analysis, Data curation.

Zhongwang Wei:

Supervision,

Software, Investigation.

Yanqi Liu:

Methodology, Investigation, Data

curation.

Naikang Xu:

Formal analysis, Data curation.

Jingyi Wang:

Methodology, Investigation.

Anne Gobin:

Resources, Methodology,

Investigation.

Qingliang Li:

Investigation, Formal analysis, Conceptu-

alization.

Wei Shangguan:

Writing

–

review

&

editing, Supervision,

Resources, Project administration, Investigation, Funding acquisition.

Yongjiu Dai:

Supervision, Resources, Project administration, Funding

acquisition, 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

The study was partially supported by the Fundamental and Inter-

disciplinary Disciplines Breakthrough Plan of the Ministry of Education

of China (JYB2025XDXM902), National Natural Science Foundation of

China (42375144, 42275155), Guangdong Basic and Applied Basic

Research Foundation (2021B0301030007) and China Meteorological

Administration Youth Innovation Team (CMA2024QN01).

G. Shi et al.

Geoderma

473

(2026)

117998

14

Appendix A. Supplementary data

Supplementary data to this article can be found online at

https://doi.

https://doi.org/10.1016/j.artint.2021.103502

Data availability

Data will be made available on request.

References

Aas, K., Jullum, M., L

ø

land, A., 2021. Explaining individual predictions when features

are dependent: more accurate approximations to Shapley values. Artif. Intell. 298,

103502.

https://doi.org/10.1111/rssb.12377

Abderrahmane, H.B., Oulad-Naoui, S., Ziani, B., 2025. Geospatial Soil Quality Analysis: A

Roadmap for Integrated Systems. https://doi.org/10.48550/arXiv.2512.09817.

Apley, D.W., Zhu, J., 2020. Visualizing the effects of predictor variables in black box

supervised learning models. J. R. Stat. Soc. Ser. B Stat Methodol. 82, 1059

–

1086.

.

Arrouays, D., McBratney, A., Bouma, J., Libohova, Z., Richer-de-Forges, A.C., Morgan, C.

L.S., Roudier, P., Poggio, L., Mulder, V.L., 2020. Impressions of digital soil maps: the

good, the not so good, and making them ever better. Geoderma Reg. 20, e00255.

https://doi.org/10.1016/j.geodrs.2020.e00255

.

Baig, M.H.A., Zhang, L., Shuai, T., Tong, Q., 2014. Derivation of a tasselled cap

transformation based on Landsat 8 at-satellite reflectance. Remote Sens. Lett. 5,

423

–

431.

https://doi.org/10.1080/2150704X.2014.915434

.

Batjes, N.H., 2016. Harmonized soil property values for broad-scale modelling (WISE30

s) with estimates of global soil carbon stocks. Geoderma 269, 61

–

68.

https://doi.

org/10.1016/j.geoderma.2016.01.034

https://doi.org/10.1016/j.geoderma.2025.117452

Bonfanti, N., Choler, P., Khedim, N., Cl

´

ement, J.-C., Barr

´

e, P., Goury, R., Baudin, F.,

C

´

ecillon, L., Saillard, A., Thuiller, W., Poulenard, J., 2025. Drivers of soil organic

carbon stocks and stability along elevation gradients. Geoderma 461, 117452.

.

Carbajal-Llosa, C., Barja, A., Pizarro, S., 2025. Ensemble machine learning for digital

mapping of soil pH and electrical conductivity in the andean agroecosystem of Peru.

Front. Soil Sci. 5, 1673628.

https://doi.org/10.3389/fsoil.2025.1673628

.

Carvalhais, N., Forkel, M., Khomik, M., Bellarby, J., Jung, M., Migliavacca, M.,

М

u, M.,

Saatchi, S., Santoro, M., Thurner, M., Weber, U., Ahrens, B., Beer, C., Cescatti, A.,

Randerson, J.T., Reichstein, M., 2014. Global covariation of carbon turnover times

with climate in terrestrial ecosystems. Nature 514, 213

–

217.

https://doi.org/10.5194/gmd-7-1247-2014

10.1038/nature13731

.

Chai, T., Draxler, R.R., 2014. Root mean square error (RMSE) or mean absolute error

(MAE)?

–

arguments against avoiding RMSE in the literature. Geosci. Model Dev. 7,

1247

–

1250.

https://doi.org/10.1007/s11430-

Chen, J., Chen, J., 2018. GlobeLand30: operational global land cover mapping and big-

data analysis. Sci. China Earth Sci. 61, 1533

–

1534.

018-9255-3

https://doi.org/10.1016/j.scitotenv.2022.160602

Chen, T., Guestrin, C., 2016. XGBoost: A Scalable Tree Boosting System, in: Proceedings

of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and

Data Mining. Presented at the KDD

’

16: The 22nd ACM SIGKDD International

Conference on Knowledge Discovery and Data Mining, ACM, San Francisco

California USA, pp. 785

–

794. https://doi.org/10.1145/2939672.2939785.

Chen, F., Feng, P., Harrison, M.T., Wang, B., Liu, K., Zhang, C., Hu, K., 2023. Cropland

carbon stocks driven by soil characteristics, rainfall and elevation. Sci. Total Environ.

862, 160602.

.

Chi, Y., Fan, M., Zhao, C., Yang, Y., Fan, H., Yang, X., Yang, J., Tao, J., 2022. Machine

learning-based estimation of ground-level NO

2

concentrations over China. Sci. Total

Environ. 807, 150721.

https://doi.org/10.1016/j.scitotenv.2021.150721

.

Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J.,

Wichmann, V., B

¨

ohner, J., 2015. System for automated geoscientific analyses

(SAGA) v. 2.1.4. Geosci. Model Dev. 8, 1991

–

2007.

https://doi.org/10.5194/gmd-8-

1991-2015

.

Cui, J., Xu, Y., Wang, M., Liu, A., Sun, L., Feng, X., Yang, Q., Wang, S., Liu, H., Lv, Y.,

Liu, K., 2025. Nonlinear threshold responses and spatial heterogeneity of soil organic

carbon under contrasting pedoclimatic regimes. Front. Plant Sci. 16, 1703663.

https://doi.org/10.3389/fpls.2025.1703663

.

Dai, Y., Shangguan, W., Wei, N., Xin, Q., Yuan, H., Zhang, S., Liu, S., Lu, X., Wang, D.,

Yan, F., 2019. A review of the global soil property maps for Earth system models.

Soil 5, 137

–

158.

https://doi.org/10.5194/soil-5-137-2019

https://doi.org/10.1016/j.jenvman.2022.116890

El Bilali, A., Abdeslam, T., Ayoub, N., Lamane, H., Ezzaouini, M.A., Elbeltagi, A., 2023.

An interpretable machine learning approach based on DNN, SVR, Extra tree, and

XGBoost models for predicting daily pan evaporation. J. Environ. Manage. 327,

116890.

https://doi.org/10.1016/bs.

Evangelista, S.J., Field, D.J., McBratney, A.B., Minasny, B., Ng, W., Padarian, J.,

Dobarco, M.R., Wadoux, A.M.C., 2024. Soil security

—

strategizing a sustainable

future for soil. Adv. Agron. 183, 1

–

70.

agron.2023.10.001

.

Friedman, J.H., 2001. Greedy function approximation: a gradient boosting machine. The

Annal. Stat. 29.

https://doi.org/10.1214/aos/1013203451

.

Greenwell, B., M., 2017. pdp: An R Package for Constructing Partial Dependence Plots.

The R Journal 9, 421. https://doi.org/10.32614/RJ-2017-016.

Grinand, C., Maire, G.L., Vieilledent, G., Razakamanarivo, H., Razafimbelo, T.,

Bernoux, M., 2017. Estimating temporal changes in soil carbon stocks at ecoregional

scale in Madagascar using remote-sensing. Int. J. Appl. Earth Obs. Geoinf. 54, 1

–

14.

https://doi.org/10.1016/j.jag.2016.09.002

.

Hengl, T., Mendes De Jesus, J., Heuvelink, G.B.M., Ruiperez Gonzalez, M., Kilibarda, M.,

Blagoti

´

c, A., Shangguan, W., Wright, M.N., Geng, X., Bauer-Marschallinger, B.,

Guevara, M.A., Vargas, R., MacMillan, R.A., Batjes, N.H., Leenaars, J.G.B.,

Ribeiro, E., Wheeler, I., Mantel, S., Kempen, B., 2017. SoilGrids250m: Global

gridded soil information based on machine learning. PLoS One 12, e0169748.

https://doi.org/10.1371/journal.pone.0169748

.

Huang, S., Eisner, S., Haddeland, I., Tadege Mengistu, Z., 2022. Evaluation of two new-

generation global soil databases for macro-scale hydrological modelling in Norway.

J. Hydrol. 610, 127895.

https://doi.org/10.1016/j.jhydrol.2022.127895

https://doi.org/10.1111/ejss.70071

Kakhani, N., Taghizadeh-Mehrjardi, R., Omarzadeh, D., Ryo, M., Heiden, U.,

Scholten, T., 2025. Towards explainable AI : interpreting soil organic carbon

prediction models using a learning-based explanation method. European J Soil

Science 76, e70071.

.

Keesstra, S.D., Bouma, J., Wallinga, J., Tittonell, P., Smith, P., Cerd

`

a, A.,

Montanarella, L., Quinton, J.N., Pachepsky, Y., van der Putten, W.H., Bardgett, R.D.,

Moolenaar, S., Mol, G., Jansen, B., Fresco, L.O., 2016. The significance of soils and

soil science towards realization of the United Nations sustainable development goals.

Soil 2, 111

–

128.

https://doi.org/10.5194/soil-2-111-2016

Lagacherie, P., 2025. Operational digital soil mapping: achievements, challenges and

future strategies to go beyond. Eur. J Soil Sci. 76, e70139.

https://doi.org/10.1111/

ejss.70139

.

Lawrence I-Kuei Lin, 1989. A Concordance Correlation Coefficient to Evaluate

Reproducibility. Biometrics 45, 255

–

268. https://doi.org/10.2307/2532051.

Li, H., Cao, Y., Xiao, J., Zhang, W., Wu, Y., Ali, A., Yuan, Z., 2025. Unveiling

uncertainties in soil organic carbon modeling: the critical role of climate response

functions. Environ. Model. Software 192, 106537.

envsoft.2025.106537

Lilburne, L., Helfenstein, A., Heuvelink, G.B.M., Eger, A., 2024. Interpreting and

evaluating digital soil mapping prediction uncertainty: a case study using texture

from SoilGrids. Geoderma 450, 117052.

geoderma.2024.117052

.

Lundberg, S., Lee, S.-I., 2017. A Unified Approach to Interpreting Model Predictions.

https://doi.org/10.48550/arXiv.1705.07874.

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/

ejss.12790

.

Ma, M., Zhao, G., He, B., Li, Q., Dong, H., Wang, S., Wang, Z., 2021. XGBoost-based

method for flash flood risk assessment. J. Hydrol. 598, 126382.

10.1016/j.jhydrol.2021.126382

https://doi.org/10.1016/j.geoderma.2009.10.007

Malone, B.P., McBratney, A.B., Minasny, B., Laslett, G.M., 2009. Mapping continuous

depth functions of soil carbon storage and available water capacity. Geoderma 154,

138

–

152.

.

Malone, B.P., Minasny, B., Odgers, N.P., McBratney, A.B., 2014. Using model averaging

to combine soil property rasters from legacy soil maps and from point data.

Geoderma 232

–

234, 34

–

44.

https://doi.org/10.1016/j.geoderma.2014.04.033

https://doi.org/10.1016/S0016-7061(03)00223-4

McBratney, A.B., Mendonça Santos, M.L., Minasny, B., 2003. On digital soil mapping.

Geoderma 117, 3

–

52.

.

Meng, Y., Yang, N., Qian, Z., Zhang, G., 2020. What makes an online review more

helpful: an interpretation framework using XGBoost and SHAP values. JTAER 16,

466

–

490.

https://doi.org/10.3390/jtaer16030029

Minasny, B., McBratney, Alex.B., 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., Bandai, T., Ghezzehei, T.A., Huang, Y.C., Ma, Y., McBratney, A.B., Ng, W.,

Norouzi, S., Padarian, J., Sharififar, A., Styc, Q., 2024. Soil science-informed

machine learning. Geoderma 452, 117094.

geoderma.2024.117094

.

Moreira De Sousa, L., Poggio, L., and Kempen, B, n.d. Comparison of FOSS4G Supported

Equal-Area Projections Using Discrete Distortion Indicatrices, ISPRS Int. J. Geo-Inf.,

8, 351, 2019.

Morley, S.K., Brito, T.V., Welling, D.T., 2018. Measures of model performance based on

the log accuracy ratio. Space Weather 16, 69

–

88.

https://doi.org/10.1002/

2017SW001669

https://doi.org/10.5194/soil-5-79-2019

Padarian, J., Minasny, B., McBratney, A.B., 2019. Using deep learning for digital soil

mapping. Soil 5, 79

–

89.

.

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/10.5194/soil-7-217-

10.5194/soil-6-389-2020

.

Poggio, L., De Sousa, L.M., Batjes, N.H., Heuvelink, G.B.M., Kempen, B., Ribeiro, E.,

Rossiter, D., 2021. SoilGrids 2.0: producing soil information for the globe with

quantified spatial uncertainty. Soil 7, 217

–

240.

2021

.

Reddy, N.L., Gopinath, M.P., 2025. Advanced deep learning framework for soil texture

classification. Sci. Rep. 15, 34407.

https://doi.org/10.1038/s41598-025-17384-5

Sahoo, S., Singha, C., Govind, A., 2025. Enhanced soil salinity index prediction using

hybrid stacking ensemble machine learning with explainable artificial intelligence

(XAI) technique: a case study of the Nile Delta. Egypt. Discov Appl Sci 8, 36.

doi.org/10.1007/s42452-025-08038-2

Schweng, S., Bernardini, L., Keiblinger, K., Kaul, H.-P., Jr., I.F., Luka

ˇ

c, N., Del Ser, J.,

Holzinger, A., 2026. What can artificial intelligence do for soil health in agriculture?

Computer Science Review 59, 100832. https://doi.org/10.1016/j.

cosrev.2025.100832.

Shangguan, W., Dai, Y., Liu, B., Zhu, A., Duan, Q., Wu, L., Ji, D., Ye, A., Yuan, H.,

Zhang, Q., Chen, D., Chen, M., Chu, J., Dou, Y., Guo, J., Li, H., Li, J., Liang, L.,

G. Shi et al.

Geoderma

473

(2026)

117998

15

Liang, X., Liu, H., Liu, S., Miao, C., Zhang, Y., 2013. A China data set of soil

properties for land surface modeling. J. Adv. Model. Earth Syst. 5, 212

–

224.

doi.org/10.1002/jame.20026

Shangguan, W., Dai, Y., Duan, Q., Liu, B., Yuan, H., 2014. A global soil data set for earth

system modeling. J. Adv. Model. Earth Syst. 6, 249

–

263.

https://doi.org/10.1002/

2013MS000293

https://doi.org/10.1016/j.geoderma.2024.116912

Shi, G., Shangguan, W., Zhang, Y., Li, Q., Wang, C., Li, L., 2024. Reducing location error

of legacy soil profiles leads to improvement in digital soil mapping. Geoderma 447,

116912.

.

Shi, G., Sun, W., Shangguan, W., Wei, Z., Yuan, H., Li, L., Sun, X., Zhang, Y., Liang, H.,

Li, D., Huang, F., Li, Q., Dai, Y., 2025. A China dataset of soil properties for land

surface modelling (version 2, CSDLv2). Earth Syst. Sci. Data 17, 517

–

543.

doi.org/10.5194/essd-17-517-2025

Song, X., Chen, H., Chen, T., Huang, Q., Deng, S., Yang, N., 2024. Spatial and temporal

variations of spring drought in Southwest China and its possible teleconnection with

the global climate events. J. Hydrol.: Reg. Stud. 51, 101655.

10.1016/j.ejrh.2024.101655

Sylvain, J.-D., Anctil, F., Thiffault,

´

E., 2021. Using bias correction and ensemble

modelling for predictive mapping and related uncertainty: a case study in digital soil

mapping. Geoderma 403, 115153.

https://doi.org/10.1016/j.jag.2016.01.005

geoderma.2021.115153

.

Vaudour, E., Gilliot, J.M., Bel, L., Lefevre, J., Chehdi, K., 2016. Regional prediction of soil

organic carbon content over temperate croplands using visible near-infrared

airborne hyperspectral imagery and synchronous field spectra. Int. J. Appl. Earth

Obs. Geoinf. 49, 24

–

38.

.

Wadoux, A.-M.-J.-C., 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

https://doi.org/10.1016/j.earscirev.2020.103359

Wadoux, A.-M.-J.-C., Minasny, B., McBratney, A.B., 2020. Machine learning for digital

soil mapping: applications, challenges and suggested solutions. Earth Sci. Rev. 210,

103359.

.

Wang, M., Zhang, S., Guo, X., Wang, G., Xia, J., Xiao, L., Luo, Z., 2025. Whole-profile soil

carbon responses to concurrent warming and precipitation changes across global

biomes. Glob. Chang. Biol. 31, e70105.

https://doi.org/10.1111/gcb.70105

.

Wu, K., Su, W., Ye, S., Li, W., Cao, Y., Jia, Z., 2023. Analysis on the geographical pattern

and driving force of traditional villages based on GIS and Geodetector: a case study

of Guizhou. China. Sci Rep 13, 20659.

https://doi.org/10.1038/s41598-023-47921-

z

Wu, Y., Li, H., Liang, X., Jiang, M., He, S., He, Y., 2025. Mechanisms behind the soil

organic carbon response to temperature elevations. Agriculture 15, 1118.

doi.org/10.3390/agriculture15111118

Xu, C., Torres-Rojas, L., Vergopolan, N., Chaney, N.W., 2023. The benefits of using state-

of-the-art digital soil properties maps to improve the modeling of soil moisture in

land surface models. Water Resour. Res. 59, e2022WR032336.

10.1029/2022WR032336

Yamazaki, D., Ikeshima, D., Sosa, J., Bates, P.D., Allen, G.H., Pavelsky, T.M., 2019.

MERIT hydro: a high-resolution global hydrography map based on latest topography

dataset. Water Resour. Res. 55, 5053

–

5073.

https://doi.org/10.1029/

2019WR024873

Zeng, T., Jin, B., Glade, T., Xie, Y., Li, Y., Zhu, Y., Yin, K., 2024. Assessing the imperative

of conditioning factor grading in machine learning-based landslide susceptibility

modeling: a critical inquiry. Catena 236, 107732.

catena.2023.107732

.

Zhang, Y., Zhao, Y.C., Shi, X.Z., Lu, X.X., Yu, D.S., Wang, H.J., Sun, W.X., Darilek, J.L.,

2008. Variation of soil organic carbon estimates in mountain regions: a case study

from Southwest China. Geoderma 146, 449

–

456.

https://doi.org/10.1016/j.agrformet.2023.109565

geoderma.2008.06.015

.

Zhang, Z., Laux, P., Baade, J., Arnault, J., Wei, J., Wang, X., Liu, Y., Schmullius, C.,

Kunstmann, H., 2023. Impact of alternative soil data sources on the uncertainties in

simulated land-atmosphere interactions. Agri. Forest Meteorol. 339, 109565.