# Digital mapping of the soil thickness of loess deposits over a calcareous bedrock in central France

Songchao Chen, Anne C Richer-De-Forges, Vera Leatitia Mulder, Guillaume

Martelet, Thomas Loiseau, Sébastien Lehmann, Dominique Arrouays

► To cite this version:

Songchao Chen, Anne C Richer-De-Forges, Vera Leatitia Mulder, Guillaume Martelet, Thomas Loiseau, et al.. Digital mapping of the soil thickness of loess deposits over a calcareous bedrock in central France. CATENA, 2021, 198, pp.105062. 10.1016/j.catena.2020.105062. hal-03467117

HAL Id: hal-03467117

https://hal.inrae.fr/hal-03467117

Submitted on 3 Feb 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
**Title:** Digital mapping of the soil thickness of loess deposits over a calcareous bedrock in central France

**Authors:**

Songchao Chenª. songchao.chen@inrae.fr

Anne C. Richer-de-Forgesª. anne.richer-de-forges@inrae.fr

Vera Leatitia Mulderᵇ. titia.mulder@wur.nl

Guillaume Marteletᶜ. g.martelet@brgm.fr

Thomas Loiseauª. thomas.loiseau@inrae.fr

Sébastien Lehmannª. sebastien.lehmann@inrae.fr

Dominique Arrouaysª. dominique.arrouays@inrae.fr

**Affiliations:**

ª INRAE, Unité InfoSol, 45075 Orléans, France

ᵇ Soil Geography and Landscape Group, Wageningen University, PO Box 47 6700

AA Wageningen, The Netherlands

ᶜ BRGM, UMR 7327, 45060 Orléans, France

**Corresponding author:**

Songchao Chen: songchao.chen@inrae.fr

Postal address: INRAE, Unité InfoSol, 2163 Avenue de la Pomme de Pin, CS 40001

Ardon, 45075 Orléans, France
1 Abstract:

2 Soil thickness (ST) plays an important role in regulating soil processes,
3 vegetation growth and land suitability. Therefore, it has been listed as one of twelve
4 basic soil properties to be delivered in *GlobalSoilMap* project. However, ST prediction
5 has been reported with poor performance in previous studies. Our case study is
6 located in the intensive agriculture Beauce area, central France. In this region, the ST
7 mainly depends on the thickness of loess (TOL) deposits over a calcareous bedrock.

8 We attempted to test the TOL prediction by coupling a large soil dataset (10978
9 sampling sites) and 117 environmental covariates. After variable selection by
10 recursive feature elimination, quantile regression forests (QRF) was employed for
11 spatial modelling, as it was able to directly provide the 90% prediction intervals (PIs).
12 Averaging a total of 50 models, generated by repeated stratified random sampling,
13 showed a substantial model performance with mean R² of 0.33, RMSE of 30.48 cm
14 and bias of -1.20 cm. The prediction interval coverage percentage showed that 86.70%
15 of the validation samples fall within the predefined 90% PIs, which also indicated the
16 prediction uncertainty produced by QRF was reasonable. The relative variable
17 importance indicated the importance of airborne gamma-ray radiometric data and
18 Sentinel 2 products in TOL prediction. The produced TOL map with 90% PIs makes
19 sense from a soil science and physiographic point of view. The final product can
20 guide evidence-based decision making for agricultural land management, especially
21 for irrigation in our case study.
# 1. Introduction

Soil thickness (ST) is an important soil property due to its influence as a controlling factor of numerous surface and subsurface soil processes. Through its influence on the plant rootable depth (Leenaars et al., 2018), ST is a major controlling factor of vegetation growth and land suitability, and is a key variable to estimate available water capacity (AWC). As a consequence, ST has been retained as a mandatory soil attribute to be delivered in GlobalSoilMap products (Arrouays et al., 2014). Previously, some attempts have been made to map ST at national, continental or global levels (e.g., Grundy et al., 2015; Lacoste et al., 2016; Mulder et al., 2016; Hengl et al., 2017; Chen et al., 2019). Grundy et al. (2015) mapped ST in Australia using about 300,000 points of observations and environmental covariates as inputs for a Cubist model. Lacoste et al. (2016) tested three digital soil mapping (DSM) approaches, based on regression tree modelling, gradient boosting modelling, and multi-resolution kriging for a dataset of ca 14,000 observations in France. Hengl et al. (2017) mapped ST at a global scale. Predictions were based on ca. 150,000 soil profiles used for training and a stack of 158 covariates which were used to fit an ensemble of machine learning methods—random forest and gradient boosting and/or multinomial logistic regression. Mulder et al. (2016) used Cubist predictions for mainland France using ca 16,000 observation points and a set of 20 spatially exhaustive covariates. Chen et al. (2019) further demonstrated how right-censored data can be accounted for in the ST modelling of mainland France. Using random survival forest, the probability of exceeding a given depth was modelled using freely available spatial data representing the main soil-forming factors. However, most of these results gave rather poor prediction performances compared to other soil properties such as soil organic carbon or clay content and pH (e.g., Mulder et al., 2016; Hengl et al., 2017). In many cases, ST prediction proved to be hampered either 1) by the lack of data measurements (Leenaars et al., 2018) or 2) by the fact that the collected ST data is often right-censored data (i.e. the observed ST is less than true ST, Chen et al., 2019), or 3) because ST has a high short-range spatial variability in specific pedological contexts (e.g., Bourennane et al., 1996; Lacoste et al., 2016). Moreover, most of the examples taken from the literature were produced using digital soil mapping (DSM, McBratney et al., 2003) and in some cases, one may expect that no relevant covariate was available to improve the performance of the predictions.
For example, while topography-related covariates such as elevation or slope often explain large part of the variability of soils in a given areas, these covariates might not improve the performance of predictions in flat areas.

Mapping ST enables several applications among which are agronomy and agricultural practices (plant rootable depth, drainage, irrigation), crop growth modelling, geotechnical engineering, water balance modelling at catchment to global scale, and Quaternary science studies.

The Beauce area, located in central France is a limestone plateau irregularly covered with Quaternary loessic silt (Macaire, 1971; Lorain, 1973; Ménillet, 1974). Soil classification varies from Luvisols to Calcic and Calcaric Cambisols. The region is rather flat, and has a noticeable proportion of rather thin soils. Intensive agriculture in this region often utilizes irrigation and most of the surface area is occupied by cereal crops (mainly maize and wheat) and sugar beet. Also, the Beauce area is home of the largest aquifer of France in the underlying calcareous rock. Upper horizons were affected by peri-glacial winds that redistributed loess deposits (Macaire, 1971; Bourennane et al., 1996), resulting in a rather homogeneous particle-size distribution of the fine earth (i.e., silt, clay, and sand). Therefore, most of the soils consist now of silt, silt loam or silty clay loam layers derived from this aeolian deposit developed on a lacustrine limestone substrate. In general, the illuviation process occurred when the thickness of loess was the largest, resulting in less clayey topsoil textures.

