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

HAL
open science

# National-scale mapping topsoil organic carbon of cropland in China using multitemporal Sentinel-2 images

Jie Xue, Xianglin Zhang, Songchao Chen, Zhongxing Chen, Rui Lu, Feng Liu, Bas van Wesemael, Zhou Shi

► To cite this version:

Jie Xue, Xianglin Zhang, Songchao Chen, Zhongxing Chen, Rui Lu, et al.. National-scale mapping topsoil organic carbon of cropland in China using multitemporal Sentinel-2 images. Geoderma, 2025, 456, pp.117272. 10.1016/j.geoderma.2025.117272. hal-05202899

HAL Id: hal-05202899

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

Submitted on 7 Aug 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# National-scale mapping topsoil organic carbon of cropland in China using multitemporal Sentinel-2 images

Jie Xue<sup>a,b</sup>, Xianglin Zhang<sup>c</sup>, Songchao Chen<sup>a,d,e,f</sup> <sup> ID </sup>, Zhongxing Chen<sup>d,e</sup> <sup> ID </sup>, Rui Lu<sup>d</sup>, Feng Liu<sup>f</sup>, Bas van Wesemael<sup>b</sup>, Zhou Shi<sup>a,d</sup> <sup> ID </sup>

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

## ARTICLE INFO

Handling Editor: Naoise Nunan

Keywords:

Soil organic carbon

Bare soil composite

Machine learning

Sentinel-2

Multitemporal remote sensing

## ABSTRACT

Precise monitoring of soil organic carbon (SOC) is urgently needed in agricultural regions to tackle global challenges like food security, water regulation, land degradation, and climate change. Remote sensing technology has emerged as a powerful method for detecting variations in SOC at localized scales. However, its application on a broader, national scale faces limitations, especially in countries like China, where soil landscapes exhibit significant diversity. This study aimed to couple bare soil reflectance and conventional environmental covariates to map Chinese cropland SOC content at a 10-m spatial resolution. First, a new time-series bare soil extraction method, the Two-Dimensional Bare Soil Separation Algorithm, was applied, utilizing Sentinel-2 images from 2018 to 2022 to generate a continuous spectral reflectance composite. Then, nine indices with the strongest correlation to SOC were selected. Additionally, a list of environmental covariates was prepared based on SCORPAN model. Finally, bootstrapping random forest models were fitted using the covariates selected through forward recursive feature selection (FRFS), and the spatial prediction SOC map was created. The results indicated that the framework was suitable for mapping SOC in croplands of China, with the best model using remote sensing indices and environmental covariates selected through FRFS achieving an $R^2$ of 0.62, an RMSE of 4.84 $g kg^{-1}$, and an uncertainty depicted by a 90 % prediction interval range of 17.88 $g kg^{-1}$. The final map showed that the Northeast China had the highest SOC content in cropland. Climatic conditions, position, and remote sensing indices are key covariates in national-scale SOC estimation. This study can be easily implemented across broad areas for the prediction of SOC with computational efficiency. The 10-m spatial resolution SOC map of China contributes to land management and the development of policies for precision agriculture.

## 1. Introduction

Soil organic carbon (SOC), the predominant carbon pool in terrestrial ecosystems, plays a critical role in enhancing soil fertility and agricultural productivity (Chen et al., 2024; Shi et al., 2024). However, improper management and use of agricultural land can trigger a vicious cycle of SOC loss, soil erosion, and deterioration (Van Oost et al., 2007; Wiesmeier et al., 2019). Therefore, accurate SOC monitoring is crucial, particularly for cropland soils in China, where there is an urgent challenge to feed 18 % of the global population with only 9 % of the world's cropland (Xue et al., 2024a; Zuo et al., 2018).

Digital soil mapping (DSM) utilizing the SCORPAN theory has been a popular framework for predicting spatial SOC content over the recent decade (McBratney et al., 2003; Zhang et al., 2023a). However, its dependence on static biophysical variables (e.g., parent material and relief) and soil legacy data limits the capacity to reflect current soil conditions. The agricultural practices, soil erosion, and climate shifts can change SOC contents rapidly (Chen et al., 2022). Besides, DSM relies on empirical models limiting generalizability across different areas. SOC prediction affected by the scale effect of the dominant factors, which causes the large-scale SOC map may not accurately depict variation at local scale (e.g., field scale) (Zhou et al., 2021). The existing national

<footnote><sup>a</sup> State Key Laboratory of Soil Pollution Control and Safety, Zhejiang University, Hangzhou 310058, China</footnote>

<footnote><sup>b</sup> Earth and Life Institute, Université catholique de Louvain, Louvain-la-Neuve 1348, Belgium</footnote>

<footnote><sup>c</sup> INRAE, UMR ECOSYS, Université Paris-Saclay, Palaiseau 91120, France</footnote>

<footnote><sup>d</sup> College of Environmental and Resource Sciences, Zhejiang University, Hangzhou 310058, China</footnote>

<footnote><sup>e</sup> ZJU-Hangzhou Global Scientific and Technological Innovation Center, Zhejiang University, Hangzhou 311215, China</footnote>

<footnote><sup>f</sup> State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008, China</footnote>

<footnote>* Corresponding author at: State Key Laboratory of Soil Pollution Control and Safety, Zhejiang University, Hangzhou 310058, China.</footnote>

<footnote>E-mail address: chensongchao@zju.edu.cn (S. Chen).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Spatial distribution of 1607 topsoil samples across Chinese agricultural zones.</figcaption></figure>

and global SOC maps, with coarse spatial resolutions (e.g., 90, 250, 1000 m), can only represent the spatial distribution pattern and are hardly useful for guiding precision agricultural policymaking (Liu et al., 2022a; Poggio et al., 2021; Nachtergaele et al., 2023). Hence, developing quantitative and transferable approaches for a broad-scale SOC estimation is of vital significance to sustainable land management.

Remote sensing offers an effective and economical approach for mapping and monitoring SOC at a large scale (Castaldi et al., 2023; Vaudour et al., 2019). Advances in computational processing capabilities and increasing available remote sensing data, enhancing the possibility to obtain fine-resolution global and national SOC map. Previous studies have shown that changes in SOC contents can cause differences in spectral reflectance and soil spectral modelling is considered to be a good method for predicting SOC (Shi et al., 2014; Viscarra Rossel et al., 2016). Vaudour et al. (2022) overviewed SOC prediction results from purely spectral models at national to continental scales, demonstrating that error magnitudes could reach 15 $g kg^{-1}$ or greater. Some researchers have incorporated spectral information as covariates in mixed models, effectively enhancing model performance (Simbahan et al., 2006; Pouladi et al., 2019), Li (2010) proved that using bare soil data strengthens the correlation between soil attributes and satellite data, compared to using vegetation indices. However, it is difficult to obtain a continuous soil reflectance image of the soil over a large area due to the vegetation and cloud cover. Although various authors have suggested methodologies for generating continuous bare soil composites using

time-series remote sensing data, there is a risk that the complex signals of land surface induce errors in subsequent spectral modelling (Dematté et al., 2018; Diek et al., 2017; Rogge et al., 2018). Some researchers used thresholds for spectral indices like Normalized Difference Vegetation Index (NDVI) and Normalized Burn Ratio 2 (NBR2), which can be used to mask crop residues or moist soils (Dvorakova et al., 2023). However, the single index could be prone to identify bare soil as other land cover types, such as impervious surface (Nguyen et al., 2021). Meanwhile, there are diverse crop systems and climates among the Chinses nine agricultural zones. Although Castaldi et al. (2023) tested the bare soil extraction method performance at 10 sites around the world, each site was covered by at most a single Sentinel-2 image, which is not sufficient to cover any of the nine agricultural zones in China. Besides, the cropland patterns in China are diverse and complex due to the agricultural sector is dominated by smallholder farms. Thus, using purely bare soil information combined with DSM model at a high spatial resolution could be an effective way to capture the cropland SOC variation, especially in small and fragmented field areas of Southern China.

The quantitative SOC estimation at a high spatial resolution and large scale remains a knowledge gap. Here, our work aimed to improve the accuracy of national SOC estimation by coupling remote sensing techniques and DSM framework with machine learning algorithms, and produce a SOC map at 10-m spatial resolution for Chinese cropland. The specific objectives were to: (1) establish SOC indices related to the corresponding spectral band; (2) compare the model performances by
<figure><img src="image_5.png" /><figcaption>Fig. 2. Distribution of external validation samples.</figcaption></figure>

different variables (i.e., single spectral bands, SOC indices constructed by bands combination, and environmental covariates); (3) generate a SOC map at 10-m spatial resolution for Chinese cropland and evaluate its accuracy and uncertainty.

## 2. Materials and methods

### 2.1. Study area and soil data

