# Does digital soil mapping prediction performance of soil texture improve when adding uncertain field texture estimates? A study based on clay content

Anne C Richer-De-Forges, Songchao Chen, Gerard B M Heuvelink, Stephan van der Westhuizen, Thomas G Orton, Hocine Bourennane, Dominique Arrouays

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

## ► To cite this version:

Anne C Richer-De-Forges, Songchao Chen, Gerard B M Heuvelink, Stephan van der Westhuizen, Thomas G Orton, et al.. Does digital soil mapping prediction performance of soil texture improve when adding uncertain field texture estimates? A study based on clay content. Geoderma, 2025, 456, pp.117277. 10.1016/j.geoderma.2025.117277. hal-05014984

HAL Id: hal-05014984

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

Submitted on 1 Apr 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
# Does digital soil mapping prediction performance of soil texture improve when adding uncertain field texture estimates? A study based on clay content

Anne C. Richer-de-Forges<sup>a, *, ID</sup>, Songchao Chen<sup>b, c ID</sup>, Gerard B.M. Heuvelink<sup>d, e ID</sup>, Stephan van der Westhuizen<sup>f</sup>, Thomas G. Orton<sup>g ID</sup>, Hocine Bourennane<sup>a</sup>, Dominique Arrouays<sup>a ID</sup>

## ARTICLE INFO

Handling Editor: Jingyi Huang

Keywords:

Hand-feel soil texture

Uncertainty

Quantile regression forest

Prediction maps

Particle-size distribution

Laboratory measurements

Measurement error

## ABSTRACT

Hand-feel soil texture observations (HFST) are less accurate, yet much more numerous, than laboratory measurements of soil texture (LAST). Therefore, it is tempting to incorporate both LAST and HFST information as calibration data in digital soil mapping (DSM) of particle-size distribution. We used about 1000 LASTs and 15,000 HFSTs over an area of about 6,800 km². We incorporated the uncertainties of HFST and LAST calibration data in DSM and compared it with a case where measurement errors were ignored. We added progressively HFST calibration data to LAST data and ran predictions and k-fold validations keeping the same validation set for all experiments. We added HFSTs according to different strategies: either based on the most uncertain predicted areas from the LAST-only model, or those from the preceding LAST + nHFST model, or randomly. We discuss the pros and cons of these different strategies. Adding HFST data brought useful information for model calibration, but only if the uncertainty was accounted for. Various strategies for adding HFSTs led to different unbalanced samplings, maps and prediction intervals. We explain how these various unbalanced samplings sharpened or enlarged the predictive distribution of various clay content ranges. Adding a too large number of HFSTs led to an over-optimistic estimation of the 90 % prediction interval and to large homogeneous patterns, smoothing the spatial variation of clay content. Adding HFSTs using weights to acknowledge their uncertainty substantially improved DSM predictions, but the number of HFSTs and the strategy to add them must be carefully adapted.

## 1. Introduction

Soil texture (ST) is the numerical proportion (weight in g 100 g⁻¹ or g kg⁻¹) of the sand, silt, and clay separates of the fine earth fraction (≤2 mm) (NRCS-USDA, 2012). The size thresholds for separating the sand, clay and size fractions vary among international classifications (e. g., IUSS Working Group WRB, 2015; Takahashi et al., 2020; USDA-NRCS, 2014) and countries (e.g., Minasny et al., 2007; Minasny and McBratney, 2001; Nemes et al., 1999; Rousseva, 1997). Whatever the thresholds of the fractions are, ST can either be measured in the

laboratory (laboratory analysis soil texture: LAST), most frequently by the pipette method (AFNOR, 2003), or estimated in the field by hand-feel (HFST) from the soil surveyor according to practical guidelines, such as those from Thien (1979) and Ritchey et al. (2015). When ST is estimated by HFST, the soil surveyor allocates the observations to ST classes defined by limits of sand, silt, and clay relative proportions. These classes can be represented in a texture triangle, either rectangular or equilateral. Many ST triangles have been published and used (e.g., Richer-de-Forges et al., 2008a) and their number of classes range from three to thirty-seven (Richer-de-Forges et al., 2022). ST triangles are also

<footnote><sup>a</sup> INRAE, Info&Sols Unit, 45075 Orléans, France</footnote>

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

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

<footnote><sup>d</sup> Soil Geography and Landscape Group, Wageningen University, Wageningen, the Netherlands</footnote>

<footnote><sup>e</sup> ISRIC – World Soil Information, Wageningen, the Netherlands</footnote>

<footnote><sup>f</sup> Stellenbosch University, Department of Statistics & Actuarial Science, Stellenbosch, South Africa</footnote>

<footnote><sup>g</sup> The University of Queensland, School of Agriculture and Food Sciences, St Lucia, Queensland 4072, Australia</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: anne.richer-de-forges@inrae.fr (A.C. Richer-de-Forges).</footnote>
used to plot measured values of particle-size distribution.

In conventional soil surveys, because of cost and time constraints, HFSTs are generally much more numerous than LASTs. In practice, when collecting information during survey campaigns, the soil surveyor does a large number of borings to extract soil samples on which he/she estimates HFST classes. Meanwhile, or in a second step, the soil surveyor selects much fewer locations, to describe soil horizons from a soil pit and sample them for LAST determinations. A HFST of these samples is also often determined in the field, prior to sending them to the lab.

In digital soil mapping (DSM; e.g., McBratney et al., 2003; Minasny and McBratney, 2016), point-scale information is often used for calibrating predictive models of soil attributes. As pointed out by numerous authors (e.g., Arrouays et al., 2020; Grunwald, 2009; Hartemink et al., 2008; Loiseau et al., 2021; Malone et al., 2018; McBratney et al., 2003; Somarathna et al., 2017), a major source of error in DSM is the sparseness of soil calibration data, which justifies the worldwide, continental and national efforts to obtain new data or to rescue legacy data (e.g., Armas et al., 2023; Arrouays et al., 2017; Batjes et al., 2017; Cooper et al., 2005; de Sousa et al., 2020; Guevara et al., 2020, 2018; Hugelius et al., 2013; Leenaars, n.d; Medeiros et al., 2024; Orgiazzi et al., 2018; Panagos et al., 2012; Pfeiffer et al., 2020; Rasaei et al., 2020; Simões, 2015; Sulaeman et al., 2013). One issue that arises in DSM using machine learning is how to incorporate uncertain soil information/measurements in the modelling process (Chen et al., 2023; Hengl et al., 2018; Takoutsing and Heuvelink, 2022; van der Westhuizen et al., 2022; Wadoux et al., 2020). Indeed, for DSM, we usually assume that soil calibration data are error-free or disregard measurement errors in soil calibration data assuming that they are identically distributed. But often that is not a realistic assumption, particularly when we have calibration data from different sources, such as high-quality laboratory data and low-quality field data (e.g., Malone and Searle, 2021a, 2021b). This is also the case for the soil texture example that we study here: we have both LAST and HFST data for calibration of a DSM model. We aim at using both calibration data sources for mapping soil texture (i.e. clay, sand and silt percentages) while accounting for the fact that the LAST data are more accurate than the HFST data.

Indeed, it is tempting to incorporate both LAST and HFST

information as calibration data (Malone and Searle, 2021a, 2021b). In that case, we must take into account measurement uncertainty of both LAST and HFST. As HFST is less accurate than LAST, incorporating HFST might result in adding more uncertainty without substantially improving the predictions. Not accounting for measurement errors would consider each individual HFST measurement as important as each individual LAST measurement. In practice, this would mean that LAST measurements play almost no role because they are much less numerous than HFST; but that would be wrong. Indeed, there are still many open questions when facing the problem of DSM using data having different sources of uncertainties. Some of them are listed in the Challenge 3 of the article"Ten challenges for the future of pedometrics" (Wadoux et al., 2021). In this study, we will concentrate on incorporating measurement uncertainty in LAST and HFST calibration data and on four objectives. These objectives are to: i) quantify measurement uncertainty of HFST and LAST calibration data and incorporate this in DSM model training; ii) evaluate whether adding HFST calibration data improves DSM prediction; iii) compare different sampling strategies for adding progressively more HFST calibration data to LAST data; and iv) analyse the possible effect of not accounting for measurement uncertainty and treating LAST and HFST calibration data as equally important. In this paper, we consider the example of predicting solely clay content (0—2 µm) of A horizons, for the sake of simplicity.

## 2. Material and methods

### 2.1. Study area

The study area is the French department of Loiret (Fig. 1), which covers 6,775 km² in the French Region Centre - Val de Loire. Loiret is one of the flattest departments in France, with a range in altitude from 70 m above sea level in the northeastern area to 270 m in the southeast (Richer-de-Forges et al., 2008b). It has a continental-oceanic climate with a mean annual temperature of 11.5 °C and a mean annual rainfall of 700 mm. It is crossed from South-East to West by the Loire Valley. The soils of Loiret show high spatial variability owing to the great diversity of lithology and derived parent materials (Lacquement et al., 2015; Liard

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Location of Loiret and location of the soil observations. LAST: laboratory analysis soil texture; HFST: hand-feel soil texture. All observations are A horizons. LAST and HFST data include 1,078 and 15,182 samples respectively.</figcaption></figure>
et al., 2017; Richer-de-Forges et al., 2008b; Tissoux et al., 2017). The primary soils observed during a 1:250,000 scale soil survey were Cam-bisols, Luvisols, Planosols, Fluvisols, Colluviosols, Arenosols, and Ver-tisols (Richer-de-Forges et al., 2008b). More locally, impervious clayey subsoils located in depressions can lead to the development of Gleysols and Histosols.

## 2.2. Soil data

Our target variable is the clay content of the A horizon. Soil data were extracted from the French national database DoneSol (Grolleau et al., 2004; InfoSol, 2021). All topsoil horizons were defined as organo-mineral A horizons (excluding O horizons when present). For the LASTs, we selected A horizons starting at 0 cm (or just below the O horizon) and presenting a particle-size analysis without carbonate removal carried out in the laboratory. Samples were air-dried, then gently crushed and passed through a 2-mm sieve. Particle-size distribution was measured according to the pipette method, which is the French standard (AFNOR, 2003). For the HFSTs, we selected A horizons starting at 0 cm and ending before 30 cm. All soil surveyors were experienced and used the same protocols and a standard texture triangle to describe the soils in the field (HFST). The ‘Aisne’ triangle defined by Jamagne (1967) was used, and encoded in the soiltexture package in R software (version 1.5.1., Moeys et al., 2018). In total, 1,078 A horizons both with LAST and HFST, and 15,182 A horizons without LAST but with HFST were used (Fig. 1).

## 2.3. Converting HFST observations to LAST proxies by taking uncertainty into account

### 2.3.1. Converting HFST observations to LAST proxies by using a set of paired observations

To determine LAST proxies for HFST data in the approach that accounts for their measurement errors, we used all soil samples from A horizons available in the Loiret dataset having both HFST and LAST information to extract a set of paired observations. From this paired dataset, we determined the centroid of all LAST observations that had a corresponding HFST value of that same class. We then allocated to each horizon having solely a HFST class a particle-size distribution (or LAST proxy) equal to the centroid of LAST observations for the class. These centroids of LAST values were calculated with the centroid function of the R package rearrr (version 0.3.3., Olsen, 2020). Fig. 2 displays the LAST centroids allocated to each HFST class of the Aisne texture triangle using the R package soiltexture (version 1.5.1., Moeys et al., 2018). From these centroids, we derived a clay content proxy for each HFST class. Note that for two HFST classes (LMS and LM; Fig. 2) the centroid was located outside the HFST class, though their clay content proxies remained within the range of theoretical values of the LMS and LM classes defined by the clay-content boundaries of those classes.

Fig. 3 displays density plots of measured clay content (LAST) for each HFST from the paired data. In general, the sandiest fractions exhibited rather peaked distributions (e.g., S, SL, LS, LSA, SA) whereas the other HFST had flatter LAST distributions representing larger uncertainties, with a noticeable exception for silty clay (AL).

### 2.3.2. Estimation of laboratory measurement error (LAST)