The available water capacity (AWC) is the maximum amount of available water that can be stored for crop growth, therefore it is an important soil information for agricultural management. Therefore, the thickness of loess (TOL) deposit is also one of the primary factors influencing the calculation of the soil AWC. Although Tetegan et al. (2015) demonstrated that the percentage of rock fragments was also one of the controlling factors of AWC in this region. Overall, in this region, irrigation management is of upmost importance in order to maintain crop yields, while preserving the underlying water table and water quality. Knowing the TOL is essential for determining a water balance and for piloting irrigation management. In terms of agronomy and environment, the TOL is a determining factor (Nicoullaud et al., 1995; Ould-Mohamed et al., 1997). Therefore, the TOL should be known accurately and cheaply over the study area. Several traditional soil maps have been produced in this
region, with scales ranging from 1:50,000 to 1:250,000, resulting in a various density of point scale soil information.

The objective of this study is to assess to which extent using this legacy data and environmental covariates (from existing geological maps, digital elevation model derivatives, airborne gamma-ray radiometry, and remote sensing data) in a DSM model allows to reach acceptable performances for TOL prediction. In this study, we decided to model the TOL up to a depth of 120 cm using Quantile Regression Forests (QRF) because the TOL was deemed useful for agricultural practices. Maize cropping is especially of interest because it is known for its high water requirement (Doorenbos et al., 1978) and thus typically requires the largest amounts of irrigation. The average rooting depth of maize is equal to 120 cm (British Standards Institution, 1988; Tetegan et al., 2015). Therefore, we only mapped the TOL up to a depth of 120 cm, as there is no difference of soil water management between soils with a TOL deeper than 120 cm and soils with a TOL of 120 cm.

102

# 2. Material and methods

## 2.1. Study area

This study was conducted in the Beauce area located at the middle Loire catchment, central France (Figure 1). It covers a total area of 4835 km², of which agriculture is the dominant land use (88.5%, Inglada et al., 2017). It has a continental-oceanic climate with a mean annual temperature of 11.5°C and a mean annual rainfall of 700 mm (Paroissien et al., 2014). Most of the soils in this study area are developed from periglacial loess deposits which covered a limestone bedrock. Cambisol (48.3%) and Luvisol (25.6%) are the major soil groups observed in this region (IUSS Working Group WRB, 2006). At the southern border of the Beauce region, some other soil groups (not developed from loess) are observed.

## 2.2. Soil data

We used available soil data from the French Soil Inventory Program (IGCS). The thickness of loess derived horizons (TOL) was determined by several criteria: 1) digging soil pits down to the calcareous material and 2) by auger borings. The presence of small rock fragments could in some case lead to an underestimation of TOL done by augering. Therefore, if a TOL of 120 cm was not reached, two other
augerings were made randomly 0.5 m apart from the first one and the maximum TOL
reached was recorded. TOL should have a texture of silt, silt loam or silty clay loam
(Bertran et al., 2016; Borderie et al., 2017). The deeper the TOL the more the
illuviation processes are pronounced and the lighter the topsoil texture.

In total, 10978 sites were used in this study to map the TOL up to a depth of 120 cm. The TOL for sites with a TOL deeper than 120 cm (n=14) was set to 120 cm before modelling to eliminate the effect of extreme values in modelling.

2.3. Environmental covariates

The environmental covariates used in this study and their data sources are listed
in Table 1. These covariates provide information on the environmental factors
assumingly controlling TOL, based on the Scorpan conceptual model (McBratney et
al., 2003). For illustrative purposes, several covariates are shown in Figure 2.

12.3.1. Relief

The Digital Elevation Model of mainland France was derived from BD TOPO 3 of
the French National Geographical Institute (IGN, 2011), at 25 m resolution. SAGA
GIS (Conrad *et al.*, 2015) was used to calculate its derivatives (relief factors),
including channel network base level (CNBL), multiresolution index of valley bottom
flatness (MrVBF), plan curvature (PIC), profile curvature (PrC), slope (SI), slope
position (SIP), slope length (SIL), terrain wetness index (TWI), valley depth (VD), and
vertical distance to channel network (VDCN). As the relief factor at neighbouring
locations is able to provide additional useful information in modelling soil patterns
(McBratney *et al.*, 2003), some previous studies investigated the potential of
incorporating local neighbourhood information into the training pixels, using
convolution filtering operations (e.g., Grinand *et al.*, 2008; Loiseau *et al.*, 2019).
Filtering can be achieved by passing a moving window over the variable to calculate
a value of the processing cell (central pixel) using the values of its neighbouring cells.
In this study, we used mean convolution circular windows to calculate the focal
means for these relief factors with radius at 200, 500 and 1000 m (Grinand *et al.*,
2008), which resulted in three raster layers derived from each original relief factor (25
m).

2.3.2. Soil
151

The soil type information were extracted from the French national soil type map at 1:1 M scale (King et al., 1994). The soil types in this study area were mainly Cambisols and Luvisols. However, some Podzols, Gleysols, Fluvisols, Arenosols and Vertisols were rarely present, mainly at the southern border of the region.

152

153

154

155

### 2.3.3. Parent material

156

The map of parent material was extracted from the French national parent material map (King et al., 1994). Undifferentiated alluvial deposits, calcareous rocks, clayey materials, sandy materials and loamy materials are the main parent materials in the study area. Note that the loamy materials are nearly always located over underlying calcareous rocks.

157

158

159

160

161

The gamma radiometric data, including Potassium (K), Thorium (Th) and Uranium (U), and total count (TC), was derived from an airborne high-resolution magnetic and radiometric survey over the Région Centre, flown by Terraquest Ltd, Canada, under the supervision of BRGM between 2008 and 2009 (Martelet et al., 2014). The line-spacing of the survey was 1 km and, along the flight lines the footprint of each gamma radiometric measurement was an ellipse of 150 × 250 m²; accordingly the data were interpolated on 250 m grids using a standard minimum curvature interpolation.

162

163

164

165

166

167

168

169

### 2.3.4. Organisms

170

A land use map was extracted from the French land use map, which was produced from Sentinel 2 data at 10 m resolution, for year 2016 (Inglada et al., 2017). This land use map was aggregated to 25 m resolution by majority sampling and the proportions (0~100%) of the nine main land-use classes within each 25×25 m pixel (which contained 6 10×10 pixels) were also included as covariates.

171

172

173

174

175