The study area is the entire Chinese territory, based on the principles of agricultural production conditions, characteristics, development directions, major issues, key measures, and the integrity of administrative units, the Chinese Society of Agricultural Resources and Regional Planning divides the country into nine agricultural zones (Wei et al., 2020). To optimize the SOC modelling, these nine zones are aggregated into three areas, mainly based on the agricultural situation of China (Fig. 1). Modeling in these three areas can accurately capture the primary variations in agricultural conditions across China, and would be more computationally efficient for SOC modeling than using nine agricultural zones. The first area, covering the north and east, includes the Northeast China Plain (NCP), the Loess Plateau (LP), and the Huang-Huai-Hai Plain (HHHP), which constitute the northern dryland agricultural area. The second area comprises the Middle-lower Yangtze Plain (MLYP) and Southern China (SC), representing the southern paddy field agricultural area and the southern and eastern coastal agricultural zones. Another area consists of the western agricultural zones, primarily the northwest animal husbandry and irrigation agriculture areas, i.e., the Northern arid and semiarid region (NASR). Meanwhile, due to sparse cropland, other areas such as the Qinghai-Tibet Plateau (QTP), the Sichuan Basin (SCB), and the Yunnan-Guizhou Plateau (YGP) are aggregated. Based on the Chinese Soil Taxonomy, there are about 14 primary soil orders (e.g., Ferralsols, Luvisols, Alisols, and Calcisols), which can be found throughout the three areas (Liu et al., 2022a). The soil samples were derived from a systematic national soil survey project of China, spanning from 2009 to 2019 (Liu et al., 2022a). It ensures uniform methods of soil sampling and measurement. Representative soil samples were selected from various soil-forming environments based on

the Chinese Soil Taxonomy (Cooperative Research Group on Chinese Soil Taxonomy, 2001). Prior to measuring SOC content using the Walkley-Black wet oxidation method, soil samples were air-dried at room temperature and sieved through a 2 mm mesh (Zhang and Gong, 2012).

We utilized 1607 topsoil samples (0–20 cm) from cropland, as shown in Fig. 1. The points were divided into three datasets based on three zones to construct ground SOC datasets, with 663 points belonging to zone 1, 419 points in zone 2, and 525 points in zone 3.

An external validation dataset of 1017 SOC content records was derived from the latest version of Chinese Soil Spectral Library (Shi et al., 2014) (Fig. 2). The soil samples were collected to a depth of 20 cm in cropland with sampling period between 2009 and 2018. Samples were air-dried, crushed and sieved at 2 mm. SOC content was analyzed using potassium dichromate volumetric-external heating method. Fig. 2 showed the descriptive statistics and the spatial distribution of these samples. In general, the SOC values varied from 7.60 to 25.46 g·kg⁻¹ with the mean value of 14.58 g·kg⁻¹ and the median value of 14.47 g·kg⁻¹ in the entire study area. The standard deviation (SD) was 5.41 g·kg⁻¹, with a coefficient of variation (CV) of 37.11 %, which exhibits high variability.

### 2.2. Sentinel-2A image acquisition and treatment

The Sentinel-2 Multi-Spectral Imager constellation includes 2A and 2B satellites, with a five-day revisiting period. It provides 13 spectral bands with varying spatial resolutions of 10, 20, and 60 m. For the study area, Level-2A products from 1 January 2018 to 31 December 2022 were obtained through Google Earth Engine. These products have been geometrically and radiometrically corrected. To acquire the high-quality images, the Sen2Cor algorithm was applied to detect and exclude pixels with cloud cover, cloud shadow, terrain shadow and thin cirrus. We chose 10 bands that span the visible (B2-B4), red-edge (B5-B7), NIR (B8 and B8A), and SWIR (B11-B12) spectral regions. All bands were adjusted to a 10-m spatial resolution through bilinear resampling and projected into the WGS84 Geographic Latitude-Longitude coordinate system.
<figure><img src="image_6.png" /><figcaption>Fig. 3. The workflow of this study.</figcaption></figure>

## 2.3. Cropland extent and auxiliary data

For the SOC prediction and bare soil extraction, the cropland throughout China of GlobeLand30 was used as the study area (Chen et al., 2015). The National Geomatics Center of China produced GlobeLand30, a worldwide land cover dataset featuring a spatial resolution of 30 m. It employs a pixel-object-knowledge classification approach and utilizes imagery from the Chinese Environmental and Disaster satellite as well as Landsat TM/ETM+ images. Cropland maps exist for 2000, 2010, and 2020. The accuracy of this product's cropland classification has been proven to be reliable compared with many cropland datasets (Xue et al., 2024a). Therefore, the cropland map of 2020 was selected.

## 2.4. Extraction of bare soil

We developed a new bare soil extraction method, namely *Two-Dimensional Bare Soil Separation* (TDBSS) (Xue et al., 2024b). The steps

are as follows: (1) establish a time-series database of Sentinel-2 from 2018 to 2022; (2) obtain Soil Adjusted Vegetation Index (SAVI; Eq. (1)) and Normalized Difference Red/Green Redness Index (NDRGI; eq. (2)) these two-dimensional indicators yield a scatter plot of three LULC types (i.e., “bare soil” or “crop cover” state of cropland, and impervious surface) for nine agricultural areas (Fig. 1); (3) generate a confidence ellipse to scatter plot distribution of each LULC type by density analysis; (4) introduce the Parallelepiped classification method to separate each LULC type by different confidence intervals; (5) eliminate the effects of crop residue and improve soil masking, we took the same measures as Demattê et al. (2018), using NBR2 < 0.15, B3-B2 > 0 and B4-B3 > 0; (6) eliminate open surface water by Modified Normalized Difference Water Index (mNDWI; eq. (4)) and Enhanced Vegetation Index (EVI; eq. 5) and NDVI indices with the condition if (mNDWI > EVI or mNDWI > NDVI) and (EVI < 0.1) (Zou et al., 2017); (7) composite the spectral reflectance of each soil pixel by applying a median value filter; (8) overlay the cropland extent to obtain the cropland bare soil composites (eqs. 1—5).
<table><caption>**Table 1**

Descriptive statistics for SOC (g kg⁻¹) of the sub-datasets in the three zones.</caption><thead><tr><th>Area</th><th>Dataset</th><th>N</th><th>Max</th><th>Min</th><th>Mean</th><th>Median</th><th>SD</th><th>CV (%)</th></tr></thead><tbody><tr><td rowspan="2">Zone 1</td><td>Calibration</td><td>530</td><td>60.97</td><td>0.54</td><td>11.66</td><td>9.70</td><td>7.95</td><td>68.15</td></tr><tr><td>Validation</td><td>133</td><td>43.23</td><td>1.47</td><td>11.98</td><td>9.69</td><td>8.04</td><td>67.12</td></tr><tr><td rowspan="2">Zone 2</td><td>Calibration</td><td>335</td><td>60.20</td><td>0.52</td><td>12.61</td><td>10.63</td><td>8.41</td><td>66.69</td></tr><tr><td>Validation</td><td>84</td><td>49.94</td><td>3.31</td><td>14.12</td><td>11.85</td><td>8.28</td><td>58.60</td></tr><tr><td rowspan="2">Zone 3</td><td>Calibration</td><td>420</td><td>44.95</td><td>0.33</td><td>14.50</td><td>13.32</td><td>6.47</td><td>44.64</td></tr><tr><td>Validation</td><td>105</td><td>40.58</td><td>3.04</td><td>15.13</td><td>14.38</td><td>6.78</td><td>44.78</td></tr><tr><td>Whole China</td><td>All</td><td>1607</td><td>60.97</td><td>0.33</td><td>13.27</td><td>11.54</td><td>7.79</td><td>58.72</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** (a) True color RGB 321 of the bare soil spatial patterns by TDBSS; (b) the scatter plot of the threshold definition.</figcaption></figure>

<table><caption>**Table 2**

The correlation coefficients between SOC and the selected SOC indices.</caption><thead><tr><th>Index</th><th>band i</th><th>band j</th><th>band k</th><th>band h</th><th>r</th><th>r1</th><th>r2</th><th>r3</th></tr></thead><tbody><tr><td>2-B SOC Index 7</td><td>B11</td><td>B7</td><td>/</td><td>/</td><td>0.11</td><td>0.42</td><td>0.15</td><td>-0.12</td></tr><tr><td>2-B SOC Index 8_1</td><td>B11</td><td>B6</td><td>/</td><td>/</td><td>0.25</td><td>0.50</td><td>0.38</td><td>-0.11</td></tr><tr><td>2-B SOC Index 8_2</td><td>B11</td><td>B12</td><td>/</td><td>/</td><td>0.23</td><td>0.14</td><td>0.25</td><td>0.24</td></tr><tr><td>3-B SOC Index 2</td><td>B4</td><td>B8</td><td>B3</td><td>/</td><td>-0.33</td><td>-0.43</td><td>-0.43</td><td>-0.11</td></tr><tr><td>3-B SOC Index 4</td><td>B12</td><td>B11</td><td>B4</td><td>/</td><td>0.13</td><td>0.35</td><td>0.31</td><td>-0.22</td></tr><tr><td>4-B SOC Index 1</td><td>B11</td><td>B3</td><td>B4</td><td>B12</td><td>0.14</td><td>0.36</td><td>0.32</td><td>-0.20</td></tr><tr><td>4-B SOC Index 8</td><td>B6</td><td>B4</td><td>B5</td><td>B8A</td><td>0.24</td><td>0.25</td><td>0.40</td><td>-0.02</td></tr><tr><td>4-B SOC Index 9</td><td>B8A</td><td>B5</td><td>B2</td><td>B7</td><td>-0.06</td><td>0.21</td><td>-0.17</td><td>-0.15</td></tr><tr><td>4-B SOC Index 10</td><td>B6</td><td>B2</td><td>B3</td><td>B8A</td><td>0.21</td><td>0.34</td><td>0.39</td><td>-0.07</td></tr></tbody></table>

$$SAVI = \frac{R_{nir} - R_{red}}{R_{nir} + R_{red} + 0.5} \times 1.5 \qquad (1)$$

$$NDRGI = \frac{R_{red} - R_{green}}{R_{red} + R_{green}} \qquad (2)$$