The INRAE Arras soil sample analysis laboratory quantifies the LAST measurement errors in terms of the lower and upper limits of a 95 % confidence interval. For the data used in this study these limits were specified as the measured clay fraction minus and plus 4 % of the measured value, respectively (Baize, 2018). Assuming these measurement errors to be normally distributed, the measurement error variance ($MEV_{LAST}$) is thus given by Equation (1).

$$MEV_{LAST} = \left( \frac{0.04 \times CLAY}{1.96} \right)^2 \#(1)$$

<figure><img src="image_4.png" /><figcaption>Fig. 2. Centroids of the particle-size distributions (LASTs) of each hand-feel soil texture (HFST) class as projected in the soil texture triangle from Jamagne (1967). ALO: heavy clay, A: clay, AL: silty clay, AS: sandy clay, LA: clayey silt, LAS: sandy clayey silt, LSA: clayey sand silt, SA: clayey sand, S: sand, SL: silty sand, LL: silt, LS: sandy silt, LMS: sandy medium silt, LM: medium silt, LLS: sandy silt. Note that LLS and LL textures were never observed.</figcaption></figure>

where CLAY is the measured clay content in %.

### 2.3.3. Quantification of measurement errors in HFST data

We used the numerous paired laboratory and field measurements to derive a variance of clay content (Table 1) within each HFST class. We assumed that this variance is constant for all measurements in the same HFST class, and we estimated this variance by computing the variance of all measurements that are in the same class. Note that this variance only quantifies the random component of the HFST measurement error, not the systematic error component. However, using the centroid of the LAST measurements within a HFST class as the estimate, we basically corrected for systematic errors. Note also that this is valid under the assumption that LAST measurement errors are negligible when compared to HFST measurement errors.

For HFST data, the $MEV_{HFST}$ is the variance of the differences between LAST and HFST observations. $MEV_{HFST}$ was on average much bigger than $MEV_{LAST}$. In addition, the distribution of $MEV_{LAST}$ showed a much narrower interquartile range than that of $MEV_{HFST}$. This result confirms that LAST measurement errors are negligible when compared to HFST measurement errors.

## 2.4. Converting HFST observations to LAST proxies without taking uncertainty into account

For the approach that disregards measurement errors, we wanted to evaluate if and how much damage is done by ignoring measurement uncertainty and treating LAST and HFST calibration data equally important. Therefore, we also considered a case in which the uncertainty of LAST and HFST calibration data is unknown and no weights are assigned. In particular, this case may happen when the producer does not have enough paired data to recalculate a centroid of measured LAST values for each HFST class as we did in Section 2.3.1. In this case, we retained a widely used strategy, i.e. using the geometric centroids of the HFST class as defined in the soil texture triangle (e.g., Levi, 2017). We calculated the geometric centroids of the classes with the soiltexture
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Density plots of measured clay content (LAST) for each hand-feel soil texture (HFST) class. ALO: heavy clay, A: clay, AL: silty clay, AS: sandy clay, LA: clayey silt, LAS: sandy clayey silt, LSA: clayey sand silt, SA: clayey sand, S: sand, SL: silty sand, LL: silt, LS: sandy silt, LMS: sandy medium silt, LM: medium silt, LLS: sandy silt. Note that LLS and LL textures were never observed.</figcaption></figure>

<table><caption>**Table 1**
Clay content (LAST, in %, g 100 g⁻¹) summary statistics by HFST class in the Loiret.</caption><thead><tr><th>HFST</th><th>Number of LAST (N)</th><th>Min.</th><th>1st Qu.</th><th>Median</th><th>Mean</th><th>3rd Qu.</th><th>Max.</th><th>MEV<sub>HFST</sub> = VarianceS<sup>2</sup><sub>HFST</sub></th><th>centroid</th></tr></thead><tbody><tr><td>ALO</td><td>64</td><td>14.30</td><td>44.90</td><td>50.50</td><td>49.84</td><td>56.23</td><td>81.60</td><td>136.45</td><td>49.84</td></tr><tr><td>A</td><td>65</td><td>16.14</td><td>30.80</td><td>35.50</td><td>35.29</td><td>40.30</td><td>57.10</td><td>61.31</td><td>35.29</td></tr><tr><td>AL</td><td>110</td><td>8.60</td><td>32.20</td><td>34.05</td><td>33.89</td><td>37.20</td><td>59.20</td><td>53.09</td><td>33.89</td></tr><tr><td>AS</td><td>41</td><td>17.19</td><td>27.00</td><td>30.50</td><td>30.14</td><td>34.87</td><td>42.70</td><td>34.70</td><td>30.14</td></tr><tr><td>LA</td><td>110</td><td>10.70</td><td>27.05</td><td>29.50</td><td>28.71</td><td>31.67</td><td>42.40</td><td>31.58</td><td>28.71</td></tr><tr><td>LAS</td><td>36</td><td>15.70</td><td>20.10</td><td>24.85</td><td>24.20</td><td>27.25</td><td>36.80</td><td>24.87</td><td>24.20</td></tr><tr><td>LM</td><td>36</td><td>6.00</td><td>13.01</td><td>17.30</td><td>17.61</td><td>21.57</td><td>33.20</td><td>42.00</td><td>17.61</td></tr><tr><td>LMS</td><td>21</td><td>11.10</td><td>14.70</td><td>16.00</td><td>16.15</td><td>17.30</td><td>22.50</td><td>7.19</td><td>16.15</td></tr><tr><td>LS</td><td>29</td><td>6.82</td><td>11.10</td><td>13.30</td><td>13.99</td><td>15.82</td><td>26.50</td><td>16.06</td><td>13.99</td></tr><tr><td>LSA</td><td>81</td><td>5.46</td><td>19.00</td><td>20.50</td><td>20.91</td><td>23.00</td><td>29.80</td><td>20.16</td><td>20.91</td></tr><tr><td>S</td><td>157</td><td>0.00</td><td>3.33</td><td>5.17</td><td>5.30</td><td>6.81</td><td>15.20</td><td>7.14</td><td>5.30</td></tr><tr><td>SA</td><td>110</td><td>5.73</td><td>14.53</td><td>16.75</td><td>17.07</td><td>19.15</td><td>29.29</td><td>18.58</td><td>17.07</td></tr><tr><td>SL</td><td>135</td><td>1.35</td><td>7.20</td><td>9.00</td><td>8.88</td><td>10.94</td><td>17.20</td><td>8.34</td><td>8.88</td></tr></tbody></table>

package in R (version 1.5.1., Moeys et al., 2018).

## 2.5. Environmental covariates

The environmental covariates used included a lithology map, airborne gamma-ray spectrometry data and relief covariates derived from a digital elevation model (DEM) at 25 m resolution of the French National Geographical Institute (IGN, 2014) (Table 2 and Fig. S1). The relief covariates were extracted (elevation) or calculated (others) at 25 m resolution and then rescaled to 90 m resolution. Elevation and slope were also extracted/calculated at 25 m resolution and their focal means and standard deviations were calculated over circles of increasing diameter of 100, 500 and 1,000 m. These focal means and standard deviations were allocated to the original 25 m pixels, and then rescaled to 90 m. The airborne gamma spectrometry survey was conducted in 2009 (Terraquest LTD, 2009). Equipment characteristics are described in detail in Martelet et al. (2014). The flights were undertaken at an elevation of 80 m along a regular flight path consisting of flight lines spaced 1 km apart in the N-S direction and perpendicular lines spaced 10 km apart. The footprint of each measurement was approximately 150 × 300 m. The total count (absorbed dose rate in nanogray/h) and concentrations of uranium (ppm), thorium (ppm), and potassium (g 100 g⁻¹) were derived from the raw spectra following International Atomic Energy Agency (IAEA) processing guidelines (IAEA, 2003, 1991). The corrected gamma spectrometric data were then interpolated by the Bureau de Recherches Géologiques et Minières (BRGM) on regular 250 m cell grids using a standard minimum curvature interpolation. The interpolated gamma-ray data showed some gaps, in particular above a couple of built-up areas and a military airport. To fill these gaps, we modeled gamma-ray data using RF and other available covariates described in Table 2. We acknowledge that imputing missing data in this way adds some uncertainty to gamma-ray data and may degrade DSM predictions as the filled areas contain some LAST and HFST data points.

However, very few observation points (about 3 %, either LAST or HFST) were present in these areas that were mainly urban or had been disturbed by bombing during World War II.

The lithology of Loiret (Martelet et al., 2014) was grouped into 7 classes (“alluvium + gravel + silty clay sand”, “flint clay + flint + flint chalk + chalk + clay”, “limestone marl”, “limestone”, “sand”, “silt” and “clay sand”); (Fig. S2). The grouping was done taking into account the number of pixels in each class and their geological similarity. This covariate was binarized to convert this categorical variable to a set of seven binary variables (0,1) each indicating presence/absence of a single class.

A first model was run using LAST values to observe the importance of the covariates in the prediction of clay. We used RF default hyperparameters, i.e., mtry = p/3 (where p is number of covariates); node-size = 5; ntree = 500. We used recursive feature elimination (RFE) to select the most relevant covariates and their optimal number. The 22 selected covariates are presented in Table 2.

## 2.6. Modelling

### 2.6.1. Modelling for obtaining the weights

We used a RF model (Breiman, 2001) with the methodology proposed in van der Westhuizen et al. (2022) to model uncertain clay content. This method uses an iterative procedure to obtain an estimate of the residual variance (RV) of the model. Weights (wt) are then obtained with:

$$wt = \frac{1}{(RV + MEV)^{\#(2)}}$$

Weights are used with the calibration data so that a measurement with a larger error variance receives a smaller weight. Note that in Equation (2) MEV is the measurement error variance, that is MEV<sub>LAST</sub> for LAST data and MEV<sub>HFST</sub> for HFST data. We implemented this methodology with the
<table><caption>**Table 2**
List of environmental covariates after recursive feature elimination.</caption><thead><tr><th>name of covariable</th><th>Description</th><th>Resolution</th><th>Unit</th><th>References</th></tr></thead><tbody><tr><td>mnt</td><td>elevation</td><td>90 × 90</td><td>m</td><td rowspan="2">DEM (IGN, 2014)</td></tr><tr><td>mnt100m</td><td>focal length mean of elevation on a circle of 100 m diameter</td><td>90 × 90</td><td>m</td></tr><tr><td>mnt500m</td><td>focal length mean of elevation on a circle of 500 m diameter</td><td>90 × 90</td><td>m</td><td></td></tr><tr><td>mnt500e</td><td>focal length standard deviation of elevation on a circle of 500 m diameter</td><td>90 × 90</td><td>m</td><td></td></tr><tr><td>mnt1kmm</td><td>focal length mean of elevation on a circle of 1 km diameter</td><td>90 × 90</td><td>m</td><td></td></tr><tr><td>mnt1kme</td><td>focal length standard deviation of elevation on a circle of 1 km diameter</td><td>90 × 90</td><td>m</td><td></td></tr><tr><td>dppr</td><td>hydrological distance from each pixel to its nearest outlet</td><td>90 × 90</td><td>m</td><td></td></tr><tr><td>beven</td><td>Topographic wetness index (Beven and Kirkby, 1979)</td><td>90 × 90</td><td></td><td></td></tr><tr><td>courb_mixte</td><td>mixed curvature</td><td>90 × 90</td><td>m</td><td></td></tr><tr><td>pente100m</td><td>focal length mean of slope on a circle of 100 m diameter</td><td>90 × 90</td><td>degree</td><td></td></tr><tr><td>pente100e</td><td>standard focal length standard deviation of slope on a circle of 100 m diameter</td><td>90 × 90</td><td>degree</td><td></td></tr><tr><td>pente500m</td><td>focal length mean of slope on a circle of 500 m diameter</td><td>90 × 90</td><td>degree</td><td></td></tr><tr><td>pente500e</td><td>focal length standard deviation of slope on a circle of 500 m diameter</td><td>90 × 90</td><td>degree</td><td></td></tr><tr><td>pente1kmm</td><td>mean focal length mean of slope on a circle of 1 km diameter</td><td>90 × 90</td><td>degree</td><td></td></tr><tr><td>pente1kme</td><td>focal length standard deviation of slope on a circle of 1 km diameter</td><td>90 × 90</td><td>degree</td><td></td></tr><tr><td>litho_sable</td><td>lithology = sand</td><td>90 × 90</td><td>binary<br>(0; 1)</td><td rowspan="4">BRGM (Martelet et al., 2014)</td></tr><tr><td>litho.marne</td><td>lithology = limestone marl</td><td>90 × 90</td><td>binary<br>(0; 1)</td></tr><tr><td>litho_limon</td><td>lithology = lœss deposit</td><td>90 × 90</td><td>binary<br>(0; 1)</td></tr><tr><td>u</td><td>uranium (gamma-ray)</td><td>90 × 90</td><td>ppm<br>(mg kg<sup>-1</sup>)</td></tr><tr><td>k</td><td>potassium (gamma-ray)</td><td>90 × 90</td><td>g 100<br>g<sup>-1</sup></td><td></td></tr><tr><td>th</td><td>thorium (gamma-ray)</td><td>90 × 90</td><td>ppm<br>(mg kg<sup>-1</sup>)</td><td></td></tr><tr><td>tot</td><td>total gamma-ray count</td><td>90 × 90</td><td>nGy h<sup>-1</sup></td><td></td></tr></tbody></table>

ranger package (version 0.15.1., Wright and Ziegler, 2017) in R with the case weights functionality. Note that this argument assigns weights to observations to be selected for the bootstrap sample to construct the trees. However, in van der Westhuizen et al. (2022) it is explained that this is similar to the proposed methodology.

Using this method, we obtained for LAST and HFST the range of weights described in Table 3. Note that summary statistics suggest that the weights of HFST data are about half of the weights of LAST data.

Then, we applied weights in a quantile regression forest (QRF) model to produce prediction intervals, with the aim of comparing whether

<table><caption>**Table 3**
Summary statistics of weights values for clay content (LAST) and hand-feel soil texture (HFST) estimates.</caption><thead><tr><th>Weight</th><th>Min.</th><th>1st Qu.</th><th>Median</th><th>Mean</th><th>3rd Qu.</th><th>Max.</th></tr></thead><tbody><tr><td>LAST data</td><td>0.011</td><td>0.021</td><td>0.028</td><td>0.030</td><td>0.036</td><td>0.059</td></tr><tr><td>HFST data</td><td>0.006</td><td>0.012</td><td>0.016</td><td>0.017</td><td>0.018</td><td>0.031</td></tr></tbody></table>

incorporating weights would reduce the uncertainty in the predictions. It should be noted that the quantiles produced from the QRF model will reflect the conditional distribution of uncertain clay and not that of the error-free clay variable. This is because the weights defined in Equation (2) will only correct for the trend (van der Westhuizen et al., 2022).

### 2.6.2. Quantile regression forest models with weights

To analyse the contribution of HFST data, we ran several clay content prediction models.

#### 2.6.2.1. Last-only model.
A first model, referred to here as the LAST-only model, was created using solely LAST data (1,078 points; in red in Fig. 1, section 2.1.). This model used Quantile Regression Forest (QRF) using the R package ranger (version 0.15.1., Wright and Ziegler, 2017) and the case weights argument, with covariates listed in Table 2. We performed a 10-fold cross validation for model performance evaluation. Therefore, we ran the LAST-only model using 10 independent calibration and validation sets. We also used these 10 LAST-only calibration and validation sets for training (LAST + additional HFSTs) models and evaluating the performance of their predictions.

#### 2.6.2.2. Sampling additional HFSTs.
We compared several approaches to evaluate whether adding weighted HFST calibration data improves DSM prediction and to compare different sampling strategies for adding progressively more HFSTs.

In a first approach for selecting additional HFSTs, we first used the LAST-only model (with weights) to generate an uncertainty map corresponding to the extent of the 90 % prediction interval width (PIW90) by applying the difference between the 95 % upper quantile and the 5 % lower quantile (Q95 – Q5 = absolute uncertainty) of predictions. An uncertainty value of the LAST-only model prediction (PIW90) was allocated to each HFST point. Then we progressively added HFST data to the calibration datasets by selecting the points that had highest PIW90 values. We increased the number of HFST points from 1,000 to 5,000 by increments of 500, then from 5,000 to 14,000 by increments of 1,000. For each set of added HFST data, we ran 10 new models (with weights), each of them using for calibration the selected HFST points plus nine of ten LAST-only training subsets, with one of the LAST-only subsets omitted for validation. Similarly, a model (All HFST) was created by adding all HFST data (15,182 points with weights) which was trained and validated as described above.

A second approach for selecting additional HFSTs consisted of sampling HFSTs based on a sequence of uncertainty maps (PIW90 maps). This sequence of maps started with that from the LAST-only model (as per the first approach), then added to the calibration dataset the 1,000 most uncertain HFST points in the uncertainty map of the LAST-only model. A new uncertainty map was then produced based on this augmented calibration dataset (the LAST + 1000 HFST model), and that uncertainty map was used to select a further 500 (most uncertain) HFST observations to add to the calibration dataset. This process was continued to give a sequence of calibration datasets and models. In all cases, the validation set remained identical to that of the LAST-only model.

A third approach for selecting additional HFSTs used random sampling to select the HFST observations to add to LAST observations. HFST points (with weights) were randomly sampled from those HFST points that were not used in the previous model. For greater robustness, for each given number *n* of added HFST, we carried out 10 independent
samplings each feeding a (LAST + n HFST) model. We therefore produced 10 models for each n HFST.

In summary, all three approaches used weights for the LAST and HFST data, but they differed in which HFST data were used. From the LAST-only models to the All model, the calibration set contained more and more HFSTs. Note that the validation sets remained the same. The approach is summarised in Fig. 4.

### 2.6.3. Quantile regression forest models without weights

In order to evaluate if and how much damage is done by ignoring measurement uncertainty and treating LAST and HFST calibration data as equally important, we tested a fourth approach described in section 2.4. We considered a case in which the uncertainty on LAST and HFST calibration data is unknown and no weights are assigned. We also considered that the producer does not have enough paired data to recalculate a centroid of LAST values for each HFST class as we did in Section 2.3.1. Then, we used as default values of HFST measurements the centroids of the theoretical corresponding HFST classes (see Section 2.4.).

Then, we proceeded as in Section 2.6.2., except that we did not introduce weights in the modelling processes. Briefly, after having produced a first model based on LAST calibration data only (without weights), we tested the three strategies described in Section 2.6.2. and Fig. 4 to progressively add HFST calibration data without weights.

### 2.6.4. Evaluation of the prediction performances

Prediction performances were evaluated with the following standard measures.

The bias was calculated as the difference between the average of the

values predicted by the model and the average of the observed data (Equation (3)). The model is said to be unbiased when the residuals (differences between the measured values and the values predicted by the model) have a mean of zero.

$$Bias = \frac{1}{n} \sum_{i=1}^{n} (y_{i}^{*} - y_{i}) \#(3)$$

where $y_i$ is the value of the observation for point i, $y_i^*$ is the value of the prediction for point i and n is the number of validation observations.

The RMSE (Root Mean Squared Error), presented in the same units as the predicted variable, provides information on the dispersion of the predictions in relation to the observations; it is given as Equation (4).

$$RMSE = \sqrt{\sum_{i=1}^{n} \frac{(y_i - y_i^*)^2}{n}} \#(4)$$

The concordance correlation coefficient (CCC) measures the agreement between two series of measurements of the same variable (Lin, 1989), (e. g., predicted data and observed data). It is a measure of the deviation between the observed and predicted values of a variable from the bisector (1:1 line). It is defined according to Equation (5).

$$CCC = \frac{2 \times covariance_{measured}^{predicated}}{\sigma_{measured}^2 + \sigma_{predicted}^2 + (\mu_{measured} - \mu_{predicted})^2} \#(5)$$

The percentage of cross-validation observations contained in the PIW90 was calculated (prediction interval coverage probability, PICP) (Shrestha and Solomatine, 2006). The PICP assesses how well the

<figure><img src="image_6.png" /><figcaption>Fig. 4. Schematisation of the three sampling methods used to progressively increase the number of additional hand-feel soil texture (HFST) added to the laboratory measurement soil texture (LAST) for predicting topsoil clay content 10-fold by random forest modelling. In red: unchanged validation set; in blue: additional sampling based on the most uncertain LAST model predictions, using (Q95 – Q5) as a measure of uncertainty; in orange: additional sampling based on the most uncertain LAST + HFST previous model predictions, using (Q95 – Q5) as a measure of uncertainty; in green: progressive additional random sampling of HFSTs without replacement. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
prediction intervals represent prediction errors and uncertainty. If prediction intervals are providing a fair representation of uncertainty, then the PICP for a PIW90 should be close to 0.90. The PICP is given as Equation (6).

$$PICP = \frac{1}{n} \sum_{i=1}^{n} q_i \#(6)$$

where n is the number of observations in the dataset and $q_i$ is defined as following:

$$q_i = \begin{cases} 1, & siy_i \in Q^\alpha(x_i), \\ 0, & siy_i \notin Q^\alpha(x_i) \end{cases} \#(7)$$

where $Q^\alpha(x_i)$ represents the range of PI values.

### 2.6.5. Mapping

The maps (Q5, median, Q95) were exported in raster and the uncertainty map (PIW90) was calculated in R by subtracting Q5 from Q95 with the raster package (version 3.6–23, Hijmans, 2010). Final maps of predicted median clay content and PIW90 were created for all trials using ArcGIS Pro®.

## 3. Results

### 3.1. Models with weights

#### 3.1.1. Last-only model

Fig. 5 displays the covariate importance for the LAST-only model. As expected from previous research (e.g., Martelet et al., 2014; Loiseau et al., 2020; Chen et al., 2021), the four airborne gamma-ray spectrometry variables ranked among the five most important covariates, whereas two lithology classes, limestone marl (litho.marne) and sand (litho.sable), ranked third and sixth, respectively. These six top covariates were followed by three focal means of elevation, ranking from the largest to the smallest focal. The loess deposit (litho.limon) lithology class ranked tenth, followed by various relief covariates exhibiting less

importance. Note that the raw elevation (mnt) was less important than the focal means calculated on larger areas.

Fig. 6 displays the maps of predicted median values of clay and the PIW90 values. The predicted median clay content values exhibited rather contrasted and compact patterns, covering a large range of clay contents. Note however that values of clay > 45 % were very rare and predicted only for a very small number of cells that are imperceptible at the figure scale. The uncertainty map displayed a large range of values, the lowest values being located in the less clayey areas (in dark green (<7.5 %) and purple for median and uncertainty maps, respectively) whereas high uncertainty values (in red) were often located in areas with the highest clay content. Fig. 7 displays predicted vs measured values on the validation set. It is clear from Fig. 7 that values of clay over about 45 % were underestimated by the model.

#### 3.1.2. Model performance for LAST + n HFST models

Fig. 8 displays the indicators of clay content prediction performance for the three approaches to select HFST data and with progressive adding of HFST data. In each figure, the first boxplots on the left (in red) correspond to the LAST-only model (with solely the 1,078 LAST points). The last boxplot on the right (in grey) corresponds to the LAST + all HFST model. Intermediate models correspond to models with LAST and an increasing number of HFST data.

For the two strategies using uncertainty maps for HFST selection, all indicators showed that the performances of the predictions improved when adding 1,000 HFSTs points in the most uncertain areas. Note that, by construction, indicators for the LAST-only and LAST + 1000 HFST models were exactly the same in the left and middle columns of Fig. 8. This tendency of improvement of the performances remained until 1,500 HFSTs for the sampling strategy based on the absolute uncertainty of the LAST-only model, whereas the maximum improvement seemed to be reached at 1,000 HFSTs for the sampling strategy based on the absolute uncertainty of the previous model (in this case, the LAST + 1000 HFST model). From 1,500 up to about 3,000 HFSTs the performances remained rather stable for both sampling strategies. When the number of added HFSTs became higher than 5,000 to 6,000, the most striking result was the continuous and strong decrease in PICP to below 0.90, showing that both models had an over-optimistic estimate of PIW90.

For the models using progressive random sampling, changes in indicator values were very progressive. One striking behaviour is that we observed no improvement of the performance indicators until about 3,500 to 4,000 HFSTs. For HFSTs larger than about 4,000, some performance indicators suggested a regular improvement, leading to better values than those from the LAST-only model (Concordance, RMSE). The bias showed an asymptotic curve, beginning at -1 (LAST-only) and tending to zero, whereas the PICP dropped continuously from LAST-only to large numbers of HFST but at a decreasing rate. Note also that when adding the first 1,000 to 4,000 HFSTs by random sampling, the PICP dropped immediately below the theoretical optimum value of 0.90, which was not the case for the sampling based on uncertainty.

#### 3.1.3. Covariates importance for LAST + n HFST models

The importance of covariates for all the LAST-only and LAST + HFST models are displayed in Fig. S3. In supplementary material S1 we describe the main changes that were observed when increasing the number of HFST according to the sampling strategies. In summary, Th and K were always among the two most important covariates, most of the time followed by total count, or by the limestone marl lithology (litho.marne) for some models with smaller numbers of HFST data. The importance of some relief covariates increased with an increasing number of HFST data.

#### 3.1.4. Maps for LAST + n HFST models

All prediction maps and the uncertainty maps (PIW90) are displayed in supplementary materials (Figs. S2 and S3). In this section, we selected a few of these maps to illustrate their changes with the increasing

<figure><img src="image_7.png" /><figcaption>Fig. 5. Importance of covariates for the LAST-only model with weights. Covariate full names and definitions are provided in Table 2.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Predicted median clay value by the LAST-only model with weights (a) and uncertainty map (Q95-Q5) value (b).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Predicted/measured plot of the LAST-only model with weights. In red: the 1:1 line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

number of added HFST.

3.1.4.1. Maps of predicted median values of clay content. Fig. 9 displays a selection of maps of predicted median values of clay content. The LAST-only map did not capture heavy clay (clay > 45) patterns, whereas the LAST + 1000 HFST maps with sampling based on uncertainty enabled it to capture large areas of heavy clay in the north and more scattered heavy clay values in the western part of the study area. The large areas of heavy clay in the north were not completely captured by random additional sampling, whereas almost no heavy clay was predicted in the western part. In the north and north-western parts, both maps using uncertainty-based sampling exhibited some rather detailed patterns differentiating the 25–30 and the 30–35 clay% classes, whereas the maps based on random sampling were more dominated by areas of 25–30 clay%. Note also that the eastern part of the Loiret was predicted as less clayey by random sampling than by sampling based on uncertainty. Over about 6,000 HFSTs, all three maps showed rather similar

patterns except for the eastern part where topsoil was predicted to be less clayey for the random-based model. Besides, the latter model did not fully retrieve the small heavy clay hot spots that were predicted by uncertainty-based models in the western part. Note that, when using the previous map uncertainty based strategy, only very high HFST sampling densities (about 10,000 HFST data, see Fig. S4) enabled it to capture some finer-scale patterns of clayey silt in the northern area, while these patterns were captured by using 7,000 HFSTs with the other strategies and remained stable until the All HFST model.

3.1.4.2. Maps of the uncertainty of the predicted clay content values. Fig. 10 shows selected maps of the uncertainty of the predicted clay content values as estimated by the PIW90 (see Fig. S5). When compared to predicted values (see Fig. 9) it is obvious that the largest uncertainties of the LAST-only model were observed in areas with the highest clay content. Conversely, lower uncertainties were observed for the less clayey areas. When using the two HFST uncertainty-based additional sampling strategies, the most uncertain areas completely disappeared (the highest clay content uncertainty class, in red), or were strongly reduced (the class ranking second in the highly uncertain areas, in orange) as soon as 1,000 HFSTs were added. Note that the uncertainty in the most uncertain areas decreased much more progressively when using the random strategy.

If we compare the two uncertainty-based strategies, the LAST-only uncertainty-based maps lead to a smaller proportion of the most uncertain predictions than the previous uncertainty-based maps. The most certain areas were mainly located on sandy soils located north to the Loire valley and previously identified by Richer-de-Forges et al. (2023, 2008b) and Eymard et al. (2024) and were predicted with a higher accuracy when using the LAST-only uncertainty based strategy. Interestingly, up to about 3,000 additional HFSTs the two uncertainty maps were comparable, though the latter cited sandy areas remained less uncertain when using the LAST-only uncertainty-based strategy than when using the previous uncertainty-based maps. Then, when increasing the number of additional samples, the locations of the least uncertain areas progressively changed in both maps. For the LAST uncertainty-based map, the areas of smallest uncertainty (dark blue) decreased in the sandy area, whereas new areas of very low uncertainty appeared in the northern loamy part of the Loiret. Conversely, with the previous-map uncertainty-based strategy, the sandy area became less uncertain and large and compact areas of low uncertainty appeared in the northern loamy part. Note that, when using random sampling, the less uncertain areas had a tendency to remain stable or grow with additional sampling, and that some large patterns of low uncertainty were already identified when adding only 1,000 HFSTs.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Boxplots of selected indicators of prediction performance of the models predicting clay content and adding progressively more uncertain hand-feel soil texture data (HFST) to the laboratory analysis data (LAST) with weights. The horizontal dotted line corresponds to the median of the LAST-only model performance and the horizontal red line corresponds to the theoretical optimum (for Bias and PICP). Left column: sampling of HFST based on the uncertainty of the LAST-only model. Middle column: sampling of HFST based on the uncertainty of the previous model. Right column: progressive random sampling of HFST. In each column, the additional HFST observation number increases from left to right. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 3.2. Models without weights

In this section, we mainly focus on the differences observed between modelling with weights or ignoring measurement uncertainty and treating LAST and HFST calibration data as equally important. We recall that ignoring measurement uncertainty included removing weights and using the centroids of the soil texture classes according to the soil texture triangle to derive a clay content proxy value of each HFST observation.

### 3.2.1. Last-only model

Fig. 11 compares the relative importance of covariates for the LAST-only model using weights (Fig. 11a) or without weights (Fig. 11b).

The six most important covariates were the same and followed the same order, except for an inversion between uranium and total count. K

and Th were less dominant when weights were used than when they were not. This difference may be interpreted as follows. High K and Th values are positively related to clay content. Meanwhile, the MEV is positively related to clay content. Therefore, the higher the clay content observed, the lower the weight of this observation. As a consequence, the importance of the gamma-ray covariates directly related to clay content was a bit less pronounced when using weights than when considering all observations as equally accurate. However, the relative importance of the covariates remained rather comparable between Fig. 11a and Fig. 11b because the weights had a moderate effect, as LAST measurements were rather precise, even for the highest clay contents.

Fig. 12 displays the predicted median clay contents and the (Q95–Q5) clay content when using the LAST-only model with weights (top) and without weights (bottom). The maps appear almost identical.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Selected predicted maps of median clay content according to three strategies for sampling additional hand-feel soil texture (HFST) with weights.</figcaption></figure>

However, small differences are worth noting. The map of predicted clay contents without weights captured some small patterns of heavy clay that were not visible when applying weights (Fig. 12c, arrow 1). The PIW90 maps suggested that some areas with very low clay contents (sandy soils) were predicted with more certainty when using weights (Fig. 12b, arrow 3), whereas some areas with medium clay contents were predicted with more certainty when not using weights (Fig. 12d, arrow 2).

Overall, the maps with weights smoothed the variations in some areas with high clay content and delineated very low clay content sandy areas with more certainty than the maps without weights. Conversely, the map without weights showed some patterns of lower certainty in some medium clay content areas.

### 3.2.2. Model performance for LAST + n HFST models

Fig. 13 displays the statistical indicators of clay content prediction
<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Uncertainty maps of clay content predictions (Q95 – Q5) according to three strategies for sampling additional hand-feel soil texture (HFST) with weights.</figcaption></figure>

performance for the LAST + n HFST models without weights. Models from both approaches that sampled additional HFSTs based on uncertainty exhibited rather similar changes in all performance indicators. The bias reached zero when adding 1,000 HFSTs and then remained rather constant until about 8,000 HFSTs, yet behaved a bit more erratically when sampling was based on the uncertainty of the previous models. Note that random sampling led to different and regular trends in

the bias, i.e., an initial negative decrease until about 1,500 HFSTs, followed by a regular increase towards but never reaching zero. For sampling based on uncertainty, an initial increase in concordance was observed with an increasing number of HFSTs, followed by a decrease. Note that the decreasing trend was observed over 2,500 HFSTs for sampling based on the uncertainty of the previous models whereas it was observed for larger HFST numbers for sampling based on the uncertainty
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Importance of covariates for the model using solely laboratory measurements of soil texture (LAST). a: with weights; b: without weights. Covariate full names and definitions are provided in Table 2.</figcaption></figure>

of the LAST-only model. For random sampling, the concordance first decreased until about 1,500 HFSTs and remained rather stable thereafter. The RMSE displayed changes somewhat mirroring the changes in concordance for the uncertainty-based sampling, while for random sampling it exhibited a regular and asymptotic increase. Both PICPs of uncertainty-based samplings exhibited comparable and rather regular changes, showing that PIW90 became slightly more over-pessimistic when adding HFSTs. Over about 5,000 HFSTs the PICPs began to decrease, reaching the 0.90 value at about 11,000 HFSTs, and then continued to decrease showing that very numerous HFSTs led to a slightly over-optimistic PIW90.

We observed some striking differences between performance indicators of the predictions using weights (see Fig. 8) or considering all observations as equally accurate (Fig. 13). When progressively adding HFST observations all models using weights reached (at varying rates) a bias equal to zero, whereas models without weights never reached a zero bias value, even for very dense samplings. At the same time, the boxplots indicated that for highly dense sampling the bias estimates had a much larger statistical dispersion without using weights than when using weights. In other words, the results suggest that neglecting weights led to more bias than using weights, and that the confidence interval of the bias increased when neglecting weights. The concordance and the RMSE of models with weights both suggest that adding any number of HFST data resulted in equal or better performances than the LAST-only model. This was not the case for the model without weights, for which high HFST densities resulted in a marked decrease of the prediction performances. Models using weights exhibited PICPs rather close to 0.90 until adding about 4,000 HFST, then followed by much lower values suggesting that the predicted PIW90 was over-optimistic. The PIW90 of models without weights were over-pessimistic for rather low HFST densities and became over-optimistic only for the highest HFST sampling densities, though their PICP values did not become as low as for the models with weights.

### 3.2.3. Covariates importance for LAST + n HFST models

A selection of covariates importance of the models without weights is shown in Fig. S6. One striking trend was that Th and K exhibited a much more important contribution than the other covariates for rather low HFST sampling densities and when sampling the most uncertain points according to the PIW90 of previous models. The dominance of these

gamma-ray covariates was observed in all models, but the difference between sampling strategies tended to decrease with increasing densities of HFSTs.

### 3.2.4. Maps for LAST + n HFST models

Fig. 14 displays a selection of predicted median clay content maps obtained when progressively adding HFST without weights. For the two sampling strategies based on uncertainty, large areas of heavy clay were evident from the addition of the 1,000 most uncertain HFSTs. Similar areas of heavy clay were obtained for all strategies when reaching about 3,500 added HFSTs. Over about 4,000 HFSTs the maps tended to be more and more similar. They also exhibited much larger patterns of some classes, which suggests that they tended to have rather homogeneous predictions over large areas when the number of HFST measurements increased. Note that for 7,000 HFSTs the random strategy and the strategy based on the uncertainty of the previous maps exhibited very similar patterns and were also almost equivalent to the “All HFSTs” map. This was not the case for the strategy based on the uncertainty of the LAST-only map. In particular, the western part of the Loiret exhibited spatial patterns linked to geomorphology that were not detected by the LAST-only uncertainty-based strategy. Note also that these latter patterns became visible with the random strategy even for only 1,000 additional HFSTs.

Fig. 15 displays a selection of the PIW90 maps not applying weights. In all cases, very large areas of high PIW90 were observed from the addition of 1,000 HFSTs. These areas tended to grow with increasing addition of HFSTs. For the random strategy, note that some low uncertainty areas progressively appeared in the less clayey areas, especially in the very sandy soils located north of and close to the Loire valley, that were particularly identified and studied by Eymard et al. (2024). Some of the areas of lowest uncertainty also appeared progressively with random sampling in sandy soils south of the Loiret (Eymard et al., 2024) and in the western and northern part of the Loiret that are characterised by rather homogeneous loessic plateaus on which medium to clayey silt topsoil texture Luvisols developed (e.g., Nédélec, 2015; Richer-de-Forges et al., 2023).
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Maps using solely laboratory measurements of clay content (LAST-only) as calibration data; (a) predicted median clay content with applying weights to observations; (b) uncertainty map (Q95 – Q5) of clay content with applying weights to observations; (c) predicted median clay content without applying weights to observations; (d) uncertainty map (Q95 – Q5) of clay content without applying weights to observations. (1) small pattern of heavy clay in map c; (2) less uncertain patterns in map d than in map b; (3) high accuracy pattern in map b.</figcaption></figure>

# 4. Discussion

## 4.1. Limitations of the study: Looking for causes and possible improvements

We observed a rather large bias in the LAST-only model predictions. This happened even with no weights. Two possible causes for the bias may be (i) using the median instead of the mean, and (ii) using unbalanced training data. Both causes could occur. To check the impact of using the median instead of the mean we ran RF on LAST-only, the predicted mean bias was three times lower when using RF mean prediction (Bias = -0.34) than when using QRF median prediction (Bias = -1.08) suggesting that using the median instead of the mean led to systematic under-prediction of soil clay content (Fig. 16).

Another possible explanation for the bias is the clustering of LAST locations, both geographically and in feature space. Fig. 1 shows that the sampling density was highly heterogeneous within the study area, both for LAST and HFST information. This unbalanced spatial sampling is due to the use of data rescue in the studied area, which led to some areas being very densely sampled, and some other areas having much more sparse information. As we used machine learning, the geographically unbalanced sampling is not necessarily a big issue for model calibration

if the target property and covariates feature spaces are well covered. However, this unbalanced spatial sampling may have missed or under-sampled some feature space combinations. To check the feature space coverage of calibration samples, we compared the probability density functions of the most important covariates for the LAST-only, the LAST + all HFST calibration sets and the whole Loiret area (Fig. S7). Fig. S7 shows that the LAST-only calibration data sample was far from optimal as it over-sampled extreme values of some of covariates (e.g., the highest values of U, Th, K total count, and the lowest values of elevation) and under-sampled the highest values of elevation. Therefore, in these high elevation areas, the LAST-only model was not able to capture variations linked to geomorphology and slopes surrounding the plateaus. Concerning LAST + all HFST calibration data, the modes of some density plots show that rather high values U and Th were over sampled, whereas some peaked density plots clearly show that some covariates values were sampled with a narrower distribution than for LAST-only and the whole Loiret (K, total count, elevation). The most striking difference between LAST + all HFST calibration data and LAST-only (and the whole Loiret) was the over-sampling of the limestone marl by the LAST + all HFST model. Therefore, the LAST + all HFST calibration data had much higher density of high clay content locations in these areas having the highest clay content. In addition, the East and the south-east of the
Models sampling progressively
additional HFST using the absolute
uncertainty map of the LAST model -
**without weights**

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

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

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

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

Models sampling progressively
additional HFST using the absolute
uncertainty map of the previous model
- **without weights**

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

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

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

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

Models sampling progressively
additional HFST using random
sampling - **without weights**

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

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

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

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

**Fig. 13.** Boxplots of selected indicators of prediction performance of the models predicting clay content and adding progressively more uncertain hand-feel soil texture data (HFST) to the laboratory analysis data (LAST) without weights. The horizontal dotted line corresponds to the median of the LAST-only model performance and the horizontal red line corresponds to the theoretical optimum (for Bias and PICP). Left column: sampling of HFST based on the uncertainty of the LAST-only model. Middle column: sampling of HFST based on the uncertainty of the previous model. Right column: progressive random sampling of HFST. In each column, the additional HFST observation number increases from left to right. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Loiret were also under-sampled by LAST-only. For these regions, when HFSTs were progressively added to calibration data, the LAST + n HFST models with weights corrected for these biases, better captured high clay content areas located on limestone marl lithology, and gradients of clay content in the eastern and south-eastern plateaus. As high clay content areas were the most uncertain when using the LAST-only model, this explains why the uncertainty-based strategies for adding HFSTs were very efficient to capture the areas with the highest clay contents. Unbalanced spatial sampling is a common issue in DSM, especially at broad-scale where obtaining new samples may be cost-prohibitive. Therefore, though imperfect, our data corresponded to a real-life situation and to the common strategy, which is "*do the best with what you*

have". Most importantly, our results show that there was a clear need to extend the LAST-only calibration data with HFST data while accounting for HFST bigger measurement errors, which demonstrates why we need the methods presented in this article.

Brus et al. (2011) recommended collection of additional independent data and validation of soil maps by a design-based sampling strategy involving probability sampling and design-based estimation of quality measures. However, we could not afford establishing an independent validation sampling which may require a lot of data to be efficient (Lagacherie et al., 2020). Therefore, we used k-fold cross-validation, which is the most widely used technique to mimic an independent validation in broad-scale DSM (Chen et al., 2022). In addition, as stated
<figure><img src="image_27.png" /><figcaption>Fig. 14. Selected predicted maps of median clay content according to three strategies for sampling additional hand-feel soil texture (HFST) without weights.</figcaption></figure>

before, the LAST-only and the HFST observations were both spatially unbalanced. The cross-validation folds which came from LAST-only data were unbalanced too. Note that we could have reduced this effect of spatially unbalanced data by assigning weights to the test data based on local sampling density, as recommended by de Bruin et al. (2022). This could lead to an improved validation scheme, to keep in mind for further developments. However, this improvement is out of the scope of the present study.

The HFST observations were made by well-trained soil surveyors. Therefore, we obtained rather positive results. Our conclusions should not be extrapolated to conditions where soil surveyors are not skilled enough. Recently, Eymard et al. (2024) demonstrated that systematic biases in HFST measurements could occur between teams being more or less skilled to HFST determination. In this study, we eliminated this bias by assigning HFST measurements to the centroid of LAST observations. Salley et al. (2018) also showed a major influence of training on the
<figure><img src="image_28.png" /><figcaption>**Fig. 15.** Uncertainty maps of clay content predictions (Q95 – Q5) according to three strategies for sampling additional hand-feel soil texture (HFST) without weights.</figcaption></figure>

accuracy of HFST determination. In studies where this bias correction is not done or possible, we may get a strong bias effect from including HFST data and this might impair the results.

Finally, we want to point out that the quantiles produced with the QRF model represent the conditional distribution of the observed, error-infected response, not the true, error-free response. In other words, the quantiles reflect the combined variability from both the true response and the measurement errors. This is because the model was fitted on

error-infected data. While the weights in our methodology helped correct for the trend by reducing the influence of noisy data, they did not remove the measurement error. As a result, the predicted quantiles depend on the overall distribution of the error-infected response. This could be a cause of the over-optimistic prediction intervals that resulted when using a large number of uncertain data. Correcting quantiles to fully account for measurement error is a complex task and was outside the scope of this study.
<figure><img src="image_29.png" /><figcaption>**Fig. 16.** Boxplots of the bias for last-only models using bias with qrf median prediction (red) and rf mean prediction (green). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## 4.2. Importance of taking calibration data uncertainty into account

In this study, we clearly showed that allocating weights to calibration points according to their uncertainty has a major effect when HFSTs are used to increase the number of calibration data to model and predict clay content. The most striking effect was observed on the point prediction intervals (Fig. 15). Ignoring calibration data uncertainty (i.e., adding HFST data without using weights) when using HFST observations to infer clay content values led to a very large increase of the PIW90 values (see Fig. 10 with weights and Fig. 15 without weights). In other words, adding HFST data without weights confused the modelling and led to higher prediction uncertainties, which is logical. However, according to the validation statistics, the mapped uncertainty in Fig. 10 (with weights) is too optimistic (PICP around 0.8 for the LAST + all HFST model), while the mapped uncertainty in Fig. 15 (without weights) is also optimistic (PICP around 0.85 for the LAST + all HFST model), but closer to fair than the weighted approach. Surprisingly, this observation suggests a worse representation of uncertainty from the models that explicitly address the uncertainty than from the models which ignore it. One may have expected the opposite (i.e., closer to PICP = 0.9 from the model that included weights). Part of this counter-intuitive result might come from the fact that the validation data (10-fold cross validation of the LAST data) is far from optimal, i.e. it does not properly cover the feature space (Fig. S7) and it is spatially unbalanced (Fig. 1). However, it might also be a result of adding in a large number of equal values (of HFST data) to the calibration data, which led to overly smooth predictions and over-optimistic prediction intervals.

The issue of using uncertain data in prediction models has been of interest for a long time (e.g., Benjamin and Cornell, 1970; Cornell, 1972). However, taking explicitly into account the uncertainty of calibration data in DSM is an issue that has received increasing attention in the DSM community only recently (Chen et al., 2023; Hengl et al., 2018; Takoutsing and Heuvelink, 2022; van der Westhuizen et al., 2022; Wadoux et al., 2020). Many reasons explain this growing interest, among which are the growing pool of techniques (e.g., chemometrics, proximal sensing, remote sensing) to obtain proxies of soil parameters, the raising interest in global databases needing transfer functions between analytical results obtained utilising various standards, and the temporal changes of methods and accuracy of the laboratory analyses. When considering HFST harmonisation issues, one should also take into account the fact that the HFST triangles may vary widely among countries (Richer-de-Forges et al., 2022, 2008a) and that possible wrong

allocations to HFST classes should be considered (Richer-de-Forges et al., 2024). However, attempts to incorporate “soft” data obtained by field observations such as HFST into DSM of particle-size distributions remain rather scarce (Malone and Searle, 2021a, 2021b). This scarcity of the use of HFST may be attributed to the fact that HFST is often considered as conveying too much inherent unreliability (inexactness, and practically difficult to measure) as defined by van Asselt and Rotmans (2002). We show in this study that one should not ignore uncertain HFST data because they carry information. However, this can properly be done if one can estimate and incorporate the HFST measurement errors in the modelling process, as we did, using paired observations.

## 4.3. Usefulness of paired observations to infer proxies and to allocate weights of clay content for DSM

We first used existing laboratory tests to derive an uncertainty of LAST measurements and to allocate a weight to each LAST calibration measurement. Then, we used a set of paired observations of HFST and LAST to infer a proxy of clay content for each HFST, and to calculate associated variances which were used to allocate weights to HFST calibration data. These inferences enabled us to increase by 16 times the soil calibration dataset. Our results are consistent with previous studies recalculating or refining HFST classes' LAST values. Minasny et al. (2007) established the means and ranges of clay, silt, and sand contents from HFST in Australia and investigated the differences in HFST and LAST. They pointed out that this work was useful to estimate the particle size distribution and its uncertainty from HFSTs as input to pedotransfer functions, as Richer-de-Forges et al. (2024, 2022) further did in France. Malone and Searle (2021a) re-calibrated the Australian HFST centroids that were first determined by Minasny et al. (2007). Then, they used the centroids and their uncertainty to generate acceptable proxies of LAST for all soil profile descriptions including HFST in the Australian soil database. Malone and Searle (2021b) further integrated HFST with LAST into a DSM workflow to predict LAST at 90 m resolution and the six depths corresponding to the GlobalSoilMap specifications (Arrouays et al., 2014). They demonstrated a significant improvement of LAST prediction accuracy.

## 4.4. Using HFST data in DSM

In a study over Germany, Vos et al. (2016) argued that: "Field-based soil-texture estimates could replace laboratory analysis". However, they pointed out that this was possible because the HFSTs were estimated by very experienced soil scientists trained to utilise a very detailed soil texture triangle (37 classes). It is also worth noting that they had very few samples in the largest HFST classes, especially for large clay contents (>45 %) and silt contents (>50 %). Therefore, their excellent results were observed for a rather close range of soil textures, i.e., sandy and loamy soils. Conversely, our study showed that HFST data have much larger measurement error than LAST data, and that the HFST errors were pretty large and definitely not negligible. Vos et al. concluded that: "For many purposes it is sufficient to estimate the soil texture manually by feel". Though they did not detail the above-cited purposes, we strongly believe their conclusion is valid in the pedogenic and technical context of their experiment and that it supports the usefulness of HFST observations. However, for more diverse soil texture contexts, our results show that for DSM of particle-size fractions, assigning different weights to HFST classes and to LAST data is necessary.

Another big challenge is linked to the increasing need for broad-scale and global maps (Arrouays et al., 2021) which means that data from different countries and organisations are merged prior to modelling and mapping (e.g., Hengl et al., 2017; Poggio et al., 2021). One important issue then is to harmonise databases and DSM products that come from different LAST methods. The threshold values of particle sizes for determining texture components (clay, silt, sand) classes are far from homogeneous between international and national classifications (e.g.,
Arrouays et al., 2012; IUSS Working Group WRB, 2015; Minasny and McBratney, 2001; Morvan et al., 2008; NRCS-USDA, 2012; Rousseva, 1997). The use of different threshold values has led to the development of multiple transfer functions to transform particle-size distributions from one system to another (e.g., Buchan, 1989a, 1989b; Minasny et al., 2007; Nemes et al., 1999; Rousseva, 1997; Shirazi et al., 1988; Shirazi and Boersma, 1984; Takahashi et al., 2020; Yaalon, 1989). The methods used to determine LAST may differ between laboratories (e.g. pipette method versus laser-induced spectroscopy). The data may be more-or-less old and more-or-less hampered by various sources of errors. Therefore, the interest in taking into account the uncertainties when using such functions to transform LAST results to comparable values is likely to increase in the future as these uncertainties may have big consequences on the reliability of the predictions at broad-scale. This issue of merging data from heterogeneous laboratories and periods while taking account for their uncertainty is a very big challenge for the production of reliable maps at global scale and this challenge will likely be uneasy to solve if full *meta*-data about the analytical values is not available.

### 4.5. Consequences of adding many new HFST information

We tested several strategies to progressively add HFST observations to the calibration set, either (i) based on the PIW90 uncertainty map of clay content predictions based on LAST-only calibration data, or (ii) based on a sequence of PIW90 uncertainty maps, adding the HFST data from the most uncertain locations of the preceding map in the sequence of LAST + *n* HFST maps, or (iii) by random sampling of the HFST observations without replacement. We observed that sampling based on uncertainty maps was more efficient for rather low numbers of additional HFST data and that, unexpectedly, it seemed to be more efficient when based on the LAST-only PIW90 map of uncertainty than when based on the sequential uncertainty mapping approach. However, these two strategies based on uncertainty maps led to a strong decrease of the PICP when too many HFST observations were added. We assume that this result comes from an undesirable side-effect of adding too many equal values of clay content allocated to the same HFST from the cen- troids of its LAST values. These numerous equal values in the calibration set led to overfitting clay content predictions to very similar predictions when HFST data were numerous and located in similar regions of covariates and geographical space. This would also explain why, when progressively adding HFST data, the maps exhibited larger homoge- neous patterns. Note this latter trend was observed for all maps, whether using weights or not, and whether using uncertainty-based or random sampling. In other words, we believe that adding too many HFSTs led to oversimplifying the reality and to over-smoothing local patterns.

### 4.6. Uncertainty-based versus random strategies. Which uncertainty do we want to decrease? Practical considerations

In this study, we used all available HFSTs over the studied area, though we tested different strategies to add them progressively to LAST data for clay content prediction by DSM. We clearly showed that, if the number of HFSTs has to be limited, locating them in the most uncertain areas is more efficient to increase the global performance indicators than random sampling. These results are consistent with other studies rec- ommending “adaptive” sampling based on previous maps of uncertainty obtained using geostatistics (e.g., Huang et al., 2020; Marchant and Lark, 2006) or machine learning spatial predictions (e.g., Blackford et al., 2022; Stumpf et al., 2017; Szatmári et al., 2019).

When using weights and the uncertainty-based strategy, most of the global indicators of performance became better when adding the first subsets of HFST data and then remained rather stable, except for PICP which dropped below the optimal value of 0.9 when adding too many HFSTs. Note that when using weights and random sampling, much more progressive changes were observed, with a slight yet regular

improvement of bias, concordance and RMSE and a more progressive decrease of PICP. In other words, uncertainty-based strategies improved global performance by reducing the largest uncertainties by adding a rather small number of HFST data in the most uncertain areas, which were also the areas with the highest clay content. This had no significant consequence on the PICP when a small number of HFST data were used. However, as shown in Figs. 10 and 15, this improvement in the highest clay content areas increased the uncertainty in the areas with lowest clay content, where fewer HFST data were added. On the contrary, when using random sampling, the PIW90 maps show that the predictions became progressively more certain in the low clay content areas, espe- cially on sands. When not taking into account weights, the PIW90 values remained very large and the only predictions with low uncertainty were located in the least clayey areas. These considerations have two practical consequences. It may be efficient to preferentially add some HFST data in the most uncertain areas if the aim is to improve the global indicators of performance of the map and to better predict the high clay content values, but it might be detrimental to the performance of the prediction of the clay content in the low clay content areas. Conversely, if a random strategy is adopted, the predictions of the low clay content areas become more and more certain with the addition of HFST data, especially in the sandiest areas, but reducing the prediction intervals in the areas with highest clay content will require more data. Therefore, the strategy depends on the objective of the end-user. If an insurance company wants an accurate map of the high clay content areas to estimate the risk of clay shrink-swell for houses, it might choose the uncertainty-based strategy. A sedimentologist interested in aeolian erosion hazards might opt for the random strategy.

Overall, both should not forget that adding too many HFSTs may smooth the prediction maps and give a false impression of accuracy, as many equal calibration values may lead to an underestimation of un- certainty as represented in the PIW90 maps.

## 5. Conclusion

In this study, hand-feel soil texture field observations (HFST data) were progressively added to laboratory analyses of soil texture (LASTs) to increase the number of calibration data to predict topsoil clay content using machine learning digital soil mapping. We tested these progressive additions of HFSTs by using weights, taking into account the uncertainty of measurements, and compared this with a case where measurement errors are ignored and all measurements are equally important. We also tested sampling strategies based on focusing on the most uncertain areas or adding HFST randomly. Our main findings are the following:

* Considering all measurements as equally important led to a large increase in the uncertainty of clay content predictions.

* Adding progressively HFSTs in the most uncertain areas enabled the resulting model to map with better accuracy and improved indicators of prediction performance, though it led to a small decrease in ac- curacy of the least uncertain areas.

* Conversely, when a random strategy was adopted for adding HFSTs to the calibration data, the predictions of the least uncertain areas became sharper with the addition of HFSTs, but reducing the pre- diction intervals in the most uncertain areas required more data.

* Adding a very large number of HFSTs led to an optimistic estimation of the 90 % prediction interval and to larger homogeneous patterns, which smoothed the mapped spatial variation of clay content.

* Adding weighted HFSTs substantially improved DSM predictions, but the number of HFSTs and the strategy to add them had to be carefully chosen so as not to over-smooth the maps nor under- estimate the uncertainty of predictions.

### CRediT authorship contribution statement

Anne C. Richer-de-Forges: Writing – review & editing, Writing –
original draft, Visualization, Validation, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **Songchao Chen:** Writing – review & editing, Validation, Software, Resources, Methodology. **Gerard M. Heuvelink:** Writing – review & editing, Validation, Supervision, Methodology, Formal analysis, Conceptualization. **Stephan van der Westhuizen:** Writing – review & editing. **Thomas G. Orton:** Writing – review & editing, Validation, Methodology. **Hocine Bourennane:** Writing – review & editing, Supervision, Project administration. **Dominique Arrouays:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Investigation, 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.

Acknowledgments

Most of the soil data were gathered during soil surveys conducted under the programme "Inventaire, Conservation et Gestion des Sols" partly funded by the French Ministry in charge of Agriculture and other programmes (Carte des sols de la Région Centre) and fundings (INRAE, Régions, Départements, Agences de bassin, Chambres d'Agriculture...). We thank all these organisations and all the soil surveyors for sampling, describing soils, and providing data to the national DoneSol database. We are grateful to the "Groupement d'Intérêt Scientifique Sol" for his support of the national French Soil Information System.

The airborne gamma-spectrometric data used in this study were financed by the Région Centre Val-de-Loire, the BRGM, France, and FEDER European funds.

A.C. Richer-de-Forges is the coordinator of the Domaine d'Expertise Scientifique Theia "Cartographie Numérique des Sols" supported by the CNES, France. D. Arrouays and H. Bourennane are members of this Domaine d'Expertise Scientifique Theia "Cartographie Numérique des Sols".

Appendix A. Supplementary data

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

Data availability

The authors do not have permission to share data.

References

Afnor, 2003. Qualité du sol - Détermination de la distribution granulométrique des particules du sol - Méthode à la pipette. NF X 31-107. Association Française de Normalisation.

Armas, D., Guevara, M., Bezares, F., Vargas, R., Durante, P., Osorio, V., Jimenez, W., Oyonarte, C., 2023. Harmonized Soil Database of Ecuador (HESD): data from 2009 to 2015. Earth Syst. Sci. Data 15, 431–445. https://doi.org/10.5194/essd-15-431-2023.

Arrouays, D., Leenaars, J.G.B., Richer-de-Forges, A.C., Adhikari, K., Ballabio, C., Greve, M., Grundy, M., Guerrero, E., Hempel, J., Hengl, T., Heuvelink, G., Batjes, N., Carvalho, E., Hartemink, A., Hewitt, A., Hong, S.-Y., Krasilnikov, P., Lagacherie, P., Lelyk, G., Libohova, Z., Lilly, A., McBratney, A., McKenzie, N., Vasquez, G.M., Mulder, V.L., Minasny, B., Montanarella, L., Odeh, I., Padarian, J., Poggio, L., Roudier, P., Saby, N., Savin, I., Searle, R., Solbovoy, V., Thompson, J., Smith, S., Sulaeman, Y., Vintila, R., Rossel, R.V., Wilson, P., Zhang, G.-L., Swerts, M., Oorts, K., Karklins, A., Feng, L., Ibelles Navarro, A.R., Levin, A., Laktionova, T., Dell'Acqua, M., Suvannang, N., Ruam, W., Prasad, J., Patil, N., Husnjak, S., Pásztor, L., Okky, J., Hallett, S., Keay, C., Farewell, T., Lilja, H., Juilleret, J., Marx, S., Takata, Y., Kazuyuki, Y., Mansuy, N., Panagos, P., Van Liedekerke, M., Skalsky, R., Sobocka, J., Kobza, J., Eftekhari, K., Alavipanah, S.K., Moussadek, R., Badraoui, M., Da Silva, M., Paterson, G., Gonçalves, M. da C., Theocharopoulos, S., Yemefack, M., Tedou, S., Vrscaj, B., Grob, U., Kozák, J., Boruvka, L., Dobos, E., Taboada, M., Moretti, L.,

Rodriguez, D., 2017. Soil legacy data rescue via GlobalSoilMap and other international and national initiatives. Geo. Res. J. 14, 1–19. Doi: 10.1016/j.grj.2017.06.001.

Arrouays, D., Marchant, B.P., Saby, N.P.A., Meersmans, J., Orton, T.G., Martin, M.P., Bellamy, P.H., Lark, R.M., Kibblewhite, M., 2012. Generic issues on broad-scale soil monitoring schemes: A review. Pedosphere 22, 456–469.

Arrouays, D., McKenzie, N., Hempel, J., Richer-de-Forges, A.C., McBratney, A.B. (Eds.), 2014. GlobalSoilMap: Basis of the global spatial soil information system, 1st Edition. ed. CRC Press, London. Doi: 10.1201/b16500.

Arrouays, D., Mulder, V.L., Richer-de-Forges, A.C., 2021. Soil mapping, digital soil mapping and soil monitoring over large areas and the dimensions of soil security – A review. Soil Secur. 5, 100018. https://doi.org/10.1016/j.soisec.2021.100018.

Arrouays, D., Richer-de-Forges, A.C., Héliès, F., Mulder, V.L., Saby, N.P.A., Chen, S., Martin, M.P., Román Dobarco, M., Follain, S., Jolivet, C., Laroche, B., Loiseau, T., Cousin, I., Lacoste, M., Ranjard, L., Toutain, B., Le Bas, C., Eglin, T., Bardy, M., Antoni, V., Meersmans, J., Ratié, C., Bispo, A., 2020. Impacts of national scale digital soil mapping programs in France. Geoderma Reg. 23, e00337. https://doi.org/10.1016/j.geodrs.2020.e00337.

Baize, D., 2018. Guide des analyses en pédologie, 3e éd. revue et augmentée. ed. Savoir-faire. Éditions Quae, Versailles.

Batjes, N.H., Ribeiro, E., van Oostrum, A., Leenaars, J., Hengl, T., Mendes de Jesus, J., 2017. WoSIS: providing standardised soil profile data for the world. Earth Syst. Sci. Data 9, 1–14. https://doi.org/10.5194/essd-9-1-2017.

Benjamin, J.R., Cornell, C.A., 1970. Probability, statistics, and decision for civil engineers. McGraw-Hill, New York.

Beven, K.J., Kirkby, M.J., 1979. A physically based, variable contributing area model of basin hydrology / Un modèle à base physique de zone d'appel variable de l'hydrologie du bassin versant. Hydrol. Sci. Bull. 24, 43–69. https://doi.org/10.1080/02626667909491834.

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

Breiman, L., 2001. Random forests. Mach. Learn. 45, 5–32. https://doi.org/10.1023/A:1010933404324.

Brus, D.J., Kempen, B., Heuvelink, G.B.M., 2011. Sampling for validation of digital soil maps. Eur. J. Soil Sci. 62, 394–407. https://doi.org/10.1111/j.1365-2389.2011.01364.x.

Buchan, G.D., 1989a. Comments on "A Unifying Quantitative Analysis of Soil Texture: Improvement of Precision and Extension of Scale.". Soil Sci. Soc. Am. J. 53, 594. https://doi.org/10.2136/sssjaa1989.03615995005300020052x.

Buchan, G.D., 1989b. Applicability of the Simple Lognormal Model to Particle-size Distribution in Soils. Soil Sci. 147, 155–161. https://doi.org/10.1097/00010694-198903000-00001.

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. https://doi.org/10.1016/j.geoderma.2021.115567.

Chen, S., Richer-de-Forges, A.C., Mulder, V.-L., Martelec, G., Loiseau, T., Lehmann, S., Arrouays, D., 2021. Digital mapping of the soil thickness of loess deposits over a calcareous bedrock in central France. Catena 198, 105062. https://doi.org/10.1016/j.catena.2020.105062.

Chen, S., Saby, N.P.A., Martin, M.P., Barthès, B.G., Gomez, C., Shi, Z., Arrouays, D., 2023. Integrating additional spectroscopically inferred soil data improves the accuracy of digital soil mapping. Geoderma 433, 116467. https://doi.org/10.1016/j.geoderma.2023.116467.

Cooper, M., Mendes, L.M.S., Silva, W.L.C., Sparovek, G., 2005. A National Soil Profile Database for Brazil Available to International Scientists. Soil Sci. Soc. Am. J. 69, 649–652. https://doi.org/10.2136/sssjaa2004.0140.

Cornell, C.A., 1972. First order analysis of model and parameter uncertainty, in: Proceedings. Presented at the International Symposium on Uncertainties in Hydrologic and Water Resource Systems, Tucson, Arizona, pp. 1245–1272.

de Bruin, S., Brus, D.J., Heuvelink, G.B.M., Van Ebbenhorst Tengbergen, T., Wadoux, A.-M.-J.-C., 2022. Dealing with clustered samples for assessing map accuracy by cross-validation. Eco. Inform. 69, 101665. https://doi.org/10.1016/j.ecoinf.2022.101665.

de Sousa, L.M., Poggio, L., Batjes, N.H., Heuvelink, G.B.M., Kempen, B., Riberio, E., Rossiter, D., 2020. SoilGrids 2.0: producing quality-assessed soil information for the globe (preprint). Soils and the natural environment. Doi: 10.5194/soil-2020-65.

Eymard, A., Richer-de-Forges, A.C., Martelec, G., Tissoux, H., Bialkowski, A., Dalmasso, M., Chretien, F., Belletier, D., Ledemé, G., Laloua, D., Josière, O., Commagnac, L., Bourennane, H., Arrouays, D., 2024. Exploring the untapped potential of hand-feel soil texture data for enhancing digital soil mapping: Revealing hidden spatial patterns from field observations. Geoderma 441, 116769. https://doi.org/10.1016/j.geoderma.2023.116769.

Grolleau, E., Bargeot, L., Chafchafi, A., Hardy, R., Doux, J., Beaudou, A., Le Martret, H., Lacassin, J.-C., Fort, J.-L., Falipou, P., Arrouays, D., 2004. Le système d'information national sur les sols : DoneSol et les outils associés. Etude et Gestion Des Sols 11, 255–269 [In French].

Grunwald, S., 2009. Multi-criteria characterization of recent digital soil mapping and modeling approaches. Geoderma 152, 195–207. https://doi.org/10.1016/j.geoderma.2009.06.003.

Guevara, M., Arroyo, C., Brunsell, N., Cruz, C.O., Domke, G., Equihua, J., Etchevers, J., Hayes, D., Hengl, T., Ibelles, A., Johnson, K., Jong, B., Libohova, Z., Llamas, R., Nave, L., Ornelas, J.L., Paz, F., Ressl, R., Schwartz, A., Victoria, A., Wills, S., Vargas, R., 2020. Soil Organic Carbon Across Mexico and the Conterminous United
States (1991-2010). Global Biogeochem. Cycles 34. https://doi.org/10.1029/2019GB006219.

Guevara, M., Olmedo, G.F., Stell, E., Yigini, Y., Aguilar Duarte, Y., Arellano Hernández, C., Arévalo, G.E., Arroyo-Cruz, C.E., Bolivar, A., Bunning, S., Bustamante Canas, N., Cruz-Gaistardo, C.O., Davila, F., Dell Acqua, M., Encina, A., Figueredo Tacona, H., Fontes, F., Hernández Herrera, J.A., Ibelles Navarro, A.R., Loayza, V., Manueles, A.M., Mendoza Jara, F., Olivera, C., Osorio Hermosilla, R., Pereira, G., Prieto, P., Alexis Ramos, I., Rey Brina, J.C., Rivera, R., Rodríguez-Rodríguez, J., Roopnarine, R., Rosales Ibarra, A., Rosales Riveiro, K.A., Schulz, G.A., Spence, A., Vasques, G.M., Vargas, R.R., Vargas, R., 2018. No Silver Bullet for Digital Soil Mapping: Country-specific Soil Organic Carbon Estimates across Latin America (preprint). Soil Method. https://doi.org/10.5194/soil-2017-40.

Hartemink, A.E., McBratney, A.B., Mendonça-Santos, M. d. L., 2008. Digital Soil Mapping with Limited Data. Springer.

Hengl, T., Mendes de Jesus, J., Heuvelink, G.B.M., Ruiperez Gonzalez, M., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M.N., Geng, X., Bauer-Marschallinger, B., Guevara, M.A., Vargas, R., MacMillan, R.A., Batjes, N.H., Leenaars, J.G.B., Ribeiro, E., Wheeler, I., Mantel, S., Kempen, B., 2017. SoilGrids250m: Global gridded soil information based on machine learning. PLoS One 12, e0169748. https://doi.org/10.1371/journal.pone.0169748.

Hengl, T., Nussbaum, M., Wright, M.N., Heuvelink, G.B.M., Gräler, B., 2018. Random forest as a generic framework for predictive modeling of spatial and spatio-temporal variables. PeerJ 6, e5518.

Hijmans, R.J., 2010. Package raster: Geographic Data Analysis and Modeling. Doi: 10.32614/CRAN package.raster.

Huang, J., McBratney, A.B., Minasny, B., Malone, B., 2020. Evaluating an adaptive sampling algorithm to assist soil survey in New South Wales, Australia. Geoderma Regl. 21, e00284. https://doi.org/10.1016/j.geodrs.2020.e00284.

Hugelius, G., Tarnocai, C., Broll, G., Canadell, J.G., Kuhry, P., Swanson, D.K., 2013. The Northern Circumpolar Soil Carbon Database: spatially distributed datasets of soil coverage and soil carbon storage in the northern permafrost regions. Earth Syst. Data 5, 3-13. https://doi.org/10.5194/essd-5-3-2013.

IAEA, 2003. Guidelines for radioelement mapping using gamma ray spectrometry data (IAEA-TECDOC-1363.).

IAEA, 1991. Airborne gamma ray spectrometer surveying (Technical Reports No. Series No. 323.). IAEA.

IGN, 2014. BD ALTI®. https://geoservices.ign.fr/bdalti.

InfoSol, 2021. Dictionnaire de données - DoneSol version 3.15.

IUSS Working Group WRB, 2015. World Reference Base for Soil Resources 2014, update 2015 International soil classification system for naming soils and creating legends for soil maps.

Jamagne, M., 1967. Bases et techniques d'une cartographie des sols. Annales agronomiques, Numéro hors série [in French].

Lacquier, F., Prognon, F., Tourlière, B., Martelet, G., Deparis, J., Tissoux, H., Reninger, P.-A., Perrin, J., 2015. Méthodologie de cartographie du récollet à partir de données radiométriques aéroportées en région Centre - Etablissement de cartes lithologiques (Rapport BRGM Orléans en BRGM/RP-64932-FR). BRGM, Orléans.

Lagacherie, P., Arrouays, D., Bourennane, H., Gomez, C., Nkuba-Kasanda, L., 2020. Analysing the impact of soil spatial sampling on the performances of Digital Soil Mapping models and their evaluation: A numerical experiment on Quantile Random Forest using clay contents obtained from Vis-NIR-SWIR hyperspectral imagery. Geoderma 375, 114503. https://doi.org/10.1016/j.geoderma.2020.114503.

Leenaars, J.G.B., n.d. Africa Soil Profiles Database Version 1.2 166.

Levi, M.R., 2017. Modified Centroid for Estimating Sand, Silt, and Clay from Soil Texture Class. Soil Sci. Soc. Am. J. 81, 578-588. https://doi.org/10.2136/ssaaj2016.09.0301.

Liard, M., Tissoux, H., Deschamps, S., 2017. Les alluvions anciennes de la Loire en orléanais (France, Loiret), une relecture à l'aune de travaux d'archéologie préventive et d'un programme de datations ESR. Quaternaire 105-128. https://doi.org/10.4000/quaternaire.7914.

Lin, L.-I.-K., 1989. A Concordance Correlation Coefficient to Evaluate Reproducibility. Biometrics 45, 255. https://doi.org/10.2307/2532051.

Loiseau, T., Arrouays, D., Richer-de-Forges, A.C., Lagacherie, P., Ducommun, C., Minasny, B., 2021. Density of soil observations in digital soil mapping: A study in the Mayenne region, France. Geoderma Regional 24, e00358. https://doi.org/10.1016/j.geodrs.2021.e00358.

Loiseau, T., Richer-de-Forges, A.C., Martelet, G., Bialkowski, A., Nehlig, P., Arrouays, D., 2020. Could airborne gamma spectrometric data replace lithological maps as co- variates for digital soil mapping of topsoil particle-size distribution? A case study in Western France. Geoderma Reg. 22, e00295. https://doi.org/10.1016/j.geodrs.2020.e00295.

Malone, B., Searle, R., 2021a. Updating the Australian digital soil texture mapping (Part 1"): re-calibration of field soil texture class centroids and description of a field soil texture conversion algorithm. Soil Res. 59, 419-434. https://doi.org/10.1071/SR20283.

Malone, B., Searle, R., 2021b. Updating the Australian digital soil texture mapping (Part 2"): spatial modelling of merged field and lab measurements. Soil Res. 59, 435-451. https://doi.org/10.1071/SR20284.

Malone, B.P., Odgers, N.P., Stockmann, U., Minasny, B., McBratney, Alex. B., 2018. Digital Mapping of Soil Classes and Continuous Soil Properties, in: McBratney, Alex. B., Minasny, B., Stockmann, U. (Eds.), Pedometrics, Progress in Soil Science. Springer International Publishing, Cham, pp. 373-413. Doi: 10.1007/978-3-319-63439-5_12.

Marchant, B.P., Lark, R.M., 2006. Adaptive sampling and reconnaissance surveys for geostatistical mapping of the soil. Europ. J. Soil Sci. 57, 831-845. https://doi.org/10.1111/j.1365-2389.2005.00774.x.

Martelet, G., Nehlig, P., Arrouays, D., Messner, F., Tourlière, B., Laroche, B., Deparis, J., Saby, N.P.A., 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, pp. 401-408.

McBratney, A.B., Mendonça Santos, M.L., Minasny, B., 2003. On digital soil mapping. Geoderma 117, 3-52. https://doi.org/10.1016/S0016-7061(03)00223-4.

Medeiros, B.M., Rossi, L.S., Caten, A.T., Pereira, G.E., Silva, E.B.D., Daboit, K.T.U., 2024. Soil legacy data: An opportunity for digital soil mapping. Revista Brasileira De Ciência do Solo 48, e0230130. https://doi.org/10.36783/18069657rbcs20230130.

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

Minasny, B., McBratney, A.B., 2001. The Australian soil texture boomerang: a comparison of the Australian and USDA/FAO soil particle-size classification systems. Soil Res. 39, 1443. https://doi.org/10.1071/SR00065.

Minasny, B., McBratney, A.B., Field, D.J., Tranter, G., McKenzie, N.J., Brough, D.M., 2007. Relationships between field texture and particle-size distribution in Australia and their implications. Soil Res. 45, 428. https://doi.org/10.1071/SR07051.

Moes, J., Shangguan, W., Petzold, R., Minasny, B., Rosca, B., Jelinski, N., Zelazny, W., Marcondes Silva Souza, R., Safanelli, J.L., ten Caten, A., 2018. R package "soiltexture": Functions for Soil Texture Plot, Classification and Transformation.

Morvan, X., Saby, N.P.A., Arrouays, D., Le Bas, C., Jones, R.J.A., Verheijen, F.G.A., Bellamy, P.H., Stephens, M., Kibblewhite, M.G., 2008. Soil monitoring in Europe: a review of existing systems and requirements for harmonisation. Sci. Total Environ. 391, 1-12. https://doi.org/10.1016/j.scitotenv.2007.10.046.

Nédélec, H., 2015. Carte des sols de Pithiviers au 1/50 000. Chambre d'Agriculture du Loiret. [In French].

Nemes, A., Wösten, J.H.M., Lilly, A., Oude Voshaar, J.H., 1999. Evaluation of different procedures to interpolate particle-size distributions to achieve compatibility within soil databases. Geoderma 90, 187-202. https://doi.org/10.1016/S0016-7061(99)00142-8.

NRCS-USDA, 2012. Field book for describing and sampling soils, Version 3.0. National Soil Survey Center. National resources Conservation service. U.S. Department of Agriculture. Lincoln, NE, USA.

Olsen, L.R., 2020. Package rearrr: Rearranging Data. Doi: 10.32614/CRAN package.rearr.

Orgiazzi, A., Ballabio, C., Panagos, P., Jones, A., Fernández-Ugalde, O., 2018. LUCAS Soil, the largest expandable soil dataset for Europe: a review. Eur. J. Soil Sci. 69, 140-153. https://doi.org/10.1111/ejss.12499.

Panagos, P., Van Liedekerke, M., Jones, A., Montanarella, L., 2012. European Soil Data Centre: Response to European policy support and public data requirements. Land Use Policy 29, 329-338. https://doi.org/10.1016/j.landusepol.2011.07.003.

Pfeiffer, M., Padarian, J., Osorio, R., Tourmante, N., Olmedo, G.F., Guevara, M., Aburto, F., Albornoz, F., Antilén, M., Araya, E., Arellano, E., Barret, M., Barrera, J., Boeckx, P., Briceno, M., Bunning, S., Cabrol, L., Casanova, M., Cornejo, P., Corradini, F., Cuarqueo, G., Doetterl, S., Duran, P., Escudey, M., Espinosa, A., Francke, S., Fuentes, J.P., Fuentes, M., Gajardo, G., García, R., Gallaud, M., Galleguillos, M., Gomez, A., Hidalgo, M., Ivelic-Sáez, J., Mashalaba, L., Matus, F., Meza, F., Mora, M. de la L., Mora, J., Muñoz, C., Norambuena, P., Olivera, C., Ovalle, C., Panichini, M., Pauchard, A., Pérez-Quezada, J.F., Radic, S., Ramirez, J., Riveras, N., Ruiz, G., Salazar, O., Salgado, I., Seguel, O., Sepúlveda, M., Sierra, C., Tapia, Y., Tapia, F., Toledo, B., Torrico, J.M., Valle, S., Vargas, R., Wolff, M., Zagal, E., 2020. CHLSOC: The Chilean Soil Organic Carbon database, a multi-institutional collaborative effort. Earth System Science Data 12, 457-468. Doi: 10.5194/essd-12-457-2020.

Poggio, L., de Sousa, L.M., Batjes, N.H., Heuvelink, G.B.M., Kempen, B., Ribeiro, E., Rossiter, D., 2021. SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty. Soil 7, 217-240. https://doi.org/10.5194/soil-7-217-2021.

Rasaei, Z., Rossiter, D.G., Farshad, A., 2020. Rescue and renewal of legacy soil resource inventories in Iran as an input to digital soil mapping. Geoderma Reg. 21, e00262. https://doi.org/10.1016/j.geodrs.2020.e00262.

Richer-de-Forges, A.C., Arrouays, D., Chen, S., Román Dobarco, M., Libohova, Z., Roudier, P., Minasny, B., Bourennane, H., 2022. Hand-feel soil texture and particle-size distribution in central France. Relationships and Implications. CATENA 106155. https://doi.org/10.1016/j.catena.2022.106155.

Richer-de-Forges, A.C., Arrouays, D., Feller, C., Jamagne, M., 2008a. Perdus dans le triangle des textures. Etude et Gestion Des Sols 15, 97-111.

Richer-de-Forges, A.C., Arrouays, D., Poggio, L., Chen, S., Lacoste, M., Minasny, B., Libohova, Z., Roudier, P., Mulder, V.L., Nédélec, H., Martelet, G., Lemercier, B., Lagacherie, P., Bourennane, H., 2023. Hand-feel soil texture observations to evaluate the accuracy of digital soil maps for local prediction of soil particle size distribution: A case study in Central France. Pedosphere 33, 731-743. https://doi.org/10.1016/j.pedsph.2022.07.009.

Richer-de-Forges, A.C., Chen, S., Arrouays, D., Bourennane, H., Minasny, B., 2024. Hand-feel soil texture classes and particle-size distribution as predictors of soil water content at field capacity. Further insights into the sources of uncertainty. Catena 244, 108268. https://doi.org/10.1016/j.catena.2024.108268.

Richer-de-Forges, A.C., Renaux, B., Verbèque, B., Souc伙erianadin, N., Eimberck, M., 2008b. Référentiel Régional Pédologique de la région Centre : carte des pépaysages du Loiret à 1/250 000 [In french].

Ritchey, E.L., McGrath, J.M., Gehring, D., 2015. Determining Soil Texture by Feel. Agriculture and Natural Resources Publications 4.

Rousseva, S.S., 1997. Data transformations between soil texture schemes. Eur. J. Soil Sci. 48, 749-758. https://doi.org/10.1111/j.1365-2389.1997.tb00574.x.
Salley, S.W., Herrick, J.E., Holmes, C.V., Karl, J.W., Levi, M.R., McCord, S.E., van der Waal, C., Van Zee, J.W., 2018. A Comparison of Soil Texture-by-Feel Estimates: implications for the Citizen Soil Scientist. Soil Sci. Soc. Am. J. 82, 1526–1537. https://doi.org/10.2136/ssaj2018.04.0137.

Shirazi, M.A., Boersma, L., 1984. A unifying quantitative analysis of soil texture. Soil Sci. Soc. Am. J. 48, 142–147. https://doi.org/10.2136/ssaj1984.03615995004800010026x.

Shirazi, M.A., Boersma, L., Hart, J.W., 1988. A unifying quantitative analysis of soil texture: improvement of precision and extension of scale. Soil Sci. Soc. Amer. J. 52, 181–190. https://doi.org/10.2136/ssaj1988.03615995005200010032x.

Shrestha, D.L., Solomatine, D.P., 2006. Machine learning approaches for estimation of prediction interval for the model output. Neural Netw. 19, 225–235. https://doi.org/10.1016/j.neunet.2006.01.012.

Simões, M.G., 2015. Democrati-zação da informação de solos do brasil: geortal e anco de dados de solos com acesso via web 32, 16.

Somarathna, P.D.S.N., Minasny, B., Malone, B.P., 2017. More Data or a Better Model? Figuring Out What Matters Most for the Spatial Prediction of Soil Carbon. Soil Sci. Soc. Am. J. 81, 1413–1426. https://doi.org/10.2136/ssaj2016.11.0376.

Stumpf, F., Schmidt, K., Goebes, P., Behrens, T., Schönbrodt-Stitt, S., Wadoux, A., Xiang, W., Scholten, T., 2017. Uncertainty-guided sampling to improve digital soil maps. Catena 153, 30–38. https://doi.org/10.1016/j.catena.2017.01.033.

Sulaeman, Y., Minasny, B., McBratney, A.B., Sarwani, M., Sutandi, A., 2013. Harmonizing legacy soil data for digital soil mapping in Indonesia. Geoderma 192, 77–85. https://doi.org/10.1016/j.geoderma.2012.08.005.

Szatmári, G., László, P., Takács, K., Szabó, J., Bakacsi, Z., Koós, S., Pásztor, L., 2019. Optimization of second-phase sampling for multivariate soil mapping purposes: case study from a wine region, Hungary. Geoderma 352, 373–384. https://doi.org/10.1016/j.geoderma.2018.02.030.

Takahashi, T., Nakano, K., Nira, R., Kumagai, E., Nishida, M., Namikawa, M., 2020. Conversion of soil particle size distribution and texture classification from ISSS system to FAO/USDA system in Japanese paddy soils. Soil Sci. Plant Nutr. 66, 407–414. https://doi.org/10.1080/00380768.2020.1763143.

Takoutsing, B., Heuvelink, G.B.M., 2022. Comparing the prediction performance, uncertainty quantification and extrapolation potential of regression kriging and

random forest while accounting for soil measurement errors. Geoderma 428, 116192. https://doi.org/10.1016/j.geoderma.2022.116192.

Terraquest LTD, 2009. Airborne horizontal gradient magnetic and radiometric survey, Centre region, France. (Operation report No. Operation report B-268, BRGM.). BRGM, France.

Thien, S., 1979. Determining Soil Texture by the “Feel Method.” Jour. Agronomy education.

Tissoux, H., Prognon, F., Martelet, G., Tourlière, B., Despriée, J., Liard, M., Lacquement, F., 2017. Contribution of airborne gamma spectrometry to both characterisation and mapping of fluvial silico-clastic deposits in Loire valley and Sologne (Centre, France). Quaternaire 28, 87–103.

USDA-NRCS, 2014. Keys to soil taxonomy. 12th edition. United States Department of Agriculture, National Resources Conservation Service, Lincoln, NE, USA.

Van Asselt, M.B.A., Rotmans, J., 2002. Uncertainty in integrated assessment modelling. Clim. Change 54, 75–105. https://doi.org/10.1023/A:1015783803445.

van der Westhuizen, S., Heuvelink, G.B.M., Hofmeyr, D.P., Poggio, L., 2022. Measurement error-filtered machine learning in digital soil mapping. Spatial Stat. 47, 100572. https://doi.org/10.1016/j.spasta.2021.100572.

Vos, C., Don, A., Prietz, R., Heidkamp, A., Freibauer, A., 2016. Field-based soil-texture estimates could replace laboratory analysis. Geoderma 267, 215–219. https://doi.org/10.1016/j.geoderma.2015.12.022.

Wadoux, A.-M.-J.-C., Minasny, B., McBratney, A.B., 2020. Machine learning for digital soil mapping: applications, challenges and suggested solutions. Earth Sci. Rev. 210, 103359. https://doi.org/10.1016/j.earscirev.2020.103359.

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

Wright, M.N., Ziegler, A., 2017. ranger: a fast implementation of random forests for high dimensional data in C++ and R. J. Stat. Softw. 77. https://doi.org/10.18637/jss.v077.i01.

Yaalon, D.H., 1989. Comments on “A Unifying Quantitative Analysis of Soil Texture.”. Soil Sci. Soc. Am. J. 53, 595. https://doi.org/10.2136/ssaj1989.03615995005300020053x.