The monthly normalized difference vegetation index (NDVI) from the MODIS (MCD43A4 16-day Version 6) in 500 m resolution and the PROBA-V 10-day product level 2B TOC (Copernicus, 2016) in 300 m resolution were used in this study. These 24 monthly NDVI data in 2003 (extreme warm and dry year) and 2016 (normal year) were collected and reduced into the first three principal components by principal component analysis to eliminate their multicollinearity. For more details, we refer to Loiseau et al. (2019).

176

177

178

179

180

181
We also included 42 covariates related to Sentinel 2 bands (year of 2016 to 2017) and indices, which were produced in an earlier study from Loiseau et al. (2019) for mainland France at 90 m resolution. The Sentinel 2 data were processed to Level-2A (atmospheric and topographic corrections) by the French National Centre for Space Studies (Hagolle et al., 2015). These covariates included 10 Sentinel 2 bands (2, 3, 4, 5, 6, 7, 8, 8A, 11 and 12), 11 spectral indices (brightness index, saturation index, hue index, coloration index, redness Index, carbonate index, ferrous iron, clay index, normalized difference 1, normalized difference 2 and grain size index) and their focal means determined by a low-pass filter with an average within a 2x2 km window. For more details, we refer to Loiseau et al. (2019).

### 2.3.5. Position

The coordinates, i.e., latitude and longitude (extracted for each recorded sampling site), were used in modelling. In addition, 10 oblique geographic coordinates were calculated at angles of 15°, 30°, 45°, 60°, 75°, 105°, 120°, 135°, 150° and 165°. The oblique coordinate (OC) at an angle of θ can be calculated as below (Møller et al., 2019):

$$OC = \sqrt{X^2 + Y^2} \times \cos(\theta - \tan^{-1}(\frac{Y}{X})) \quad (1)$$

where X and Y are the latitude and longitude.

Note that when θ is 0° or 90°, the oblique coordinate equals to latitude or longitude.

### 2.3.6. Harmonization of environmental covariates

The environmental covariates had different resolutions and scales, we therefore harmonized them at 25 m resolution using nearest neighbour interpolation for spatial predictive modelling and mapping at non-visited locations.

## 2.4. Variable selection using recursive feature elimination

Considering the large set of environmental covariates (n=117), variable selection was applied by recursive feature elimination (Kunn, 2020) prior to fitting the spatial predictive model. The recursive feature elimination (incorporating resampling) adopts a backwards selection, which includes several steps: (1) split data into training and test set by resampling (i.e., k-fold cross-validation); (2) train the model on the training set using all predictors, calculate the model performance on the test set, and rank
predictors using their model importance; (3) for each predictor subset size $S_i$ ($i=1, 2, ..., s$), train the model on the training set using the $S_i$ most important predictors, and calculate the model performance on the test set; (4) compare the model performance profile over the $S_i$ on the test set, and determine the optimal number of predictors.

To select the important covariates and improve the mapping efficiency, the recursive feature elimination was performed on the whole data using rfe function in caret package (Kunn, 2020) in R (R Core Team, 2019). The model was set to Random Forest (default values with tree number of 500 and mtry of $p/3$ where $p$ is the size of predictors) using 5-fold cross-validation. Seven predictor subset sizes (5, 10, 15, 20, 40, 60, 80 and 100) were tested and the model performance indicated that 80 variables (Table 2) were optimal and then used for later modelling.

## 2.5. Spatial predictive modelling and model performance evaluation

Quantile Regression Forest (QRF, Meinshausen, 2006) has been growingly used in DSM for delivering soil information as it is able to provide uncertainty estimates straightforwardly with a fair model performance (e.g., Vaysse and Lagacherie, 2017; Lombardo et al., 2018; Loiseau et al., 2019). Therefore, QRF was used for modelling TOL in this study.

Since QRF is an extension of Random Forest (RF, Breiman, 2001), we start with RF. Assume $X$ and $Y$ are the predictor variables and responses, for regression, RF generates a large number ($b$) of bootstrap trees by using $m$ training samples ($X_i, Y_i$), $i=1,..., m$. Here, bootstrap refers to repeated ($b$ times) selection of a random sample with replacement of the training samples. For each node in a bootstrap tree, a random subset of the predictor variables is used for split-point selection. The prediction of a bootstrap tree for a new sample $D=X_d$ is the conditional mean estimate ($\hat{X}$) of $Y$, which can be represented by:

$$ \hat{X} = \sum_{i=1}^{m} w_i Y_i \quad (2) $$

where $w_i$ is the weight of the sample ($X_i, Y_i$) in the same leaf of the bootstrap tree.

The final prediction of the new sample $D$ is approximated by the mean predictions of $b$ bootstrap trees.
Apart from the conditional mean estimate in RF, QRF also uses the weighted samples to derive a conditional distribution. This distribution function is able to provide the probability of Y being lower than a given percentile and thus to calculate the prediction intervals. For more details about the constructions of the conditional distribution, we refer to Meinshausen (2006).

We used the quantregForest package (Meinshausen, 2017) in R (R Core Team, 2019) for implementing QRF to derive the median prediction and 90% prediction intervals (90% PIs, 5th and 95th quantiles). The default number of tree (ntree=500) and minimum size of terminal nodes (nodesize=5) were used for QRF, and the number of variables randomly sampled as candidates at each split (mtry) was optimized in the caret package (Kunn, 2020) by 5-fold cross-validation in R (R Core Team, 2019). The variable importance was determined by the increased mean square error (IncMSE, in %) between the model excluding and including a given variable, and this information was integrated in QRF model. In our case, the variable importance was calculated by the average of 50 repeated models.

Considering the highly varying soil sampling density (Brus et al., 2011), we divided the study area into 20 compact equal area geographical strata (Figure 3) using the spcosa package (Walvoort et al., 2020) in R (R Core Team, 2019), and performed stratified random sampling (5 sites for each strata) for selecting the validation set. It resulted in a set of 10878 sites for model calibration and 100 sites for model validation. To derive a robust result, we repeated this procedure 50 times and took the average as the final model performance.

Four indicators were used to evaluate the model performance in validation set: (1) modelling efficiency (R²); (2) root mean square error (RMSE); (3) bias; (4) prediction interval coverage percentage (PICP), which describes the percentage of the observed TOL falls within the estimated upper and lower 90% PIs.

$$R^2 = 1 - \frac{\sqrt{\sum_{i=1}^{n} (\hat{y}_i - y_i)^2}}{\sqrt{\sum_{i=1}^{n} (y_i - \bar{y}_i)^2}} \quad (3)$$
$$270 \quad \text{RMSE} = \sqrt{\sum_{i=1}^{n} \frac{(\hat{y}_i - y_i)^2}{n}} \quad (4)$$