$$NBR2 = \frac{R_{swir1} - R_{swir2}}{R_{swir1} + R_{swir2}} \qquad (3)$$

$$mNDWI = \frac{R_{green} - R_{swir1}}{R_{green} + R_{swir1}} \qquad (4)$$

$$EVI = \frac{R_{nir} - R_{red}}{1 + R_{nir} + 6 * R_{red} - 7.5 * R_{blue}} \qquad (5)$$

where $R_{nir}$, $R_{red}$, $R_{green}$, $R_{swir1}$, $R_{blue}$ are the NIR, Red, Green, SWIR1 and Blue bands, respectively.

## 2.5. SOC modelling covariates

### 2.5.1. Remote sensing data

We calculated a series of remote sensing indices, including three classifications, using the bare soil reflectance data from Sentinel-2. Firstly, the reflectance from Sentinel-2 bands 2 to 8A, as well as bands 11 and 12, was used. These bands include Blue (B), Green (G), Red (R), Red Edge 1 (RE1), Red Edge 2 (RE2), Red Edge 3 (RE3), Near Infrared (NIR), Red Edge 4 (RE4), Short-Wave Infrared 1 (SWIR 1), Short-Wave Infrared 2 (SWIR 2). Secondly, we constructed new indices in different bands combination. It contains eight formulations of two bands (2-B) SOC indices, nine formulations of three bands (3-B) SOC indices, and ten formulations of four bands (4-B) SOC indices and listed in Table S1–S3. $B_i$, $B_j$, $B_k$ and $B_h$ refer to a single band from Band 2 to Band 8A, Band 11, and Band 12 from Sentinel-2. A total of 3540 combinations were calculated, and then we selected the most four important indices for each combination. Thus, there were 108 band combinations left. Finally, 13 vegetation and 3 SOC indices calculated by bare soil reflectance were introduced in the SOC prediction model. Their names and formulations
<table><caption>**Table 3**
Predictive performance of SOC prediction by remote sensing indices.</caption><thead><tr><th>Area</th><th>Covariate</th><th>R<sup>2</sup></th><th>RMSE (g kg<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="5">Zone 1</td><td>Single Bands</td><td>0.24</td><td>6.98</td></tr><tr><td>New Indices</td><td><strong>0.36</strong></td><td><strong>6.39</strong></td></tr><tr><td>Existing Indices</td><td>0.27</td><td>6.86</td></tr><tr><td>Single Bands with New Indices</td><td>0.32</td><td>6.61</td></tr><tr><td>New and Existing Indices</td><td>0.36</td><td>6.39</td></tr><tr><td rowspan="5">Zone 2</td><td>All remote sensing variables</td><td>0.32</td><td>6.60</td></tr><tr><td>Single Bands</td><td>0.35</td><td>6.62</td></tr><tr><td>New Indices</td><td>0.38</td><td>6.49</td></tr><tr><td>Existing Indices</td><td>0.39</td><td>6.43</td></tr><tr><td>Single Bands with New Indices</td><td>0.40</td><td>6.37</td></tr><tr><td rowspan="5">Zone 3</td><td>New and Existing Indices</td><td><strong>0.41</strong></td><td><strong>6.35</strong></td></tr><tr><td>All remote sensing variables</td><td>0.32</td><td>6.42</td></tr><tr><td>Single Bands</td><td>0.17</td><td>6.15</td></tr><tr><td>New Indices</td><td>0.21</td><td>5.99</td></tr><tr><td>Existing Indices</td><td>0.12</td><td>6.33</td></tr><tr><td rowspan="5">Whole China</td><td>Single Bands with New Indices</td><td>0.20</td><td>6.03</td></tr><tr><td>New and Existing Indices</td><td>0.18</td><td>6.09</td></tr><tr><td>All remote sensing Indices</td><td><strong>0.22</strong></td><td><strong>5.96</strong></td></tr><tr><td>Single Bands</td><td>0.28</td><td>6.63</td></tr><tr><td>New Indices</td><td>0.35</td><td>6.29</td></tr><tr><td rowspan="5">Whole China</td><td>Existing Indices</td><td>0.29</td><td>6.58</td></tr><tr><td>Single Bands with New Indices</td><td>0.34</td><td>6.36</td></tr><tr><td>New and Existing Indices</td><td><strong>0.35</strong></td><td><strong>6.28</strong></td></tr><tr><td>All remote sensing variables</td><td>0.32</td><td>6.42</td></tr></tbody></table>

<table><caption>**Table 4**
Predictive performance of SOC prediction by remote sensing indices with environmental covariates before and after FRFS.</caption><thead><tr><th>Area</th><th>Covariate</th><th>R<sup>2</sup></th><th>RMSE (g kg<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="2">Zone 1</td><td>All</td><td>0.67</td><td>4.63</td></tr><tr><td>FRFS</td><td><strong>0.67</strong></td><td><strong>4.59</strong></td></tr><tr><td rowspan="2">Zone 2</td><td>All</td><td>0.52</td><td>5.68</td></tr><tr><td>FRFS</td><td><strong>0.59</strong></td><td><strong>5.23</strong></td></tr><tr><td rowspan="2">Zone 3</td><td>All</td><td>0.42</td><td>5.12</td></tr><tr><td>FRFS</td><td><strong>0.49</strong></td><td><strong>4.82</strong></td></tr><tr><td rowspan="2">Whole China</td><td>All</td><td>0.58</td><td>5.08</td></tr><tr><td>FRFS</td><td><strong>0.62</strong></td><td><strong>4.84</strong></td></tr></tbody></table>

were listed in the Table S4. To reduce the redundancy of variables, we calculated the importance between the indices mentioned above and SOC. The new indices among the top thirty variables in importance will be selected for the final SOC prediction.

### 2.5.2. Other geospatial data

According to the SCORPAN theory, there are a list of potential covariates related to SOC at a national scale. We selected 68 relevant predictors, including soil, climate, organisms, relief, and position (Table S5). There are 13 soil factors, two of which are derived from early work i.e. the soil erosion (SE) (Teng et al., 2019), and NASA-USDA Enhanced SMAP Global soil moisture data (SM, Mladenova et al., 2020). The other 11 soil indices at 10 m are calculated from Sentinel 2 images using GEE, including Brightness Index (BI, Abbas and Khan, 2007), Carbonate Index (Carl, Boettinger et al., 2008), Gypsum index (Gl, Boettinger et al., 2008), Hue Index (HI, Mandal, 2016), Saturation Index (SI, Mandal, 2016), Iron Oxide (IO, Hewson et al., 2001), Reflectance Absorption Index (RAI, Ghaemi et al., 2013), and the others.

The 23 climate factors are mainly derived from the 1 km WorldClim dataset (Fick and Hijmans, 2017), including 17 BIO variables (e.g., Mean diurnal range (BIO02), Precipitation of wettest month (BIO13), etc.) and 2 climatology (i.e., mean annual temperature (MAT) and mean annual precipitation (MAP)). We also obtained Day-time and Night-time Land Surface Temperature (LSTD and LSTN) from 1 km MOD11A1 dataset distributed by NASA's Land Processes Distributed Active Archive Center (LP DAAC). Another two covariates Potential Evapotranspiration and Evapotranspiration (PET and ET) of 500 m MOD17A2H dataset

produced by LP DAAC were considered as well (Running et al., 2015).

A series of 10 m organism factors were calculated by Sentinel 2 images in GEE, such as Canopy Index (CANI, Vescovo and Gianelle, 2008), Normalized Difference Green/Red Index (NDRI, Bannari et al., 1995), Enhanced Vegetation Index (EVI, Huete et al., 2002), and six other indices. Meanwhile, the calculation of Net Primary Production (NPP) and Gross Primary Production (GPP) at 500 m were from MOD17A2H (Running et al., 2015).

We obtained the digital elevation data at 30 m resolution from SRTM of NASA JPL. Using the TAGEE tool, we derived 13 relief derivatives from the elevation data (Safanelli et al., 2020), such as Elevation (ELE), Aspect (ASP), Hillshade (HS), Slope (SLO), Northness (North), Eastness (East), Horizontal curvature (HC), Vertical curvature (VC), Mean curvature (MEC), Minimal curvature (MIC), Maximal curvature (MAC), Gaussian curvature (GC), and Shape Index (SHI).

Position factors were determined using latitude and longitude within GEE, taking the form of oblique geographic coordinates (OGC) at angles of 15°, 30°, 45°, 60°, 75°, 105°, 120°, 135°, 150°, and 165° (Møller et al., 2020).

All covariates were standardized to the CGCS WGS 1984 geographic coordinate system.

## 2.6. SOC predictive modelling

### 2.6.1. Predictive model

By leveraging Sentinel-2 and other geographic data, we used the China cropland SOC sampling data to train the bootstrapping random forest (RF) model. This data was segmented into subsets for each of the agricultural zones. We independently constructed estimation models for each zone, enabling us to calculate the SOC across these three regions. Finally, we validated the accuracy and assessed the uncertainty to produce a cropland SOC map for China at a 10-m resolution. The mainly workflow is shown in the Fig. 3.

The RF model employs a bagging ensemble learning algorithm that averages many individual regression trees to derive the final prediction (Breiman, 2001). RF enhance model robustness by mitigating overfitting through bagging and demonstrating low sensitivity to multicollinearity (Hu et al., 2024b; Zhang et al., 2024). There are two main parameters decided the model performance, i.e., the number of trees (ntree), and the number of variables in each tree (mtry). The RF model was carried out in Rstudio using the “caret” package (Kuhn, 2008), and we used parameters with ntree of 500 and mtry of square root of the number of variables. Meanwhile, we calculated the %IncMSE to present the importance of the variables.

### 2.6.2. Predictor selection

To streamline the model structure without compromising performance, we employed forward recursive feature selection (FRFS, Xiao et al., 2022) to identify the most important predictors for the RF model. This method aimed to reduce redundant information and improve the efficiency of the model. FRFS adopts forward selection strategy, though it is based on the improvement of recursive feature elimination (RFE, Guyon et al., 2002), which uses the backward selection. Therefore, FRFS can fully consider the unique information of some less importance variables instead of ignore them by the controlling variables (Zhang et al., 2024).

The key stages of FRFS are outlined below: (1) fit a model using n variables, and calculate their importance; (2) begin by choosing the most significant variable to fit the model, and then assess its performance through cross-validation; (3) integrate one of the remaining variable with the initial factor to fit the model and record the best performance among each combination; (4) repeat steps 2 and 3, incrementing the number of variables from 3 to n, until the optimal variables for model fitness are identified. In this study, the internal model was RF, and the FRFS method was performed in RStudio.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Scatter plots of the RF models using the selected covariates in (a-c) three zones and (d) the whole China.</figcaption></figure>

### 2.6.3. Model evaluation

Independent validation was employed to assess the performance of the SOC predictive model. We used a random strategy to partition the soil samples (n = 1607) into a calibration dataset (80 %, n = 1285) and a validation dataset (20 %, n = 322).

The performance of RF model was assessed by coefficients of determination (R²; eq. (6)), and root mean square error (RMSE; eq. (7)). Models achieving higher R² and lower RMSE are regarded as more effective (eqs. (6) and (7)).

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

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

where y means the measured value and ŷ represents the predicted value; ŷ and ŷ̄ are the mean of measurements and predictions, respectively; n means the sample number.

### 2.6.4. Uncertainty assessment

Twenty RF models were created using non-parametric bootstrapping to assess and quantify prediction uncertainty, allowing for the determination of prediction intervals. The bootstrapping is a repeated sam-

pling technique that randomly selects samples with replacement. Both the bootstrapping dataset and the original dataset contain the same number of samples and exhibit a similar probability distribution (Zhang et al., 2022). This procedure resulted in 20 bootstrapping predictions for each pixel. The final map was generated by averaging these 20 maps, and the uncertainty was assessed using the variance from these maps and the Mean Squared Error (MSE) obtained from validation data (Malone et al., 2016; Zhou et al., 2022). The 90 % prediction interval for each pixel was determined by taking the square root of the total prediction variance and scaling it with a z-score of 1.645, corresponding to a 90 % confidence level. These indices can be calculated using equations (8) to (12).

$$VAR_{all} = VAR_{boot} + MSE_{vali} \qquad (8)$$

$$SD = \sqrt{VAR_{all}} \qquad (9)$$

$$UPL = MEAN_{boot} + SD \times 1.645 \qquad (10)$$

$$LPL = MEAN_{boot} - SD \times 1.645 \qquad (11)$$

$$PIR = UPL - LPL \qquad (12)$$

where VARall represents the variance calculated from 20 bootstrapping predictions, MSEvali is the mean square error derived from the validation
<figure><img src="image_9.png" /><figcaption>Fig. 6. Predicted and observed SOC contents of external validation dataset (n = 1017). The red line is the 1:1 line and the blue line is the fitted line.</figcaption></figure>

data, SD is the standard deviation, MEAN\_boot refers to the average of the bootstrapping predictions. The terms UPL, LPL, and PIR denote the upper prediction limit, lower prediction limit, and the range of the prediction interval, respectively. Furthermore, we employed the Prediction Interval Coverage Probability (PICP, Malone et al., 2016) to evaluate the effectiveness of our prediction interval. This metric calculates the proportion of observed validation data that falls within the specified prediction limits. Ideally, for a 90 % prediction interval, the PICP should approximate 0.9, suggesting that the uncertainty has been accurately quantified.

# 3. Results

## 3.1. Statistics of SOC

The statistical characteristics of sub-datasets were similarly distributed across the different zones (Table 1), ensuring the efficiency of model training. In general, the SOC values varied from 0.33 to 60.97 g·kg⁻¹ with the mean value of 13.27 g·kg⁻¹ and the median value of 11.54 g·kg⁻¹ in the entire study area. The SD was 7.79 g·kg⁻¹, with a CV of 58.72 %, which shows large variability.

## 3.2. Selection of constructed indices related to topsoil

The true color RGB per-pixel reflectance composite images of bare soil by the TDBSS method are displayed in Fig. 4(a). Fig. 4(b) is a scatter plot showing the thresholds of SAVI and NDRGI in cropland with vegetation, bare soil, and impervious surfaces, using the LP agricultural zone as an example. It can be observed that the ellipses for bare soil and impervious surfaces are almost tangent, easily distinguishing these three land cover types in the two-dimensional space.

As shown in Table S1 to S3, there are 27 formulations for 10 bands combinations. Therefore, to simplify the model, the four most important combinations of each formulation were maintained in the first selection. Then, we calculated the importance between bare soil reflectance related covariates and SOC. There were nine band combinations in the top thirty importance, which were selected for the SOC prediction. The selected nine band combinations and their correlation coefficients with SOC in each zone were shown in Table 2. The absolute *r* of three zones

ranked in descending order as follows: r1 (0.14–0.50) > r2 (0.15–0.43) > r3 (0.02–0.24). It indicated the bare soil indices were more suitable in north China.

## 3.3. Predicting SOC using RF model

### 3.3.1. Modelling SOC by remote sensing indices

Table 3 lists 20 bootstrapping independent validation results for SOC predictions by combining different remote sensing indices calculated from bare soil reflectance. Model performance varied with specific combinations. For zone 1, the 9 newly constructed indices either on their own or combined with existing indices performed similarly, with R² of 0.36 and RMSE of 6.39 g·kg⁻¹. It can be seen that the combination of new and existing indices performed best in zone 2, for which R² was 0.41 and RMSE was 6.35 g·kg⁻¹. In zone 3, the highest R² (0.22) and lowest RMSE (5.96 g·kg⁻¹) were shown in using a model with all remote sensing indices. For the perspective of the entire territory, SOC was predicted with the best accuracy by the combination of new and existing indices, with R² of 0.35 and RMSE of 6.28 g·kg⁻¹, which performed almost no different from using only the new indices. Therefore, the new indices can be used on their own and thus reducing the number of covariates.

### 3.3.2. Modelling SOC by remote sensing indices with environmental covariates

Table 4 and Fig. 5 display 20 bootstrapping independent validation results for SOC predictions by integrating remote sensing indices and environmental covariates, before and after applying FRFS. It can be seen that the accuracy has improved at both regional and national scale, and the best accuracies were obtained after FRFS. Compared to the models with only remote sensing indices, the relative improvements of R²s with the inclusion of environmental covariates were 48.44–56.58 %, 30.51–39.71 %, and 34.62–44.58 % and the relative decreases of RMSEs of 25.79–37.40 %, 16.92–21.90 %, and 18.01–26.54 % from zone 1 to zone 3 separately. As for the whole China, the highest prediction accuracy achieved an R² of 0.62, with the RMSE as low as 4.84 g·kg⁻¹.

Model performance also varied with the location of zones. The prediction accuracy of SOC increased substantially with increasing latitude, i.e., the R² improved from zone 3 to zone 1. Such increase in accuracy was observed in the remote sensing indices models as well, where zone 3 displayed the worst accuracy while zone 1 and zone 2 performed better (Table 3).

### 3.3.3. Model performance of external validation dataset

For the external validation dataset, the scatter plot between the observed and predicted SOC values were plotted in Fig. 6. The R² value was 0.50, which was lower than the R² of the validation data. This is reasonable given the narrower SOC distribution of the external validation dataset compared to the one used for model construction. The small RMSE (5.42 g·kg⁻¹) also indicated a good model performance.

## 3.4. Relative importance of covariates

Our findings show that the best SOC predictions were obtained by combining remote sensing indices with environmental covariates after FRFS (R² of 0.62, RMSE of 4.84 g·kg⁻¹). The climate and position covariates account for the majority, though the FRFS selected different variables in each zone. The high relative importance of climate (42.79, 44.47 and 27.87 %), and position (28.12, 20.31 and 41.64 %) covariates in zone 1 to 3 and their strong correlations with SOC (|r|>0.6 in zone 1 and 3), indicating they contribute a lot to the SOC prediction models (Fig. 7). There are still four remote sensing indices left after variable selection by FRFS. They are SWIR1 band (B11) and the 3-B SOC Index 2 (3B_2) in zone 1, 4-B SOC Index 10 (4B_10) in zone 2, and Red band (B4) in zone 3. The relative importance of remote sensing variables in zone 1 to 3 were with 14.59, 11.08, and 10.78 %.
<figure><img src="image_10.png" /><figcaption>Fig. 7. Correlation plot among SOC and selected predictors for each zone. The predictor abbreviations match the ones detailed in Tables S1-S4.</figcaption></figure>

<table><caption>**Table 5**
Optimal variables for the RF model were selected using the FRFS algorithm.</caption><thead><tr><th>Factors</th><th>Zone 1 (N = 9)</th><th>Zone 2 (N = 12)</th><th>Zone 3 (N = 8)</th></tr></thead><tbody><tr><td>Remote Sensing</td><td>SWIR1, 3B_2</td><td>4B_10</td><td>Red</td></tr><tr><td>Soil Climate</td><td>/<br>BIO08,11,14,15</td><td>FM<br>PET,<br>BIO03,04,08,19</td><td>/<br>MAP, BIO11</td></tr><tr><td>Organisms</td><td>NDRI, CANI</td><td>/</td><td>NDRI</td></tr><tr><td>Relief Position</td><td>/<br>OGC120</td><td>East, North, ASP<br>OGC45, OGC150</td><td>GC<br>OGC75, OGC120,<br>OGC165</td></tr></tbody></table>

### 3.5. SOC map of China cropland at 10-m spatial resolution

We trained 20 bootstrapping-RF models across three regional zones using selected covariates (Table 5), and then averaged the outcomes to produce the final result. The national cropland SOC map at 10-m spatial resolution is presented in Fig. 8(a). The mean SOC content was predicted to vary from 6.72 to 32.23 g kg⁻¹. The highest value occurred in the Northeast China and Tianshan Mountains of NASR. Low SOC values were primarily found in the HHHP and SCB areas. Meanwhile, we

selected three representative regions from three agricultural zones to further refine the SOC map at the field level (Fig. 8(c, e, g)), clearly illustrating the variation in SOC content among different fields.

The PICP values were 0.95, 0.95, 0.97 and 0.95 for zone 1, zone 2, zone 3 and the whole China, separately. All of them were near 90 %, suggesting that the uncertainty estimations for the SOC predictions are reliable. Fig. 8(b) showed 90 % prediction interval map produced by the RF model based on 20 bootstrapping, presenting the uncertainty for SOC estimations. High prediction uncertainty was detected in YGP areas and Greater and Lesser Khingan Ranges. We also scaled down the uncertainty map to the field level, demonstrating that the uncertainty in southern China exceeds that in northern China (Fig. 8(d, f, h)).

## 4. Discussion

### 4.1. The potential of multitemporal Sentinel-2 remote sensing for SOC prediction

The accuracy of different remote sensing indices models proved that bare soil reflectance calculated from Sentinel-2 images had the ability for SOC predictions. Previous studies have verified that vis-NIR spectroscopy is effective in predicting most soil properties as well, especially for SOC (Xu et al., 2019; Xue et al., 2023). That primarily stems from the
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Spatial patterns of (a) average SOC content and (b) relative 90% prediction limit range across China; (c) – (h) show these two values at field scale in three zoomed regions.</figcaption></figure>

direct spectral responses to carbon and nitrogen, which are attributable
to the overtones and combinations of C-H, C-O, C=O and N-H bonds
within the vis-NIR spectrum (Zhang et al., 2023c). As demonstrated by
previous researches, some characteristic bands of C-H may be around
825, 853, and 1650 nm, which are corresponding to the B8, B8A, and
B11 bands of Sentinel-2 imagery. For C-O, the presence of a character-
istic band at 2200 nm corresponds to the B12 band of the Sentinel-2
imagery. The B6 at 740 nm from Sentinel-2 imagery might be one of
the representative bands for N-H (Viscarra Rossel et al., 2006).

Meanwhile, Xu et al. (2018) has identified that SOC showed high selectivity ratio value in the visible region of 580–800 nm (i.e., bands 4–7 of Sentinel-2), which is largely influenced by iron oxides and the inherent color of the SOC.

Previous studies have found that combining different spectral bands
or transforming the spectrum effectively enhances model performance
in estimating soil properties. This improvement is attributed to the
enhanced spectral signals and selection of sensitive bands (Liu et al.,
2022b; Wang et al., 2022a). Our results align with previous findings, as
<caption>**Table 6**</caption>

Predictive performance of SOC prediction in the whole China by the global modelling and the difference with the regional modelling*.

<table><thead><tr><th>Covariate</th><th>R<sup>b</sup></th><th>RMSE (g kg<sup>-1</sup>)</th><th>ΔR<sup>b</sup></th><th>ΔRMSE</th></tr></thead><tbody><tr><td>Single Bands</td><td>0.27</td><td>6.69</td><td>0.01</td><td>-0.06</td></tr><tr><td>New Indices</td><td>0.26</td><td>6.73</td><td>0.09</td><td>-0.44</td></tr><tr><td>Existing Indices</td><td>0.24</td><td>6.80</td><td>0.05</td><td>-0.22</td></tr><tr><td>Single Bands with New Indices</td><td>0.26</td><td>6.72</td><td>0.08</td><td>-0.36</td></tr><tr><td>New and Existing Indices</td><td>0.29</td><td>6.59</td><td>0.06</td><td>-0.31</td></tr><tr><td>All remote sensing variables</td><td>0.28</td><td>6.62</td><td>0.04</td><td>-0.20</td></tr><tr><td>All covariates</td><td>0.57</td><td>5.14</td><td>0.01</td><td>-0.06</td></tr><tr><td>All covariates after FRFS</td><td>0.61</td><td>4.88</td><td>0.01</td><td>-0.04</td></tr></tbody></table>

* The ΔR² and ΔRMSE were obtained by subtracting the global modeling results from the regional modeling results.

<caption>**Table 7**</caption>

Comparison with similar studies conducted in various regions worldwide based on accuracy indices and basic information.

<table><thead><tr><th>Source</th><th>Spatial scale</th><th>Spatial resolution</th><th>RMSE (g kg<sup>-1</sup>)</th><th>R<sup>2</sup></th><th>Method</th></tr></thead><tbody><tr><td>Our predictions</td><td>National</td><td>10 m</td><td>4.84</td><td>0.62</td><td>RF</td></tr><tr><td>NSGC</td><td>National</td><td>90 m</td><td>—</td><td>—</td><td>RF</td></tr><tr><td>SoilGrids250m 2.0</td><td>Global</td><td>250 m</td><td>8.41</td><td>0.14</td><td>RF</td></tr><tr><td>HWSD</td><td>Global</td><td>1000 m</td><td>8.50</td><td>0.15</td><td>—</td></tr></tbody></table>

predictors using the indices by bands combination, the model performances were improved with R² increasing by 0.05–0.12 and RMSE decreasing by 0.27–0.59 g kg⁻¹ compared to the single bands (Table 3). The improvements mainly derived from the new constructed indices. From the most important nine bands combination indices, it can be seen that each band was used at least one time. However, modelling by remote sensing indices didn’t achieve the proximal spectrum performance. It is the soil moisture, roughness, crop residues and other field conditions affect a lot to remote sensing directly observation (Ji et al., 2015; Li et al., 2023). Meanwhile, the reflectance characteristics of the soil can be influenced by images processing errors and mixed pixels, especially in the southern China (Zhang et al., 2024).

Although remote sensing techniques have the potential, adding variables with spatial variation information can enhance the model performance in the large areas. Our findings show that the best SOC predictions were obtained by combining remote sensing indices with environmental covariates after FRFS (R² of 0.62, RMSE of 4.84 g kg⁻¹). The climate and position covariates account for the majority, though the FRFS selected different variables in each zone, which was in line with Liang et al. (2019), Luo et al. (2023) mapping soil organic matter overall or part of China. Climatic conditions, which affect both carbon inputs into the soil and SOC decomposition, play a crucial role in determining SOC content globally as well as at regional scales (Wiesmeier et al., 2019). Meanwhile, monsoon climate dominates eastern China (zone 1 and 3), causing the northeast wind from the continental direction in autumn and winter, and the southwesterly wind from the oceanic direction in spring and summer. Therefore, temperature and precipitation show regular changes in this direction, affecting the spatial patterns of SOC (Zhang et al., 2023a). Besides, in the large-scale areas, combining OGCs variables with other environmental data can provide greater accuracy and more realistic spatial patterns (Møller et al., 2020). There are still four remote sensing indices left after variable selection by FRFS, which also verified that the remote sensing contributes a lot in SOC prediction (Fig. 7). It is important to note that feature selection methods such as FRFS and RFE are effective in addressing irrelevant predictors (i.e., those with low predictive power) but do not fully resolve issues of multicollinearity. Although RF are generally insensitive to multicollinearity, caution should be exercised when interpreting variable importance metrics derived from RF in the presence of multicollinearity.

Meanwhile, estimating national-scale SOC according to the agricultural zones can reveal the mechanisms of key covariates, which improves the model performance as well. We quantified the difference of model performance between global and the regional modeling strategy. The results obviously presented that the regional modeling outperformed the global modeling among the different predictor models, with a ΔR² ranging from 0.01 to 0.09 and a ΔRMSE of 0.04 to 0.44 g kg⁻¹ (Table 6). The improvement of regional modeling is due to it can more accurately capture the soil-landscape relationships to the global model (Brungard et al., 2021; Wang et al., 2024). The dominant factors of different zones showed in Fig. 7 were also verified it.

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Surface (0–20 cm) SOC maps excerpted from this study, NSGC, SoilGrids250, and HWSD, respectively, in a 12 km × 6 km window (124.36–124.52°E and 43.25–43.33°N) located in Jilin Province.</figcaption></figure>
## 4.2. Uncertainty analysis

The 90 % prediction intervals were shown in Fig. 8. Several reasons may result in the uncertainty of our SOC estimation. Firstly, the uncertainty is from the high spatial heterogeneity. It showed that YGP and SCB had relative high uncertainty, where the patterns of fragmented and small cropland, and dramatically changeable altitude (Zhang et al., 2023b). On the contrary, the croplands of zone 1 in the east and northeast China are integrated and continuous, where the uncertainty show relatively low (Xue et al., 2024b). Meanwhile, the different use of the cropland may bring uncertainty. The high uncertainty occurred in the primary grain-producing regions extending to the farming-pastoral areas along the eastern edge of the northern and semi-arid zones. Secondly, most high uncertainty areas have sparse samples (Liu et al., 2022a). It can be seen that the uncertainty values in the zone 2 and 3 were higher than the other areas. The different soil sampling schemes may affect the model performance, which could be explored in the future and reduce the uncertainty caused by soil sampling design (Hu et al., 2024a). Thirdly, we used the Sentinel-2 images from 2018 to 2022, but ignored the time of soil sampling. Even though there are tiny changes of SOC over the survey period, considering the soil measurement time will be beneficial, which could reduce the prediction uncertainty (Helfenstein et al., 2022). Finally, the coarse original spatial resolution of the climate covariates (1 km) could bring uncertainty to current fine-resolution mapping. Thus, finer spatial resolution covariates are necessary to DSM works.

## 4.3. Comparison to others SOC maps

The spatial patterns of SOC in our study largely aligns with the estimation made by Liu et al. (2022a) using quantile regression forest and by Ding et al. (2023) based on the DNDC model. We compared our model performance (Table 7) and the SOC spatial pattern in a zoomed-in area of northeast China (Fig. 9) with three widely referenced datasets: SoilGrids250m 2.0 (Poggio et al., 2021), the Harmonized World Soil Database (HWSD, Nachtergaele et al., 2023), and the National Soil Information Grids of China (NSGC, Liu et al., 2022a). As comparing the accuracy metrics calculated using the validation dataset (N = 322), our results ($R^2$ = 0.62, RMSE = 4.84 $g kg^{-1}$) outperformed both Soil-Grids250 ($R^2$ = 0.14, RMSE = 8.41 $g kg^{-1}$) and HWSD ($R^2$ = 0.15, RMSE = 8.50 $g kg^{-1}$). Due to the use of different training and validation sets of the same soil data, we didn't compare the accuracy metrics with NSGC map. The introduction of bare soil spectra information, finer spatial resolution of environmental covariates may enhance the accuracy of our predictions. Meanwhile, more extensive soil survey across China could significantly improve the model performance (Hu et al., 2024b; Liu et al., 2022a).

Zooming into a 12 km × 6 km window (124.36–124.52E and 43.25–43.33 N) located in Jilin Province, we displayed a visual comparison of three surface (0–20 cm) SOC maps excerpted from this study, NSGC (Liu et al., 2022a), SoilGrids250m 2.0 (Poggio et al., 2021), and HWSD (Nachtergaele et al., 2023), respectively. The SOC map produced in this study displayed significant differences in spatial details, revealing spatial variability and local landscape patterns more clearly than other maps with spatial resolutions between 90 and 1000 m. Meanwhile, Fig. 8 showed the SOC contents at the field scale, clearly illustrating the variations in SOC between different fields. Thus, a fine SOC map with a 10 m spatial resolution is crucial for accurately depicting the spatial variability of SOC in Chinese croplands, offering significant benefits for precision agriculture and soil management.

## 4.4. Limitations and perspectives

Despite the progress made in this study, certain limitations remain that will be addressed in our ongoing work. Firstly, although a lot of satellite images are open accessed, they are not fully utilized. Combining

the images from multi-satellites, such as Landsat and Sentinel series, sub-weekly or even daily land surface observations could be provided at medium to high spatial resolutions (Roy et al., 2021; Zhu et al., 2022). Such dense time-series is potentially a treatment to the images' shortage due to the rainy and cloudy weather, especially in southern China, where there is relatively higher planting intensity and shorter bare soil period (Silvero et al., 2021).

Secondly, further studies are needed to develop strategies for fully integrating synthetic aperture radar (SAR; e.g., Sentinel-1) with optical time series images. The independence of SAR sensors from clouds and daytime/nighttime variations enables their capability to capture the structural characteristics of bare soil, crop residues, and other land covers (La et al., 2017). Besides, applying the soil moisture removal algorithms to remote sensing are necessary, which can improve the model performance by remote sensing directly observation.

Finally, soil management (e.g., irrigation measures, fertilizer input, straw return) is a useful information for SOC estimation. However, there is rare researches incorporate variables that characterize soil management in most cases. With the advancement of remote sensing technology, the crop residue cover mapping can be acquired (Quemada et al., 2018; Xiang et al., 2022). Meanwhile, the synthetic fertilizer and manure input data can be found in some regional panel data. Therefore, these factors can be used in SOC prediction in the future, which will improve the model performance and better reveal the SOC variation in croplands. Meanwhile, using deep learning methods, such as Long Short-Term Memory networks (Wang et al., 2022b) and Convolutional Neural Networks (Meng et al., 2023), is an effective way to improve prediction accuracy.

## 5. Conclusions

This study introduced an innovative framework for SOC estimation and generated a high-resolution (10 m) SOC map for Chinese cropland. We proved incorporating remote sensing indices calculated from bare soil composite with environment covariates in RF model after FRFS method could enhance the prediction accuracy of SOC, achieving an $R^2$ of 0.62 and an RMSE of 4.84 $g kg^{-1}$. The model performance is better than the existing SOC maps in such large area. The highest predicted SOC content was in Northeast China. The climate, position, and remote sensing indices were identified the most important covariates (zone 1–3: 85.5, 75.86, 80.29 %) for mapping SOC. Our findings provided an effective way to the high-resolution and high-accuracy SOC monitor at a large scale, which contributes to land management and the development of policies for precision agriculture.

### CRediT authorship contribution statement

**Jie Xue:** Writing – original draft, Software, Methodology, Data curation.
**Xianglin Zhang:** Validation, Data curation.
**Songchao Chen:** Writing – review & editing, Resources, Methodology, Conceptualization.
**Zhongxing Chen:** Investigation, Data curation.
**Rui Lu:** Resources.
**Feng Liu:** Writing – review & editing, Validation, Data curation.
**Bas van Wesemael:** Writing – review & editing, Supervision, Formal analysis.
**Zhou Shi:** Writing – review & editing, Supervision, 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

This study was supported by the National Key Research and Development Program of China (2022YFB3903505, 2022YFB3903503). Jie
Xue also gratefully acknowledges the China Scholarship Council foundation and the Academic Rising Star Program for Doctoral Students of Zhejiang University.

## Appendix A. Supplementary data

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

## Data availability

The authors do not have permission to share data.

## References

Abbas, A., Khan, S., 2007. Using remote sensing techniques for appraisal of irrigated soil salinity. In Proceedings of the International Congress on Modelling and Simulation (MODSIM), Christchurch, New Zealand, 10–13 December; pp. 2632–2638.

Bannari, A., Morin, D., Bonn, F., Huete, A.R., 1995. A review of vegetation indices. Remote Sens. Rev. 13, 95–120.

Boettinger, J.L., Ramsey, R.D., Bodily, J.M., Cole, N.J., Kienast-Brown, S., Nield, S.J., Saunders, A.M., Stum, A.K., 2008. Landsat Spectral Data for Digital Soil Mapping, in: Hartemink, A.E., McBratney, A., Mendonça-Santos, M.d.L. (Eds.), Digital Soil Mapping with Limited Data. Springer Netherlands, Dordrecht, pp. 193–202.

Breiman, L., 2001. Random Forests. Mach. Learn. 45, 5–32.

Brungard, C., Nauman, T., Duniway, M., Veblen, K., Nehring, K., White, D., Salley, S., Anchang, J., 2021. Regional ensemble modeling reduces uncertainty for digital soil mapping. Geoderma 397, 114998.

Castaldi, F., Halil Koparan, M., Wetterlind, J., Zydelis, R., Vinci, I., Özge Savas, A., Kivrak, C., Tunçay, T., Volungevićius, J., Obber, S., Ragazzi, F., Malo, D., Vaudour, E., 2023. Assessing the capability of Sentinel-2 time-series to estimate soil organic carbon and clay content at local scale in croplands. ISPRS J. Photogramm. Remote Sens. 199, 40–60.

Chen, J., Chen, J., Liao, A., Cao, X., Chen, L., Chen, X., He, C., Han, G., Peng, S., Lu, M., Zhang, W., Tong, X., Mills, J., 2015. Global land cover mapping at 30m resolution: A POK-based operational approach. ISPRS J. Photogramm. Remote Sens. 103, 7–27.

Chen, S., Arrouays, D., Leatitia Mulder, V., Poggio, L., Minasny, B., Roudier, P., Libohova, Z., Lagacherie, P., Shi, Z., Hannam, J., Meersmans, J., Richer-de-Forges, A.C., Walter, C., 2022. Digital mapping of GlobalSoilMap soil properties at a broad scale: A review. Geoderma 409, 115567.

Chen, S., Chen, Z., Zhang, X., Luo, Z., Schillaci, C., Arrouays, D., Richer-de-Forges, A.C., Shi, Z., 2024. European topsoil bulk density and organic carbon stock database (0–20 cm) using machine-learning-based pedotransfer functions. Earth Syst. Sci. Data 16 (5), 2367–2383.

Cooperative Research Group on Chinese Soil Taxonomy, 2001. Keys to Chinese Soil Taxonomy, 3rd ed. Press of University of Science and Technology of China, Hefei.

Dematte, J.A.M., Fongaro, C.T., Rizzo, R., Safanelli, J.L., 2018. Geospatial Soil Sensing System (GEOS): A powerful data mining procedure to retrieve soil spectral reflectance from satellite images. Remote Sens. Environ. 212, 161–175.

Diek, S., Fornallaz, F., Schaepman, M.E., De Jong, R., 2017. Barest Pixel Composite for Agricultural Areas Using Landsat Time Series. Remote Sens. (Basel) 9, 1245.

Ding, W., Chang, N., Zhang, G., Kang, J., Yi, X., Zhang, J., Zhang, J., Wang, L., Li, H., 2023. Soil organic carbon changes in China's croplands: A newly estimation based on DNDC model. Sci. Total Environ. 905, 167107.

Dvorakova, K., Heiden, U., Pepers, K., Staats, G., van Os, G., van Wesemael, B., 2023. Improving soil organic carbon predictions from a Sentinel-2 soil composite by assessing surface conditions and uncertainties. Geoderma 429, 116128.

Fick, S.E., Hijmans, R.J., 2017. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315.

Ghaemi, M., Astaraei, A.R., Sanaeinejad, S.H., Zare, H., 2013. Using satellite data for soil condition exchange capacity studies. Int. Agrophys. 27, 409–417.

Guyon, I., Weston, J., Barnhill, S., Vapnik, V., 2002. Gene Selection for Cancer Classification using Support Vector Machines. Mach. Learn. 46, 389–422.

Helfenstein, A., Mulder, V.L., Heuvelink, G.B.M., Okx, J.P., 2022. Tier 4 maps of soil pH at 25 m resolution for the Netherlands. Geoderma 410, 115659.

Hewson, R.D., Cudahy, T.J., Huntington, J.F., 2001. Geologic and alteration mapping at Mt Fitton, South Australia, using ASTER satellite-borne data, IGARSS 2001. Scanning the Present and Resolving the Future. Proceedings. IEEE 2001 International Geoscience and Remote Sensing Symposium (Cat. No.01CH37217), pp. 724-726 vol.722.

Hu, B., Xie, M., Shi, Z., Li, H., Chen, S., Wang, Z., Zhou, Y., Ni, H., Geng, Y., Zhu, Q., Zhang, X., 2024a. Fine-resolution mapping of cropland topsoil pH of Southern China and its environmental application. Geoderma 442.

Hu, B., Xie, M., Zhou, Y., Chen, S., Zhou, Y., Ni, H., Peng, J., Ji, W., Hong, Y., Li, H., Shi, Z., 2024b. A high-resolution map of soil organic carbon in cropland of Southern China. Catena 237, 107813.

Huete, A., Didan, K., Miura, T., Rodriguez, E.P., Gao, X., Ferreira, L.G., 2002. Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 83, 195–213.

Ji, W., Viscarra Rossel, R.A., Shi, Z., 2015. Accounting for the effects of water and the environment on proximally sensed vis-NIR soil spectra and their calibrations. Eur. J. Soil Sci. 66, 555–565.

Kuhn, M., 2008. Building Predictive Models in R Using the caret Package. J. Stat. Softw. 28, 1–26.

La, T.V., Khenchaf, A., Comblet, F., Nahum, C., 2017. Exploitation of C-Band Sentinel-1 Images for High-Resolution Wind Field Retrieval in Coastal Zones (Iroise Coast, France). IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 10, 5458–5471.

Li, W., Li, D., Liu, S., Baret, F., Ma, Z., He, C., Warner, T.A., Guo, C., Cheng, T., Zhu, Y., Cao, W., Yao, X., 2023. RSARE: A physically-based vegetation index for estimating wheat green LAI to mitigate the impact of leaf chlorophyll content and residue-soil background. ISPRS J. Photogramm. Remote Sens. 200, 138–152.

Li, Y., 2010. Can the Spatial Prediction of Soil Organic Carbon Contents at Various Sampling Scales Be Improved by Using Regression Kriging with Auxiliary Information? Geoderma 159, 63–75.

Liang, Z., Chen, S., Yang, Y., Zhao, R., Shi, Z., Viscarra Rossel, R.A., 2019. National digital soil map of organic matter in topsoil and its associated uncertainty in 1980's China. Geoderma 335, 47–56.

Liu, F., Wu, H., Zhao, Y., Li, D., Yang, J.L., Song, X., Shi, Z., Zhu, A.X., Zhang, G.L., 2022a. Mapping high resolution National Soil Information Grids of China. Sci. Bullet. (Beijing) 67, 328–340.

Liu, Q., He, L., Guo, L., Wang, M., Deng, D., Lv, P., Wang, R., Jia, Z., Hu, Z., Wu, G., Shi, T., 2022b. Digital mapping of soil organic carbon density using newly developed bare soil spectral indices and deep neural network. Catena 219.

Luo, C., Zhang, W., Zhang, X., Liu, H., 2023. Mapping soil organic matter content using Sentinel-2 synthetic images at different time intervals in Northeast China. Int. J. Digital Earth 16, 1094–1107.

Malone, B.P., Jha, S.K., Minasny, B., McBratney, A.B., 2016. Comparing regression-based digital soil mapping and multiple-point geostatistics for the spatial extrapolation of soil data. Geoderma 262, 243–253.

Mandal, U.K., 2016. Spectral color indices based geospatial modeling of soil organic matter in Chitwan District, Nepal. Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci. XLI-B2, 43–48.

McBratney, A.B., Mendonça Santos, M.L., Minasny, B., 2003. On digital soil mapping. Geoderma 117, 3–52.

Meng, X.T., Bao, Y.L., Luo, C., Zhang, X.L., Liu, H.J., 2023. SOC content of global Mollisols at a 30 m spatial resolution from 1984 to 2021 generated by the novel ML-CNN prediction model. Remote Sens. Environ. 300, 113911.

Mladenova, I.E., Bolten, J.D., Crow, W., Sazib, N., Reynolds, C., 2020. Agricultural drought monitoring via the assimilation of SMAP soil moisture retrievals into a global soil water balance model. Front. Big Data 3, 10. https://doi.org/10.3389/fdata.2020.00010.

Møller, A.B., Beucher, A.M., Pouladi, N., Greve, M.H., 2010. Oblique geographic coordinates as covariates for digital soil mapping. SOIL 6, 269–289.

Nachtergaele, F., van Velthuizen, H., Verelst, L., Wiberg, D., Henry, M., Chiozza, F., et al., 2023. Harmonized World Soil Database version 2.0. Food and Agriculture Organization of the United Nations.

Nguyen, C.T., Chidthaisong, A., Kieu Diem, P., Huo, L.Z., 2021. A modified bare soil index to identify bare land features during agricultural fallow-period in Southeast Asia using Landsat 8. Land 10, 231.

Poggio, L., de Sousa, L.M., Batjes, N.H., Heuvelink, G.B.M., Kempen, B., Ribeiro, E., Rossier, D., 2021. SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. SOIL 7, 217–240.

Pouladi, N., Møller, A., Tabatabai, S., Greve, M., 2019. Mapping Soil Organic Matter Contents at Field Level with Cubist, Random Forest and Kriging. Geoderma 342, 85–92.

Quemada, M., Hively, W.D., Daughtry, C.S.T., Lamb, B.T., Shermeyer, J., 2018. Improved crop residue cover estimates obtained by coupling spectral indices for residue and moisture. Remote Sens. Environ. 206, 33–44.

Rogge, D., Bauer, A., Zeidler, J., Mueller, A., Esch, T., Heiden, U., 2018. Building an exposed soil composite processor (SCMAP) for mapping spatial and temporal characteristics of soils with Landsat imagery (1984–2014). Remote Sens. Environ. 205, 1–17.

Roy, D.P., Huang, H., Houborg, R., Martins, V.S., 2021. A global analysis of the temporal availability of PlanetScope high spatial resolution multi-spectral imagery. Remote Sens. Environ. 264, 112586.

Running, S., Mu, Q., Zhao, M., 2015. MOD17A2H MODIS/Terra Gross Primary Productivity 8-Day L4 Global 500m SIN Grid V006. In: N.E.L.P. DAAC (Ed.).

Safanelli, J.L., Poppiel, R.R., Ruiz, L.F., Bonfatti, B.R., Mello, F.A., Rizzo, R., Dematte, J. A.M., 2020. Terrain Analysis in Google Earth Engine: A Method Adapted for High-Performance Global-Scale Analysis. ISPRS Int. J. Geo-Inform.

Shi, T., Collins, S.L., Yu, K., Peñuelas, J., Sardans, J., Li, H., Ye, J., 2024. A global meta-analysis on the effects of organic and inorganic fertilization on grasslands and croplands. Nat. Commun. 15, 3411.

Shi, Z., Wang, Q., Peng, J., Ji, W., Liu, H., Li, X., Viscarra Rossel, R.A., 2014. Development of a national VNIR soil-spectral library for soil classification and prediction of organic matter concentrations. Sci. China Earth Sci. 57, 1671–1680.

Silvero, N.E.Q., Dematte, J.A.M., Amorim, M.T.A., Santos, N.V.d., Rizzo, R., Safanelli, J. L., Poppiel, R.R., Mendes, W.d.S., Bonfatti, B.R., 2021. Soil variability and quantification based on Sentinel-2 and Landsat-8 bare soil images: a comparison. Remote Sens. Environ. 252, 112117.

Simbahan, G., Dobermann, A., Goovaerts, P., Ping, J., Haddix, M., 2006. Fine-Resolution Mapping of Soil Organic Carbon Based on Multivariate Secondary Data. Geoderma 132, 471–489.

Teng, H., Hu, J., Zhou, Y., Zhou, L., Shi, Z., 2019. Modelling and mapping soil erosion potential in China. J. Integr. Agric. 18, 251–264.

Van Oost, K., Quine, T.A., Govers, G., De Gryze, S., Six, J., Harden, J.W., Ritchie, J.C., McCarty, G.W., Heckrath, G., Kosmas, C., Giraldez, J.V., da Silva, J.R.M., Merckx, R.,
2007. The Impact of Agricultural Soil Erosion on the Global Carbon Cycle. Science 318, 626–629.

Vaudour, E., Gholizadeh, A., Castaldi, F., Saberioon, M., Borùvka, L., Urbina-Salazar, D., Fouad, Y., Arrouays, D., Richer-de Forges, A.C., Biney, J., et al., 2022. Satellite imagery to map topsoil organic carbon content over cultivated areas: an overview. Remote Sens 14 (12), 2917.

Vaudour, E., Gomez, C., Fouad, Y., Lagacherie, P., 2019. Sentinel-2 image capacities to predict common topsoil properties of temperate and Mediterranean agroecosystems. Remote Sens. Environ. 223, 21–33.

Vescovo, L., Gianelle, D., 2008. Using the MIR bands in vegetation indices for the estimation of grassland biophysical parameters from satellite remote sensing in the Alps region of Trentino (Italy). Adv. Space Res. 41, 1764–1772.

Viscarra Rossel, R.A., Behrens, T., Ben-Dor, E., Brown, D.J., Demattè, J.A.M., Shepherd, K.D., Shi, Z., Stenberg, B., Stevens, A., Adamchuk, V., Aïchi, H., Barthès, B.G., Bartholomeus, H.M., Bayer, A.D., Bernoux, M., Böttcher, K., Brodský, L., Du, C.W., Chappell, A., Fouad, Y., Genot, V., Gomez, C., Grunwald, S., Gubler, A., Guerrero, C., Hedley, C.B., Knadel, M., Morrás, H.J.M., Nocita, M., Ramirez-Lopez, L., Roudier, P., Campos, E.M.R., Sanborn, P., Sellitto, V.M., Sudduth, K.A., Rawlins, B.G., Walter, C., Winowiecki, L.A., Hong, S.Y., Ji, W., 2016. A global spectral library to characterize the world's soil. Earth Sci. Rev. 155, 198–230.

Viscarra Rossel, R.A., Walvoort, D.J.J., McBratney, A.B., Janik, L.J., Skjemstad, J.O., 2006. Visible, near infrared, mid infrared or combined diffuse reflectance spectroscopy for simultaneous assessment of various soil properties. Geoderma 131, 59–75.

Wang, N., Chen, S., Huang, J., Frappart, F., Taghizadeh, R., Zhang, X., Wigneron, J.-P., Xue, J., Xiao, Y., Peng, J., Shi, Z., 2024. Global Soil Salinity Estimation at 10 m Using Multi-source Remote Sensing. J. Remote Sens. 4, 0130.

Wang, N., Peng, J., Chen, S., Huang, J., Li, H., Biswas, A., He, Y., Shi, Z., 2022a. Improving remote sensing of salinity on topsoil with crop residues using novel indices of optical and microwave bands. Geoderma 422.

Wang, S., Guan, K.Y., Zhang, C.H., Lee, D.K., Margenot, A.J., Ge, Y.F., Peng, J., Zhou, W., Zhou, Q., Huang, Y.Z., 2022b. Using soil library hyperspectral reflectance and machine learning to predict soil organic carbon: Assessing potential of airborne and spaceborne optical soil sensing. Remote Sens. Environ. 271, 112914.

Wei, W., Pang, S., Wang, X., Zhou, L., Xie, B., Zhou, J., Li, C., 2020. Temperature vegetation precipitation dryness index (TVPDI)-based dryness-wetness monitoring in China. Remote Sens. Environ. 248, 111957.

Wiesmeier, M., Urbanski, L., Hobley, E., Lang, B., von Lützow, M., Marin-Spiotta, E., van Wesemael, B., Rabot, E., Ließ, M., Garcia-Franco, N., Wollschläger, U., Vogel, H.-J., Kögel-Knabner, I., 2019. Soil organic carbon storage as a key function of soils - A review of drivers and indicators at various scales. Geoderma 333, 149–162.

Xiang, X., Du, J., Jacinthe, P.-A., Zhao, B., Zhou, H., Liu, H., Song, K., 2022. Integration of tillage indices and textural features of Sentinel-2A multispectral images for maize residue cover estimation. Soil Tillage Res. 221, 105405.

Xiao, Y., Xue, J., Zhang, X., Wang, N., Hong, Y., Jiang, Y., Zhou, Y., Teng, H., Hu, B., Lugato, E., Richer-de-Forges, A.C., Arrouays, D., Shi, Z., Chen, S., 2022. Improving pedotransfer functions for predicting soil mineral associated organic carbon by ensemble machine learning. Geoderma 428, 116208.

Xu, D., Ma, W., Chen, S., Jiang, Q., He, K., Shi, Z., 2018. Assessment of important soil properties related to Chinese Soil Taxonomy based on vis-NIR reflectance spectroscopy. Comput. Electron. Agric. 144, 1–8.

Xu, D., Zhao, R., Li, S., Chen, S., Jiang, Q., Zhou, L., Shi, Z., 2019. Multi-sensor fusion for the determination of several soil properties in the Yangtze River Delta, China. Eur. J. Soil Sci. 70, 162–173.

Xue, J., Zhang, X., Chen, S., Hu, B., Wang, N., Shi, Z., 2024a. Quantifying the agreement and accuracy characteristics of four satellite-based LULC products for cropland classification in China. J. Integr. Agric. 23, 283–297.

Xue, J., Zhang, X., Chen, S., Lu, R., Wang, Z., Wang, N., Hong, Y., Chen, X., Xiao, Y., Ma, Y., Shi, Z., 2023. The validity domain of sensor fusion in sensing soil quality indicators. Geoderma 438, 116657.

Xue, J., Zhang, X., Huang, Y., Chen, S., Dai, L., Chen, X., Yu, Q., Ye, Su., Shi, Z., 2024b. A two-dimensional bare soil separation framework using multi-temporal Sentinel-2 images across China. International Journal of Applied Earth Observation and Geoinformation 134, 104181.

Zhang, G., Gong, Z., 2012. Soil survey laboratory methods. Science Press, Beijing, pp. 8–23.

Zhang, X., Chen, S., Xue, J., Wang, N., Xiao, Y., Chen, Q., Hong, Y., Zhou, Y., Teng, H., Hu, B., Zhuo, Z., Ji, W., Huang, Y., Gou, Y., Richer-de-Forges, A.C., Arrouays, D., Shi, Z., 2023a. Improving model parsimony and accuracy by modified greedy feature selection in digital soil mapping. Geoderma 432, 116383.

Zhang, X., Xue, J., Chen, S., Wang, N., Xie, T., Xiao, Y., Chen, X., Shi, Z., Huang, Y., Zhuo, Z., 2023b. Fine Resolution Mapping of Soil Organic Carbon in Croplands with Feature Selection and Machine Learning in Northeast Plain China. Remote Sens. (Basel) 15, 5033.

Zhang, X., Xue, J., Chen, S., Wang, N., Shi, Z., Huang, Y., Zhuo, Z., 2022. Digital Mapping of Soil Organic Carbon with Machine Learning in Dryland of Northeast and North Plain China. Remote Sens. (Basel) 14, 2504.

Zhang, X., Xue, J., Chen, S., Zhuo, Z., Wang, Z., Chen, X., Xiao, Y., Shi, Z., 2024. Improving model performance in mapping cropland soil organic matter using time-series remote sensing data. J. Integr. Agric. 23, 2820–2841.

Zhang, X., Xue, J., Xiao, Y., Shi, Z., Chen, S., 2023c. Towards Optimal Variable Selection Methods for Soil Property Prediction Using a Regional Soil Vis-NIR Spectral Library. Remote Sens. (Basel) 15, 465.

Zhou, Y., Chartin, C., Van Oost, K., van Wesemael, B., 2022. High-resolution soil organic carbon mapping at the field scale in Southern Belgium (Wallonia). Geoderma 422.

Zhou, Y., Chen, S., Zhu, A.X., Hu, B., Shi, Z., Li, Y., 2021. Revealing the scale- and location-specific controlling factors of soil organic carbon in Tibet. Geoderma 382, 114713.

Zhu, Z., Qiu, S., Ye, S., 2022. Remote sensing of land change: A multifaceted perspective. Remote Sens. Environ. 282, 113266.

Zuo, L., Zhang, Z., Carlson, K.M., MacDonald, G.K., Brauman, K.A., Liu, Y., Zhang, W., Zhang, H., Wu, W., Zhao, X., Wang, X., Liu, B., Yi, L., Wen, Q., Liu, F., Xu, J., Hu, S., Sun, F., Gerber, J.S., West, P.C., 2018. Progress towards sustainable intensification in China challenged by land-use change. Nat. Sustainability 1, 304–313.

Zou, Z., Dong, J., Menarguez, M.A., Xiao, X., Qin, Y., Doughty, R.B., Hooker, K.V., David Hambright, K., 2017. Continued decrease of open surface water body area in Oklahoma during 1984–2015. Sci. Total Environ 595, 451–460.