$$271 \quad \text{Bias} = \sum_{i=1}^{n} \frac{(\hat{y}_i - y_i)}{n} \quad (5)$$

where n is the sample size of observations, $y_i$ and $\hat{y}_i$ are observed value and predicted value for sample i, and $\bar{y}_i$ is the average of observed values.

In addition, we also reported the model performance by the internal validation using out-of-bag data (around 34% of data that is not used for growing each tree) in QRF. The final TOL map and its 90% PIs were produced by QRF fitted using all the sampling sites.

## 3. Results and discussion

### 3.1. Summary of TOL in the Beauce area

Table 3 displays the statistics of the TOL in the Beauce area. Among 10978 sites, TOL ranged from 0 to 120 cm, with a mean and median TOL at 36.46 cm and 30 cm, respectively. A skewness of 0.66 (larger than 0.5) indicated the data were slightly positively skewed while a kurtosis less than 3 (2.85) showed that the data were light-tailed. Though the log transformation, i.e. log(TOL+a), is able to convert soil data to normal distribution (a skewness less than 0.5 and a kurtosis close to 3), it did not improve the model performance. Therefore, instead of data transformation, we used the original TOL data for spatial modelling in this study.

Figure 4 presents the TOL located in 20 compact equal area geographical strata. It showed a large difference of number of sampling sites among these 20 geographical strata, ranging from 22 to 1351. These geographical strata with high median TOL (>60 cm) had much less sampling sites (51 to 201), and this is the main reason for evaluating the model performance by random stratified sampling.

Figure 5 shows the Pearson correlations coefficients between TOL and top 30 environmental covariates. Elevation and its focal mean derivatives had the highest positive correlations (r > 0.3). Other positively correlated covariates were related to
oblique coordinates, channel network base level (CNBL), graphic coordinates, gamma-ray radiometry and two Sentinel 2 indexes (grain size index and clay index). Negative correlations with TOL were found with valley depth, Sentinel 2 bands and its indexes. Overall, the correlations between TOL and covariates were found to be rather low (|r| < 0.35).

## 3.2. Model performance of Quantile Random Forest

The mean $R^2$, RMSE and bias from the internal validation using the out-of-bag data in QRF were 0.31, 26.88 cm and 0.23 cm (data not shown). Figure 6 indicates the model performance after we repeated 50 times the validation procedure using QRF. The mean $R^2$ and RMSE were 0.33 and 30.48 cm respectively. The mean bias of -1.20 cm indicated that the prediction was almost un-biased for 50 repeats. The mean PICP indicated that on average 86.7% of the validation samples fall within the defined 90% PIs, therefore the uncertainty estimates from the QRF model was valid for non-visited locations.

As shown in Table 4, the global soil thickness (ST) products (Hengl et al., 2017; Shangguan et al., 2017) had better model performance than those at national or regional scale (Guerrero et al., 2014; Kidd et al., 2015; Vaysse and Lagacherie, 2015; Lacoste et al., 2016; Mulder et al., 2016; Zhang et al., 2018). This may be attributed to the fact that global ST products include a substantial proportion of very thin soils (i.e., soils prone to severe erosion) and of very thick ones (i.e., Arenosols in desert dunes, Shangguan et al., 2017). There was no large difference of model performance between national and regional products, even if we used nearly 11000 sampling sites in this study. This is because the TOL is highly variable at short distances. By incorporating a large exhaustive set of environmental covariates, however, the map produced in this study performed slightly better than almost all the previous studies at regional and national scales.

The large range between upper and lower limits of 90% PIs for $R^2$ and RMSE indicated the randomness involved in data split brought a large amount of uncertainty in model evaluation. Therefore, instead of a single time data split, repeated random (stratified) sampling adopted in this study would provide more robust estimates for the model performance so as to avoid under- or over-optimistic decision making in management of soil resources.
329

## 3.3. Variable importance of environmental covariates

330

331

332

333

334

335

336

337

338

339

340

341

342

Figure 7 displays the top 30 environmental covariates in QRF model calculated as the average of 50 repeats. It indicated that the gamma radiometric data (U, Th and TC) and hue index (focal mean) calculated from Sentinel 2 images were the most important environmental covariates in modelling TOL in the study area. They were followed by longitude, NDVI PC1, slope position (with a radius of 1000 m), TWI (with a radius of 1000 m), normalized difference (focal mean) and grain size index (focal mean), representing position, organisms and relief factors in Scorpan conceptual model. For many relief (e.g., slope position, TWI, VDCN, curvature, valley depth, CNBL, elevation, slope) and organisms (e.g., hue index, normalized difference, grain size index, ferrous iron) factors, their derivatives calculated from neighbouring information performed better than original covariates. Spatial position, i.e., latitude and longitude, were identified important in Figure 7 while oblique coordinates were not listed among the top 30 covariates.

343

344

345

346

347

348

349

Interestingly, the variable importance in the QRF model was not in line with the correlations between TOL and covariates (see Figure 5). This may be due to the fact that the relations between TOL and covariates are not linear. If the relationships were linear then the most important covariates should have been those with highest |r| which is not the case in this study. Another reason may be that the importance of covariates results also from interactions between them, that are not visible using Pearson correlations but that are taken into account in QRF model.

350

351

352

353

354

355

356

357

358

359

360

361

Our results indicate a high importance of airborne gamma radiometric data in TOL modelling as they can capture soil information relevant to soil texture and to the presence of the calcareous rock at low depth. Indeed, the substrate of part of the study area (composed carbonates) is completely different from the TOL and it has been shown that calcium mitigates surface gamma-spectrometric signatures because it has a poor gamma-spectrometric response (Martelet *et al.*, 2013). Therefore, it is not surprising that gamma radiometric data plays an important role, especially for predicting thin TOL. Also, the large plateaus with deep TOL in the northern part are depleted in K (see Figure 2). This is because soils with large TOL were prone to illuviation, resulting in lower clay content in topsoil. So in this case, it is an indirect relationship with TOL. Our results also confirm the contribution of neighbouring information (e.g., focal hue index, slope position 1000m, TWI 1000m, focal normalized
difference) of relief and organism factors in spatial modelling of TOL, which implicates the multi-scale influence of covariates on soil properties. Concerning slope and TWI, the importance of this neighbouring information may be due to the gradient of TOL that shows that very large flat plateaus (mainly in the north) are characterized by a deeper TOL. These derivatives likely performed well because it is not the same geomorphological context if you have a flat location inside a very large flat plateau than if you have locally flat "pixels" in a region where the relief is more accentuated, such for instance in the southwest (Behrens et al., 2019). Other studies also have shown the potential of multi-scale covariates derivatives in improving model performance in DSM (Behrens et al., 2018b, 2019). Compared to simple convolution approach (focal mean), wavelet transforms, empirical mode decomposition, and the Gaussian scale space may even better represent the multi-scale information of environmental covariates (Behrens et al., 2010, 2018a, 2018b; Biswas et al., 2013a, 2013b; Zhou et al., 2016; Huang et al., 2017; Zhao et al., 2018; Liang et al., 2019) so as to improve model performance in DSM. The Sentinel 2 spectral bands may not always provide direct information related to soil, while a great potential has been shown from its derived indicators (e.g., NDVI, hue index, normalized difference, ferrous iron, and grain size index) in this study. Considering its high spatial and temporal resolution, Sentinel 2 has a great potential in delivering useful information of soil surface for DSM across scales (Gholizadeh et al., 2018; Castaldi et al., 2019; Loiseau et al., 2019; Vaudour et al., 2019). Some of the Sentinel 2 data we used come from a mosaic of images of bare soils built by Loiseau et al. (2019). Therefore, these Sentinel 2 data provide direct information on soil colour which may reflect thin TOL or absence of TOL through the presence of white calcareous rocks at the surface. They may indirectly reflect also texture through bright colours due to slaking that occurs mainly on very loamy topsoil soils which correspond to the deepest TOL where Luvisols have developed. The land use map produced by Sentinel 2 was not among the top 30 environmental covariates as it may be masked by the NDVI data due to their correlation or NDVI better explains the spatial variability than land use map. Therefore, relative importance of environmental covariates should be taken with caution as high contributing covariates can inadvertently bear part of the contribution of the less contributing covariates (Chen et al., 2018).
### 3.4. Maps of thickness of loess and its 90% prediction intervals

Figure 8 presents the spatial distribution of TOL and its lower and upper limits of 90% PIs. It displays the general increasing thickness of loess soils from south-west to north in the study area. Very shallow loess (<10 cm) was found in south-west of the study area, and very deep loess (>100 cm) was mainly found in the northern part. Highest TOL were mainly located in rather flat areas located on high elevation plateaus, while shallow loess was mainly located at lower elevations and in more dissected relief, especially in the vicinity of small valleys. Note that there is a border effect from the south-west to the west of the region. This border effect corresponds to the outcropping limit of the TOL, where sandy or clayey materials locally overlay the calcareous. The regions with thin soil (<10 cm with a lot of outcrops of the calcareous) correspond to the areas with black gamma-ray radiometry patterns matching on steep slopes around the drainage lines (mostly rivers).

The maps of lower and upper limits of 90% PIs clearly show different spatial structures. On the northern part with the highest elevations and high mean TOL, the 95th percentile is equal or deeper than 1.2 m, which means that high TOL are largely dominant in these plateaus. On the contrary the extreme southern part of the region exhibits TOL that rarely exceed 0.6 m. Moreover, except for some very local areas having a high mean TOL, the 5th percentile map suggests that the upper calcareous surface is undulating at very short distances and that local calcareous outcrops may be found in nearly all the southern part of the study area. The wind direction of loess deposits was from northwest (Bertran et al., 2016; Borderie et al., 2017). The Beauce area corresponds to the southern margin of the Paris basin loess deposits which show a clear gradient from North to South (Bertran et al., 2016). The gradient of loess that we observe in the Beauce region from north to south may be due to this. In addition, the northern part is characterized by large flat plateaus where no erosion occurred, except along the main deep valleys, whereas the southwestern part is characterized by a local relief that may have induced erosion and redistribution processes (Macaire, 1971). All these observations were confirmed by the expert knowledge of the soil surveyors who did some traditional reconnaissance soil mapping in this region. Interestingly, when doing reconnaissance maps at 1:250,000 the soil surveyors delineated small natural regions in order to create the broadest geographical ensembles of the legend (Richer-de-Forges, 2008; Richer-de-Forges et
427

428

429

430

431

432

433

434

435

al., 2008). Figure 9 shows these small natural regions drawn by the soil surveyors on the study area. The comparison between Figure 8 and Figure 9 clearly shows that the map of the TOL makes sense both from soil and physiography point of views. One should keep in mind that 90% PI is a very large PI. Therefore, it is normal that such wide ranges are found. Another reason for the large PI comes from the fact that our map does not have very high model performance, and there is still a large room to improve it. Useful outputs for irrigation or drainage management, however, are maps of probability of exceeding a given depth for TOL in the study area (see an example in the next section).

436

### 3.5. Example of application

437

438

439

440

441

442

443

444

One example of application is to map the probability of the TOL to exceed a given depth or, on the contrary, to map the probability of the TOL to be less than a given value. Figure 10 displays an example of these practical applications, which extracts the probability of exceeding of 30 cm from the function between the TOL and prediction quantile (from 0 to 100% with an interval of 2%) within the QRF model. The soils that have a very low probability to exceed a 30 cm TOL are unsuitable for conventional tillage and have a very low AWC. Therefore, optimizing the irrigation on these soils should greatly save water.

445

446

## 4. Conclusion

447

448

449

450

451

452

453

454

455

456

457

458

In this study, we utilized a large soil dataset (10978 sampling sites) and 117 environmental covariates relevant to soil, organisms, relief, parent material and spatial position for mapping thickness of loess at a regional scale. The 50 repeated Quantile Random Forest had an average R² of 0.33, which was slightly better than those obtained in most previous studies at regional or national scale (R² of 0.11~0.41). A PICP of 86.70% showed that around 86.70% of the validation samples fall within the predefined 90% PIs, which indicated that the prediction uncertainty produced by Quantile Random Forest was reasonable and can be properly used in decision making of land management. The relative importance of environmental covariates indicated the importance of elevation and gamma radiometry in modelling thickness of loess and also proved the necessity of incorporating neighbour information in relief and organisms for spatial modelling. The produced map of
459

thickness of loess and its 90% prediction intervals made sense from a soil science

460

perspective. This map can be further used for efficient irrigation management as well

461

as crop growth and yield modelling.

462

## Acknowledgement

463

Dominique Arrouays is coordinator, Vera Leatitia Mulder is member, and Anne

464

C. Richer-de-Forges, and Guillaume Martelet are collaborators of the Research

465

Consortium GLADSOILMAP supported by “LE STUDIUM” Loire Valley Institute for

466

Advances Research Studies. This work was funded by the French National Research

467

Agency (Soilserv program, ANR 16 CE32 0005 01). Songchao Chen has received

468

the support of China Scholarship Council for 3 years’ Ph.D. study in INRAE and

469

Agrocampus Ouest (under grant agreement no. 201606320211). We also

470

acknowledge Dr. Hélène Tissoux (BRGM, France) for sharing her knowledge on

471

loess deposits in France and the Beauce region.

472
## References

Arrouays, D., Grundy, M.G., Hartemink, A.E., Hempel, J.W., Heuvelink, G.B., Hong, S.Y., ..., Mendonca-Santos, M.d.L., 2014. GlobalSoilMap: toward a fine-resolution global grid of soil properties. In Advances in agronomy (Vol. 125, pp. 93-134). Academic Press.

Behrens, T., MacMillan, R. A., Viscarra Rossel, R.A., Schmidt, K., Lee, J., 2019. Teleconnections in spatial modelling. Geoderma, 354, 113854.

Behrens, T., Schmidt, K., MacMillan, R.A., Viscarra Rossel, R.A., 2018a. Multi-scale digital soil mapping with deep learning. Scientific reports, 8(1), 1-9.

Behrens, T., Schmidt, K., MacMillan, R.A., Viscarra Rossel, R.A., 2018b. Multiscale contextual spatial modelling with the Gaussian scale space. Geoderma, 310, 128-137.

Behrens, T., Zhu, A.X., Schmidt, K., Scholten, T., 2010. Multi-scale digital terrain analysis and feature selection for digital soil mapping. Geoderma, 155(3-4), 175-185.

Bertran, P., Liard, M., Sitzia, L., Tissoux, H., 2016. A map of Pleistocene aeolian deposits in Western Europe, with special emphasis on France. Journal of Quaternary Science, 31(8), 844-856.

Biswas, A., Cresswell, H.P., Chau, H.W., Viscarra Rossel, R.A., Si, B.C., 2013a. Separating scale-specific soil spatial variability: A comparison of multi-resolution analysis and empirical mode decomposition. Geoderma, 209, 57-64.

Biswas, A., Cresswell, H.P., Viscarra Rossel, R.A., Si, B.C., 2013b. Characterizing scale - and location - specific variation in non - linear soil systems using the wavelet transform. European Journal of Soil Science, 64(5), 706-715.

Borderie, Q., Chamaux, G., Roussaffa, H., Douard, M., Fencke, E., Rodot, M.-A., Perrichon, P., Selles, H., 2017. La couverture lœssique d'Eure-et-Loir (France): Potentiel pédo-sedimentaire et organisation spatiale. Quaternaire, 28(3), 389-400.
Bourennane, H., King, D., Chery, P., Bruand, A., 1996. Improving the kriging of a soil variable using slope gradient as external drift. European Journal of Soil Science, 47(4), 473-483.

Breiman, L., 2001. Random forests. Machine learning, 45(1), 5-32.

British Standards Institution, 1988. BSI Standards Catalogue: 1988. The Institution.

Brus, D.J., Kempen, B., Heuvelink, G.B.M., 2011. Sampling for validation of digital soil maps. European Journal of Soil Science, 62(3), 394-407.

Castaldi, F., Hueni, A., Chabrillat, S., Ward, K., Buttafuoco, G., Bomans, B., ..., van Wesemael, B., 2019. Evaluating the capability of the Sentinel 2 data for soil organic carbon prediction in croplands. ISPRS Journal of Photogrammetry and Remote Sensing, 147, 267-282.

Chen, S., Martin, M.P., Saby, N.P., Walter, C., Angers, D.A., Arrouays, D., 2018. Fine resolution map of top-and subsoil carbon sequestration potential in France. Science of the Total Environment, 630, 389-400.

Chen, S., Mulder, V.L., Martin, M.P., Walter, C., Lacoste, M., Richer-de-Forges, A.C. C., Arrouays, D., 2019. Probability mapping of soil thickness by random survival forest at a national scale. Geoderma, 344, 184-194.

Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., Böhner, J., 2015. System for Automated Geoscientific Analyses (SAGA) v. 2.1.4. Geoscientific Model Development, 8, 1991-2007.

Copernicus, 2016. Contains Modified Copernicus Service Information.
https://land.copernicus.eu/global/products/ndvi.

Doorenbos, J., Kassam, A.H., Bentvelder, C., Uittenboogaard, G., 1978. Yield response to water. FAO Irrigation and Drainage Paper 33. FAO, Rome, Italy, p. 144.

Gholizadeh, A., Žižala, D., Saberioon, M., Borůvka, L., 2018. Soil organic carbon and texture retrieving and mapping using proximal, airborne and Sentinel-2 spectral imaging. Remote Sensing of Environment, 218, 89-103.
Grinand, C., Arrouays, D., Laroche, B., Martin, M.P., 2008. Extrapolating regional soil landscapes from an existing soil map: sampling intensity, validation procedures, and integration of spatial context. Geoderma, 143(1-2), 180-190.

Grundy, M.J., Viscarra Rossel, R.A., Searle, R.D., Wilson, P.L., Chen, C., Gregory, L. J., 2015. Soil and landscape grid of Australia. Soil Research, 53(8), 835-844.

Guerrero, E., Pérez, A., Arroyo, C., Equihua, J., Guevara, M., 2014. Building a national framework for pedometric mapping: soil depth as an example for Mexico. In: GlobalSoilMap: Basis of the global spatial soil information system, Taylor & Francis, CRC Press, London 103-108.

Hagolle, O., Huc, M., Villa Pascual, D., Dedieu, G., 2015. A multi-temporal and multi-spectral method to estimate aerosol optical thickness over land, for the atmospheric correction of FormoSat-2, LandSat, VENUS and Sentinel-2 images. Remote Sensing, 7(3), 2668-2691.

Hengl, T., de Jesus, J.M., Heuvelink, G.B., Gonzalez, M.R., Kilibarda, M., Blagotić, A., Guevara, M.A., 2017. SoilGrids250m: Global gridded soil information based on machine learning. PLoS one, 12(2), e0169748.

Huang, J., Wu, C., Minasny, B., Roudier, P., McBratney, A.B., 2017. Unravelling scale-and location-specific variations in soil properties using the 2-dimensional empirical mode decomposition. Geoderma, 307, 139-149.

IGN, 2011. BD ALTI® Version 2.0 - Descriptif de Contenu. Available online at: http://professionnels.ign.fr/sites/default/files/DC_BDALTI_2-0.pdf

Inglada, J., Vincent, A., Arias, M., Tardy, B., Morin, D., Rodes, I., 2017. Operational high resolution land cover map production at the country scale using satellite image time series. Remote Sensing, 9(1), 95.

IUSS Working Group, WRB, 2006. World reference base for soil resources. World Soil Resources Report, 103.

Kidd, D., Webb, M., Malone, B., Minasny, B., McBratney, A., 2015. Eighty-metre resolution 3D soil-attribute maps for Tasmania, Australia. Soil Research, 53(8), 932-955.
King, D., Daroussin, J., Tavernier, R., 1994. Development of a soil geographic database from the Soil Map of the European Communities. CATENA, 21(1), 37-56.

Kuhn, M., 2020. caret: Classification and Regression Training. R package version 6.0-85. https://CRAN.R-project.org/package=caret

Lacoste, M., Mulder, V.L., Richer-de-Forges, A.C., Martin, M.P., Arrouays, D., 2016. Evaluating large-extent spatial modeling approaches: A case study for soil depth for France. Geoderma Regional, 7(2), 137-152.

Leenaars, J.G., Claessens, L., Heuvelink, G.B., Hengl, T., González, M.R., van Bussel, L.G., ..., Cassman, K.G., 2018. Mapping rootable depth and root zone plant-available water holding capacity of the soil of sub-Saharan Africa. Geoderma, 324, 18-36.

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.

Loiseau, T., Chen, S., Mulder, V. L., Román Dobarco, M., Richer-de-Forges, A.C., Lehmann, S., ..., Arrouays, D., 2019. Satellite data integration for soil clay content modelling at a national scale. International Journal of Applied Earth Observation and Geoinformation, 82, 101905.

Lombardo, L., Saia, S., Schillaci, C., Mai, P.M., Huser, R., 2018. Modeling soil organic carbon with Quantile Regression: Dissecting predictors' effects on carbon stocks. Geoderma, 318, 148-159.

Lorain, J.M., 1973. La géologie du calcaire de Beauce. Bulletin de Liaison des Laboratoires des Ponts et Chaussées. Special Issue, 78.

Macaire, J.M., 1971. Etude sédimentologique des formations superficielles sur le tracé de l'autoroute A10 entre Artenay et Meung sur Loire (PhD thesis). University of Orleans, France.

Martelet, G., Drufin, S., Tourlière, B., Saby, N. P. A., Perrin, J., DeParis, J., Prognon, J. F., Jolivet, C., Ratié, C., Arrouays, D., 2013. Regional regolith parameters
prediction using the proxy of airborne gamma ray spectrometry. Vadose Zone Journal, 12(4), 1-14.

Martelet, G., Nehlig, P., Arrouays, D., Messner, F., Tourlière, B., Laroche, B., Deparis, J., Saby, N., Richer-de-Forges, A.C., Jolivet, C., Ratié, C., 2014. Airborne gamma-ray spectrometry: potential for regolith-soil mapping and characterization. In: GlobalSoilMap: Basis of the global spatial soil information system. Taylor & Francis, CRC Press, London 401-408.

McBratney, A.B., Santos, M.d.L., Minasny, B., 2003. On digital soil mapping. Geoderma, 117(1-2), 3-52.

Meinshausen, N., 2006. Quantile regression forests. Journal of Machine Learning Research, 7, 983-999.

Meinshausen, N., 2017. quantregForest: Quantile Regression Forests. R package version 1.3-7.

Ménillet, F., 1974. Etude pétrographique et sédimentologique des calcaires d'Etampes et de Beauce. Formations dulcaquicoles du Stampien supérieur à l'Aquitanien dans le Bassin de Paris. Thesis, University of Paris-Sud (Orsay).

Møller, A.B., Beucher, A.M., Pouladi, N., Greve, M.H., 2019. Oblique geographic coordinates as covariates for digital soil mapping. SOIL Discussions, 1-20.

Mulder, V. L., Lacoste, M., Richer-de-Forges, A.C., Arrouays, D., 2016. GlobalSoilMap France: High-resolution spatial modelling the soils of France up to two meter depth. Science of the total environment, 573, 1352-1369.

Nicoulaud, B., Darthout, R., Duval, O., 1995. Etude de l'enracinement du blé tendre d'hiver et du maïs dans les sols argilo-limoneux de Petite Beauce. Etude et Gestion des Sols, 2, 183-200.

Ould-Mohamed, S., Bruand, A., Bruckler, L., Bertuzzi, P., Guillet, B., Raison, L., 1997. Estimating long-term drainage at a regional scale using a deterministic model. Soil Science Society of America Journal, 61, 1473-1482.

Paroissien, J. B., Saby, N., de Forges, A. R., Arrouays, D., Louis, B., 2014. Populating soil maps with legacy data from a soil testing databases. In:
GlobalSoilMap: Basis of the global spatial soil information system. Taylor & Francis, CRC Press, London 319-324.

R Core Team, 2019. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/

Richer-de-Forges, A.C., 2008. Notice explicative de la carte des pédopaysages du Loiret à 1/250 000 (Référentiel Régional Pédologique de la région Centre). INRA InfoSol. (in French)

Richer-de-Forges, A.C., Desbourdes, S., Lehmann, S., 2008. Référentiel Régional Pédologique de la région Centre : carte des pédopaysages du Loiret à 1/250 000. INRA InfoSol. (in French)

Shangguan, W., Hengl, T., de Jesus, J. M., Yuan, H., Dai, Y., 2017. Mapping the global depth to bedrock for land surface modeling. Journal of Advances in Modeling Earth Systems, 9(1), 65-88.

Tetegan, M., Richer-de-Forges, A.C., Verbeque, B., Nicoullaud, B., Desbourdes, C., Bouthier, A., Arrouays, D., Cousin, I., 2015. The effect of soil stoniness on the estimation of water retention. Catena, 129, 96-102.

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

Vaysse, K., Lagacherie, P., 2015. Evaluating digital soil mapping approaches for mapping GlobalSoilMap soil properties from legacy data in Languedoc-Roussillon (France). Geoderma Regional, 4, 20-30.

Vaysse, K., Lagacherie, P., 2017. Using quantile regression forest to estimate uncertainty of digital soil mapping products. Geoderma, 291, 55-64.

Walvoort, D., Brus, D., de Gruijter, J., 2020. Spatial Coverage Sampling and Random Sampling from Compact Geographical Strata. R package version 0.3-9. https://CRAN.R-project.org/package=spcosa

Zhang, W., Hu, G., Sheng, J., Weindorf, D.C., Wu, H., Xuan, J., ..., Gu, Z., 2018. Estimating effective soil depth at regional scales: Legacy maps versus
environmental covariates. Journal of Plant Nutrition and Soil Science, 181(2), 167-176.

Zhao, R., Biswas, A., Zhou, Y., Zhou, Y., Shi, Z., Li, H., 2018. Identifying localized and scale-specific multivariate controls of soil organic matter variations using multiple wavelet coherence. Science of The Total Environment, 643, 548-558.

Zhou, Y., Biswas, A., Ma, Z., Lu, Y., Chen, Q., Shi, Z., 2016. Revealing the scale-specific controls of soil organic matter at large scale in Northeast and North China Plain. Geoderma, 271, 71-79.
Figures

<figure><img src="image_2.png" /><figcaption>Figure 1 Study area and soil sampling sites, the Beauce area that locates at the middle Loire catchment, central France</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2 Examples of environmental covariates</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3 Compact equal area geographical strata</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4 Boxplot of TOL for each compact equal area geographical strata. The number of sampling sites is indicated in blue.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5 Pearson correlation coefficient (r) between TOL and top correlated environmental covariates (r > 0.15 or r < -0.15).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6 Model performance of 50 repeats evaluated by R² (a), RMSE (b), bias (c) and coverage of PICP (d) on validation set</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7 Relative importance of environmental covariates in Quantile Random Forest (average of 50 repeats). Only the top 30 variables are shown here.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 8 Spatial distribution of the thickness of loess and its 90% prediction intervals</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9 Natural regions delineated by soil surveyors (Richer-de-Forges, 2008; Richer-de-Forges *et al.*, 2008)</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10 The probability of exceeding 30 cm for TOL in the study area. The probability at 30 cm is extracted from the probability distribution of Quantile Random Forest for each pixel.</figcaption></figure>
# Tables

<table><caption>Table 1 Environmental covariates used for digital soil mapping</caption><thead><tr><th>Variable</th><th>Number</th><th>Resolution</th><th>Scorpan factor</th><th>Reference</th></tr></thead><tbody><tr><td>Channel network base level</td><td>4</td><td>25 m</td><td>Relief<sup>a</sup></td><td>IGN, 2011</td></tr><tr><td>Elevation</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Multiresolution index of valley bottom flatness</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Plan curvature</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Profile curvature</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Slope</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Slope position</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Slope length</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Terrain wetness index</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Valley depth</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Vertical distance to channel network</td><td>4</td><td>25 m</td><td>Relief</td><td>IGN, 2011</td></tr><tr><td>Soil type</td><td>1</td><td>1:1000000</td><td>Soil</td><td>King et al. (1995)</td></tr><tr><td>Parent material</td><td>1</td><td>1:1000000</td><td>Parent material</td><td>King et al. (1995)</td></tr><tr><td>Gamma radiometric (K, U, Th, TC)</td><td>4</td><td>200 m</td><td>Parent material</td><td>Martelet et al. (2014)</td></tr><tr><td>Land cover and probability</td><td>10</td><td>10, 20, 60 m</td><td>Organisms</td><td>Inglada et al. (2017)</td></tr><tr><td>Sentinel 2 spectral bands and indices</td><td>42</td><td>90 m</td><td>Organisms</td><td>Loiseau et al. (2019)</td></tr><tr><td>First three PCs of monthly NDVI<sup>b</sup></td><td>3</td><td>300, 500 m</td><td>Organisms</td><td>Loiseau et al. (2019)</td></tr><tr><td>Coordinates (Latitude, Longitude)</td><td>2</td><td>25 m</td><td>Position</td><td>IGN, 2011</td></tr><tr><td>Oblique coordinates<sup>c</sup></td><td>10</td><td>25 m</td><td>Position</td><td>Møller et al. (2019)</td></tr></tbody></table>

<footnote><sup>a</sup> For all the covariates in relief factor, except for the original products, their local mean values with radius at 200, 500 and 1000 m are also calculated by convolution circular windows. <sup>b</sup> PCs, principal components; NDVI, normalized difference vegetation index. <sup>c</sup> Oblique coordinates at angles of 15°, 30°, 45°, 60°, 75, 105°, 120°, 135°, 150° and 165° are produced</footnote>
<table><caption>Table 2 Random Forest model performance over the predictor subset size using recursive feature selection</caption><thead><tr><th>Subset size</th><th>RMSE</th><th>R<sup>2</sup></th><th>Selected</th></tr></thead><tbody><tr><td>5</td><td>28.70</td><td>0.1808</td><td>No</td></tr><tr><td>10</td><td>27.07</td><td>0.2712</td><td>No</td></tr><tr><td>15</td><td>26.72</td><td>0.2904</td><td>No</td></tr><tr><td>20</td><td>26.52</td><td>0.3007</td><td>No</td></tr><tr><td>40</td><td>26.50</td><td>0.3022</td><td>No</td></tr><tr><td>60</td><td>26.50</td><td>0.3019</td><td>No</td></tr><tr><td><strong>80</strong></td><td><strong>26.44</strong></td><td><strong>0.3052</strong></td><td><strong>Yes</strong></td></tr><tr><td>100</td><td>26.50</td><td>0.3025</td><td>No</td></tr><tr><td>117</td><td>26.48</td><td>0.3038</td><td>No</td></tr></tbody></table>
<table><caption>Table 3 Statistics of the thickness of loess (in cm)</caption><thead><tr><th>Variable</th><th>Number</th><th>Minimum</th><th>Q1</th><th>Median</th><th>Mean</th><th>Q3</th><th>Maximum</th><th>Skewness</th><th>Kurtosis</th></tr></thead><tbody><tr><td>TOL</td><td>10978</td><td>0</td><td>0</td><td>30</td><td>36.46</td><td>60</td><td>120</td><td>0.66</td><td>2.85</td></tr></tbody></table>

Q1, the first quartile; Q3, the third quartile.
<table><caption>Table 4 A summary of the model performance of soil thickness (or soil depth) mapping from regional to global scales</caption><thead><tr><th>Reference</th><th>Location</th><th>R<sup>2</sup></th></tr></thead><tbody><tr><td colspan="3"><strong>Regional scale</strong></td></tr><tr><td>Kidd et al. (2015)</td><td>Tasmania, Australia</td><td>0.16</td></tr><tr><td>Vaysse and Lagacherie (2015)</td><td>Languedoc-Roussillon, France</td><td>0.23</td></tr><tr><td>Zhang et al. (2018)</td><td>Xinjiang, China</td><td>0.28</td></tr><tr><td>This study</td><td>Beauce, France</td><td>0.34</td></tr><tr><td colspan="3"><strong>National scale</strong></td></tr><tr><td>Guerrero et al. (2014)</td><td>Mexico</td><td>0.41</td></tr><tr><td>Lacoste et al. (2016)</td><td>France</td><td>0.22</td></tr><tr><td>Mulder et al. (2016)</td><td>France</td><td>0.11</td></tr><tr><td colspan="3"><strong>Global scale</strong></td></tr><tr><td>Hengl et al. (2017)</td><td>Globe</td><td>0.57</td></tr><tr><td>Shangguan et al. (2017)</td><td>Globe</td><td>0.59</td></tr></tbody></table>