Philippe Lagacherie

D Arrouays

H Bourennane

Cecile Gomez

L Nkuba-Kasanda

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

Keywords: Uncertainty, sampling methods, spatial distribution indicators, Quantile Random Forest

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Digital Soil Mapping has now emerged as a credible solution for providing soil data to decision-makers acting at global or local scales. Following the GlobalSoilMap specifications [START_REF] Arrouays | The GlobalSoilMap project specifications[END_REF], a significant number of countries across the world are now covered by a high-resolution (90-m) grid documented with estimates of soil properties and their associated prediction uncertainties. To obtain this spatial soil information, Digital Soil Mapping models relating the targeted soil properties with exhaustively available covariates were built by running learning algorithms onto spatial sampling of available legacy soil profiles with measured soil properties. Numerous different algorithms have been proposed, mainly depending on the availability and type of soil data (e.g. point data (profiles, augerings), soil maps with various scales and detailed legends) and on the availability of covariates [START_REF] Minasny | Methodologies for global soil mapping[END_REF].

It has been largely acknowledged that the main limitation of such digital models is the spatial sampling of the legacy measured soil profiles. The average spacing of such soil spatial sampling used in most operational DSM applications have been very large, e.g. 28 km [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF] and 5.5 km [START_REF] Mulder | GlobalSoilMap France: High-resolution spatial modelling the soils of France up to two meter depth[END_REF]. Furthermore, these soils' spatial samplings have been constituted by the soil profiles of the different existing soil surveys within the mapped area, which means that their locations were not specifically selected for a DSM application and are often clustered at a few sub-areas of the study region, with large subareas without any soil information. Consequently, the performances of soil predictions are often severely limited, especially for soil properties whose pattern of variation is largely below the soil profiles' spacing (Vaysse and Lagacherie, 2016;[START_REF] Gomez | Importance of the spatial extent for using soil properties estimated by laboratory VNIR/SWIR spectroscopy: Examples of the clay and calcium carbonate content[END_REF].

One way to increase DSM performance is therefore to increase the average densities and improve the spatial distribution of the measured soil profiles (Voltz et al., submitted).

Because of their high cost of acquisition, such new locations should be selected with care.

For that, the prerequisites are: i) to accurately evaluate the uncertainty of the current soil property maps made from the unevenly distributed legacy measured soil profiles, ii) to determine the expected gain of an increase in the average spatial sampling densities and iii) to identify the spatial sampling characteristics that should be considered in the sampling strategies. These three prerequisites will be examined in this paper and are successively evoked in the following.

The uncertainty evaluation of the soil property maps is an intrinsic part of a DSM application.

This was theorized very early [START_REF] Mcbratney | On digital soil mapping[END_REF] and translated into specifications in the GlobalSoilMap project [START_REF] Arrouays | The GlobalSoilMap project specifications[END_REF]. Although there exist guidelines for the assessment of uncertainty of Digital Soil Maps [START_REF] Heuvelink | Uncertainty quantification of GlobalSoilMap products[END_REF], their practical application is not straightforward, which makes the uncertainty evaluation itself uncertain. The spatial sampling from which uncertainty is evaluated is a crucial point: the probability sampling that is advocated to obtain an unbiased estimate of uncertainty [START_REF] Brus | Sampling for validation of digital soil maps[END_REF]) cannot be obtained with the already fixed legacy soil data locations. Furthermore, uncertain estimations of the uncertainty of soil maps can be non-negligible even when using independent probability sampling [START_REF] Kempen | Three-dimensional mapping of soil organic matter content using soil type-specific depth functions[END_REF], Lagacherie et al., 2019). Lastly, it is uncertain that the set of locations used as evaluation data well describes the real pattern of the soil cover that is to be mapped, especially the short-range variations. For all of these reasons, determining the uncertainty of a DSM map is still an open question.

Although there has been a large consensus in the DSM community that the density of soil data is a clear limiting factor for DSM [START_REF] Lagacherie | Digital soil mapping: A state of the art, in Digital Soil Mapping with limited soil data[END_REF], few experiments showing their impacts on DSM mapping performances exist in the literature. [START_REF] Somarathna | More Data or a Better Model? Figuring what Matters Most for the Spatial Prediction of Soil Carbon[END_REF] and Wadoux et al. (2019) both observed an increase in the performances of Soil Carbon Mapping as the amount of input data increased, regardless the algorithms used to build the DSM models. However, [START_REF] Somarathna | More Data or a Better Model? Figuring what Matters Most for the Spatial Prediction of Soil Carbon[END_REF] observed that this increase in performances was lower for the highest sampling densities, which may suggest that, beyond a given threshold of density, it would not be worthwhile to add further new locations. Such a threshold was not observed by Wadoux et al. (2019), which means that it should be highly case-dependent.

More case studies should therefore be studied to gain expertise regarding fixing this threshold.

Apart from their average densities, the soil spatial samplings may vary with regard to their spatial distribution across the study areas, which may also strongly affect the prediction performances of the DSM models that use such data as input. This was acknowledged early on [START_REF] Brus | Design-based versus model-based estimates of spatial means: Theory and application in environmental soil science[END_REF], and sampling algorithms for optimizing the completeness and evenness of the input soil locations within the geographical space [START_REF] Brus | Designing spatial coverage samples using the k-means clustering algorithm[END_REF] or soil covariate space [START_REF] Minasny | A conditioned Latin hypercube method for sampling in the presence of ancillary information[END_REF] have been proposed. Similarly, Zhang and Zhu (2019) observed a slight increase in DSM performances when the spatial sampling was optimized with regard to its representativeness of the distribution of the soil covariates. [START_REF] Lark | How should a spatial-coverage sample design for a geostatistical soil survey be supplemented to support estimation of spatial covariance parameters?[END_REF] showed that geostatistical predictions can be improved by adding 10% of the closely spaced locations to regular sampling to better represent the short-range spatial structures. Finally, [START_REF] Adamchuk | Using targeted sampling to process multivariate soil sensing data[END_REF] improved the soil property predictions from soil sensors by developing a sampling strategy that considered the spread among sensor output, the local homogeneity and the physical coverage across an entire field as target spatial sampling properties. All of these works converge towards the idea that sampling strategies can be leveraged to increase DSM performances. However, a recent experiment (Wadoux et al., 2019) showed that this increase could not be obtained by some sampling strategies, which means that the benefit of sampling strategies could be highly case-dependent. This paper presents a numerical experiment that relates the characteristics of legacy spatial samplings used as input soil data (spatial density and spatial distribution characteristics) with the DSM performance and the accuracy of their ex-ante evaluation. The study used the pattern of the topsoil pseudo clay content derived from airborne Visible Near InfraRed-Short Wave InfraRed (Vis-NIR-SWIR) hyperspectral data acquired over the Cap Bon region (300 km 2 , Tunisia) at five-metre resolution [START_REF] Gomez | Using a VNIR/SWIR hyperspectral image to map topsoil properties over bare soil surfaces in the Cap Bon region (Tunisia)[END_REF]. This pattern is composed of wellpredicted clay values (R 2 = 0.75) that are free of visible artefacts and pedologically plausible, which allows it to be considered as a fair representation of the variations of a real soil property across the landscape [START_REF] Lagacherie | How far can the uncertainty on a Digital Soil Map be known?: A numerical experiment using pseudo values of clay content obtained from Vis-SWIR hyperspectral imagery[END_REF]. Such a soil dataset provided a quasiunlimited number of pseudo-measured sites, which enabled the testing of a large number (12,000) of spatial samplings. These spatial samplings were used as input data for a Quantile Random Forest algorithm that produced DSM models whose performances in predicting clay values and the associated accuracy were analysed with regard to the sampling characteristics.

Insert figure 1 here

Case study

This case study has already been described in a previous paper [START_REF] Lagacherie | How far can the uncertainty on a Digital Soil Map be known?: A numerical experiment using pseudo values of clay content obtained from Vis-SWIR hyperspectral imagery[END_REF].

Large excerpts of this paper are used in the following.

2.1.

The study area

The study area is in the Cap Bon region in northern Tunisia (36°24'N to 36°53'N; 10°20'E to 10°58'E), which is 60 km east of Tunis (Figure 1a). This 300-km² area includes the Lebna catchment, which is mainly rural (>90%). The Lebna catchment is devoted to the cultivation of cereals in addition to legumes, olive trees, vineyards and natural vegetation for animals.

The region is characterized by its rolling hills and elevations between 0 and 226 m. The climate varies from humid to semi-arid, with an inter-annual precipitation of 600 mm and an inter-annual potential evapotranspiration of 1500 mm. The soil pattern of the Lebna catchment is mainly the result of variations in lithology. The variations in the bedrock between Miocene sandstone and Marl cause large variations in the physical and chemical soil properties (Zante et al., 2005). Furthermore, the distance between successive sandstone outcrops decreases significantly as the terrain changes from the ocean to the mountains, which also causes variations in the soil property patterns [START_REF] Gomez | Using a VNIR/SWIR hyperspectral image to map topsoil properties over bare soil surfaces in the Cap Bon region (Tunisia)[END_REF]. The soil materials were redistributed laterally along the slopes during the Holocene, which adds to the complexity of the soil pattern. The numerical experiment uses an image of topsoil clay content as input. The topsoil clay contents were derived from a Vis-NIR-SWIR hyperspectral image [START_REF] Gomez | Using a VNIR/SWIR hyperspectral image to map topsoil properties over bare soil surfaces in the Cap Bon region (Tunisia)[END_REF]. The approach used to produce the data is summarized below. More details regarding the preand post-processing of the hyperspectral image can be found in [START_REF] Gomez | Using a VNIR/SWIR hyperspectral image to map topsoil properties over bare soil surfaces in the Cap Bon region (Tunisia)[END_REF].

On November 2, 2010, AISA-Dual airborne-based hyperspectral data were acquired over the study area with a spatial resolution of 5 m. The area of the image is approximately 12 km x 24 km. The AISA-Dual spectrometer measured the reflected radiance via 359 non-contiguous bands covering the 400-to 2450-nm spectral range, with 4.6-nm bandwidths between 400 and 970 nm and 6.5-nm bandwidths between 970 and 2450 nm. The radiance units were converted to reflectance units using ASD spectrometer measurements of uniform surfaces (parking lots, asphalt, concrete) that were collected at the same time during the over flight.

Topographical corrections were performed using a digital elevation model built from ASTER data and ground control points.

To isolate the bare soil areas, the study masked pixels with normalized difference vegetation index (NDVI) values greater than an expert-calibrated threshold (0.20). Water and Urban areas were also removed. Finally, the bare soil represented 46.3% of our study area, that is, 5,889,847 measured AISA-Dual 5-m x 5-m pixels.

A Partial Least Square Regression (PLSR) technique (e.g. [START_REF] Tenenhaus | La régression PLS[END_REF] was then applied to estimate the topsoil clay contents from AISA-DUAL reflectance at each location. The PLSR was calibrated from 129 couples of AISA-DUAL Vis-NIR-SWIR reflectance spectra on bare soil surfaces associated with the topsoil clay content measured on a laboratory soil sample collected from the same bare soil surfaces. Before the PLSR model was built, the reflectance was converted into "absorbance" (log [1/reflectance]). In addition, a Savitzky-Golay filter with second-order polynomial smoothing and window widths of 30 nm [START_REF] Savitzky | Smoothing and differentiation of data by simplified least squares procedures[END_REF] and a mean centring and variance scaling was applied to the spectra to reduce noise.

The calibrated PLSR model was then validated using a leave-one-out cross-validation that showed successful predictions (R² = 0.75, [START_REF] Gomez | Using a VNIR/SWIR hyperspectral image to map topsoil properties over bare soil surfaces in the Cap Bon region (Tunisia)[END_REF], figure 3). The PLSR model was then applied to all bare soil pixels to estimate the topsoil clay content, thus providing the final predicted topsoil clay properties map (Figure 1b). These treatments were implemented in R (Version 1.17) using the signal and pls packages [START_REF] Mevik | The pls Package: Principal Component and Partial Least Squares Regression in R[END_REF].

Digital Elevation Model and derivatives

= [1] = [2] = [3]
where B8, B4 and B3 are the spectral bands centred at 842 nm, 665 nm, and 490 nm, respectively.

NDVI was selected as a proxy of the vegetation health and biomass that could be in relation with topsoil clay contents through the water and nutriment retention properties. CI and RI were selected as proxies of topsoil colour that could be related to clay content through different parent materials.

Methods

DSM modelling

Three criteria were considered to select the learning algorithm used to produce the DSM models: i) the algorithm should be one of the most used and the most efficient among the recent DSM applications, ii) the algorithm had to provide local uncertainty predictions to be able to test its ability to predict the associated uncertainty, and iii) the algorithm should be run without manual intervention and repeated a great number of times in the numerical experiment. Combining these three criteria resulted in selecting the Quantile Regression Forest as the learning algorithm used in this study.

Recent performance testing has found that the Random Forest, from which the Quantile Regression Forest is derived, was among the best algorithm for obtaining predictions of soil properties [START_REF] Nussbaum | Evaluation of digital soil mapping approaches with large sets of environmental covariates[END_REF], which confirmed a test performed on a wider range of machine learning applications [START_REF] Caruana | An empirical comparison of supervised learning algorithms[END_REF]. The Quantile Regression Forest has also been demonstrated as efficient for predicting the uncertainty associated with soil property predictions [START_REF] Vaysse | Using quantile regression forest to estimate uncertainty of digital soil mapping products[END_REF]. Finally, many authors have adopted this machine learning algorithm for recent large-scale DSM applications [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF][START_REF] Dobarco | Uncertainty assessment of GlobalSoilMap soil available water capacity products: A French case study[END_REF].

Random Forests and Quantile Regression Forests

This section summarizes the main characteristics of Random Forests and Quantile Regression Forests, using excerpts of [START_REF] Meinshausen | Quantile Regression Forests[END_REF]. It was already presented in [START_REF] Lagacherie | How far can the uncertainty on a Digital Soil Map be known?: A numerical experiment using pseudo values of clay content obtained from Vis-SWIR hyperspectral imagery[END_REF]. More details on these two machine learning algorithms are given in the seminal papers by [START_REF] Breiman | Random Forests[END_REF] and [START_REF] Meinshausen | Quantile Regression Forests[END_REF], respectively.

Let Y be a real-valued response variable and X be a covariate or predictor variable that is likely high-dimensional. A standard goal of statistical analysis is to infer the relationship between Y and X. Random Forests grow a large (>500) ensemble of trees using n independent observations (Yi , Xi), i = 1, . . ., n. Each tree grows via a recursive partitioning of the source set using one predictor variable X. At each step, the source set is split into two subsets following a test on the value of X. When Y is a quantitative variable, the selected test is the one that minimizes the within-subset variance of Y [START_REF] Breiman | Classification and regression trees[END_REF]. The recursive partitioning is limited by a stopping rule, and the subsets are produced by the last split being the leaves of the tree. The ensemble of trees is produced by using a random sample of the training data and a random subset of the predictor variables for each tree.

For the regression, the prediction of a single tree ?? of a Random Forest for a new data point x can be represented as the weighted average of the original observations Yi, i = 1, . . ., n:

= ∑ , [4]
where , is the weight vector given by a positive constant that is 1 if the observation Yi is part of the same leaf and 0 otherwise.

By using Random Forests, the prediction is the average prediction of k single trees that were constructed as described above.

= !

[5] with = " ∑ , # $ [6]
One could assume that the weighted observations deliver a good approximation not only of the conditional mean but also of the full conditional distribution. This assumption is at the heart of the Quantile Regression Forest algorithm, which estimates the conditional distribution function of Y given x via:

% &| = ∑ 1 )* + ,-.

[7]

From this conditional distribution, it is possible to derive both the predicted value (the mean) and the bound of the 90% prediction interval that predicts the associated uncertainty (the 0.05 and 0.95 quantiles).

Numerous implementations of RF and QRF now exist. Ranger Package (Wright and Ziegler, 2017) was selected because it is a fast implementation of Random Forests that is suitable for multiple model building.

Tuning Random Forest Parameters

The Random Forest Algorithm has several hyperparameters that must be set by the user.

Among them, three parameters may significantly impact the results and therefore should be tuned to improve the predictions [START_REF] Probst | Hyperparameters and Tuning Strategies for Random Forest[END_REF]: i) the number of observations drawn randomly for each tree, ii) the number of variables drawn randomly for each split and iii) the minimum number of samples that a node must contain. These parameters were tuned using one of the most established tuning strategies, sequential model-based optimization [START_REF] Jones | Efficient global optimization of expensive black-boxfunctions[END_REF][START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF]. This tuning algorithm iteratively uses the results of the different already evaluated hyperparameter values and chooses future hyperparameters based on these results. It is implemented in the TuneRanger Package [START_REF] Probst | Hyperparameters and Tuning Strategies for Random Forest[END_REF].

After some trials, 100 iterations appeared to be a good compromise that ensured a fairly good convergence towards an optimized solution while being acceptable in terms of computing costs.

Insert figure 2 here 3.2. The numerical experiment 3.2.1. General approach The workflow of the numerical experiment is presented in Figure 2. First, a master evaluation set of 100,000 locations with pseudo values of clay content was randomly selected from the total set of pixels. This master set served to determine the real performances of the tested models through a set of indicators described further in 3.2.3.

The remaining locations were used to build and evaluate 12,000 models, i.e., 1,000 models for each of the 12 considered sample sizes (100, 200, 300, 400, 500, 600, 800, 1,000, 2,000, 3,000, 5,000, and 10,000 locations), each corresponding to a given average spacing ( 1732 . [8] For a given average spacing, the soil inputs of the 1000 models were selected by a specific sampling procedure (see the next section) that mimicked the more or less uneven spatial distributions observed when using legacy data as soil inputs. The soil inputs were then randomly divided into two sets, keeping 75% of the locations to produce DSM models using the QRF algorithm and 25% to obtain a so-called independent evaluation set as currently practised in DSM applications. The DSM models were obtained using QRF according to the method presented above and then evaluated from the independent datasets. The same performance indicators used for the master set were considered to enable statistical comparisons.

Insert figure 3 here

The sampling procedure

The sampling procedure was designed to randomly produce spatial samplings with contrasting degrees of unevenness. The study area was first stratified into 25 geographical strata of equal area using a K-means classification of the locations [START_REF] Walvoort | An R package for spatial coverage sampling and random sampling from compact geographical strata by k-means[END_REF].

Then, the following algorithm was applied:

1) Define a given size of sampling N 2) Select at random an integer P between 1 and 25 (the number of strata) and select at random P strata among the 25 3) Select N/2 locations by a stratified random sampling, using the P strata selected in step 2 4) Complete the spatial sampling by N/2 locations selected at random over the entire study area For a given size N, steps 2 through 4 were repeated 1,000 times to obtain 1,000 spatial sampling that differed in unevenness thanks to different random selections of P. The procedure was repeated for each selected sample size.

Figure 3 shows four examples of different spatial samplings provided by this procedure.

Model performance indicators

Four model performance indicators were considered. The first two were the Mean Square Error (MSE) and the 90% Prediction Interval Coverage probability (PICP90), which were calculated using the master evaluation set of 100,000 sites.

MSE >?E = . & G -y J [9]
with ŷi as the predicted value, yi as the observed value of Clay, and n = 100,000

PICP90 = . card ) & G ≥ y UVWXYZ and & G ≤ y UVWX]Z . J [10]
with y UVWXYZ and y UVWX]Z as the lower and upper bounds of the predicted 90% confidence interval PICP90 expressing the probability that all observed values fit within the 90% prediction limits provided by the DSM model.

The last two indicators aimed to quantify the relative errors of Mean Square Error provided either by the QRF algorithm as a by-product (MSEmod) or by removing 25% of the soil inputs taken as an independent evaluation set (MSErmv). To calculate these errors, the MSE calculated on the master evaluation (MSEref) set served as a reference. were calculated in the covariate space, in the geographical space and with regard to the target variable (Clay content). The last two (variance and spatially structured-Variance ratio)

were calculated for Clay content only. These indicators are detailed below.

Coverage index

The Coverage index (CI) [START_REF] Gunzburer | Uniformity measures for point samples in hypercubes[END_REF] expresses the degree of unevenness of the spatial distribution of the locations of the spatial sampling.

= : l m ∑ no -o ̅ r J s [13]
with o the distance of site i to its nearest neighbour and o ̅ the mean value of the o Note that, for a regular mesh, CI = 0. Then, a small value of CI means that the spatial sampling has a spatial distribution close to that of a regular grid. The R package DiceDesign [START_REF] Dupuy | DiceDesign and DiceEval: Two R Packages for Design and Analysis of Computer Experiments[END_REF] was used to calculate CI.

Kullback-Leibler Divergence (KLD)

The Kullback-Leibler divergence (KLD) [START_REF] Kullback | On information and sufficiency[END_REF] measures how close two probability distributions are. It was used here to measure the distance between the calibration or evaluation sets and the master validation set. It is a measurement of the representativeness of the calibration and validation sets with regard to the master validation set, which is assumed to represent the real distribution of the different variables across the study region.

For distributions P and Q defined in the same probability space, the Kullback-Leibler divergence between P and Qis defined to be (Wikipedia, https://fr.wikipedia.org/wiki/Divergence_de_Kullback-Leibler)

tu v|| w = -∑ v x∈z . {|3 m } x ~ x s [14]
For distributions P and Q of continuous random variables, the Kullback-Leibler divergence is

tu v||w = • 6 . log m ƒ x " x s … … [15]
The KLDs are first calculated for each variable and then averaged to obtain a unique value for the geographical and the covariate spaces.

Percentage of out of range

The percentage of out of range measures the proportion of locations of the study region having values of variables (in the covariates space, the geographical space of Clay values) that are out of the range of those of the sites included in the calibration or evaluation sets. It is also a measurement of representativeness but differs from the former due to its focus on extreme values.

% ‡ˆ‰_|‹_ /931 = OE<>`m•-+ Ž••' ' -& -+ "••- ' --s OE<>` ˜ [16]
with P the calibration or evaluation sets and R the set of locations in the study region

Clay Variance and percentage of spatially structured variance ratio

Variance of clay is a measure of dispersion within the calibration or validation dataset with regard to clay values.

A spatially structured variance ratio (SSVR) was proposed by [START_REF] Vaysse | Evaluating Digital Soil Mapping approaches for mapping GlobalSoilMap soil properties from legacy data in Languedoc-Roussillon (France)[END_REF] as the complement to 1 of the nugget-to-sill ratio [START_REF] Kerry | Determining nugget : sill ratios of standardized variograms from aerial photographs to krige sparse soil data[END_REF]. It indicates the proportion of the spatially structured variance that is captured by the model.

SSVR = variance -nugget / variance. [17]

Because it was impossible to fit variograms on 12,000 trials, the nugget could not be calculated directly. It was therefore approximated by computing the semi-variances at lags centred on the average spacing ([average spacing -100 metres, average spacing + 100 metres])

Insert figure 4 here

Results

Impact of average spacing

Figures 4a through 4d show the evolution of the different indicators of DSM performances

with the average spacing (Equation 8). The Mean Square error on the predicted Clay value (MSEref) covered a large range of values across the 12,000 models (Figure 4a). These values increased regularly with the average spacing: from MSE = 7,911 g 2 .kg -2 (68% of explained variance) to MSE = 18,882 g 2 kg -2 (22% of explained variance). The amount of variation of MSEref for a given size of sampling also increased regularly with the average spacing.

PICP90 exhibited a positive bias (overestimated uncertainty) with regard to the expected 90% value for the smallest average spacings (below 612 m) and a negative one (underestimated uncertainty) beyond this threshold (Figure 4b). However, the errors were only important (more than 1%) for the largest average spacing (-2.5% for 1732 m) and for the smallest ones (between 1.2 and 1.7 for average spacing at and below 316 m). Apart from the influence of the average sampling variations, great variabilities of performances for PICP90 estimations were observed for the largest average spacings (see the bars of Figure 4b).

The bagging procedure of the QRF algorithm and the external evaluation using 25% of removed sites both had a negative bias (Figures 4c and4d) regardless of the average spacing, which revealed systematic underestimations of the overall uncertainty of the DSM models (e.g., -17% and -10% respectively for 1732-m spacing). This bias seemed not as closely related to the average spacing of sites as was observed previously for the other indicators, although slight decreases could be observed beyond 866-m spacing. The most important variations were observed within each sampling size, as shown by the large bars of Figures 4c and4d.

Insert figure 5 here

The impacts of the spatial distributions of sites

The matrices of Figures 5a through 5d Whatever the considered DSM performance indicators, the correlation coefficients were highly variable with regard to the spatial sampling indicators and the average spacing (between 0.00 to 0.91). These coefficients tended to increase with the decrease of average spacing, which is summarized by the increase in the strength of the relation between the indicators of performance and the spatial sampling indicators (last lines of the matrices).

However, some noticeable exceptions occurred for some relations between performance indicators and spatial sampling indicators.

MSEref (Figure 5a), which expressed the ability of the DSM models to predict the correct value of clay content, was strongly positively correlated with the coverage_index and Kullback-Leibler divergence calculated in the geographical space, which means that the performances of the DSM models were better for evenly geographically distributed and representative spatial samplings. These correlations increased greatly when average spacing decreased (until -0.87 for the two indicators) but were already substantial for the largest spacing (-0.42 and -0.38). For all of the tested average spacings, MSEref was also moderately correlated with respectively). This means that the performances of the DSM models tended to increase as the clay values included in the spatial sampling were largely dispersed, and this well covered the range of clay values of the study region. The three indicators calculated in the covariate space were only significantly correlated with MSEref for the smallest average spacings. Furthermore, these correlations were always smaller than those obtained by the same indicators calculated in the geographical space. Finally, the spatially structured variance ratio (SSVR) exhibited moderate negative correlations for intermediate values of average spacing (between -0.27 and -0.42), which means that, for these values, the more the spatially structured variability (especially the short range one) was captured by the spatial sampling, the better the performances were. PICP90, ^22 _;` and ^22 >_j (Figures 5b,5c and 5d), which all expressed the ability to predict the uncertainty associated with the predicted values given by the DSM models, behaved similarly to each other with regard to the correlations with the spatial sampling indicators, with however stronger overall correlations for ^22 >_j than for PICP90 and for ^22 _;` than for ^22 >_j (see the stepwise regression coefficient, the last lines of the matrices in Figures 5b, 5c and5d). Contrary to MSEref, clear differences of correlation rankings were observed between the smallest and the largest average spacings. As far as the former are concerned, PICP90, ^22 _;` and ^22 >_j were predominantly correlated with Variance_Clay and %-outof-range_Clay (between 0.46 and 0.76 and between -0.38 and -0.62, respectively, for average spacing larger than or equal to 866 m). This means that the uncertainty was much better predicted when the clay values included in the spatial sampling were highly dispersed and covered well the range of clay values of the study region. At the smallest spacings, the strongest correlations were observed with the coverage_index and Kullback-Leibler divergence calculated in the geographical space (between -0.49 and -0.88 and between -0.48 and -0.85, respectively, for average spacing smaller than or equal to 387 m), which means that evenly geographically distributed and representative spatial samplings enabled an accurate prediction of the DSM model uncertainty. The spatial sampling indicators calculated on the covariate space only exhibited substantial correlations with PICP90, ^22 _;` and ^22 >_j at the smallest spacings. Finally, SSVR exhibited substantial correlations only for PICP90 (between 0.33 and 0.45).

Insert figure 6 here

Discussion

The impact of the average spacing

All of the results confirmed that the average spacing, related with the size of the calibration data sets used as input for the DSM approach, strongly impacted the results of a DSM approach. As already observed by [START_REF] Somarathna | More Data or a Better Model? Figuring what Matters Most for the Spatial Prediction of Soil Carbon[END_REF]Wadoux et al. (2019), we observed (Figure 4a) a clear decrease of prediction errors (MSEref) when the average spacing is decreasing. It should be noticed that, thanks to the use of pseudo-values of clay content

given by the hyperspectral image, we explored a more complete range of average spacing, which allowed us to analyse a larger range of model performances (from 22% to 68% of explained variance) that covered fairly well the ones cited in the literature. Figure 4a revealed that the decrease of prediction errors with average spacing was not linear. By substituting the average spacing X axis of Figure 4a by the spatially structured variance ratio (SSVR, see section 3.2.4), a new insight into this average proportion of clay variance was revealed (Figure 6): for the largest average spacing, that captured the least spatially structured variance (<= 800 sites, average spacing = 707 m), the average increase in performances was perfectly linear, whereas further increases of this ratio provided gains that were smaller and smaller than the previously observed linear trend. Therefore, this observed threshold separated two contrasting situations: below the threshold of average spacing of 707 m, the spacing of sampled sites was the only limiting factor, while beyond the threshold, other limiting factors, such as the precision of the covariates, also played a role in the quality of the results. This contrasting behaviour could explain why contradictory results have been obtained recently regarding the impact of the spatially structured variance ratio on DSM results, as observed by [START_REF] Vaysse | Evaluating Digital Soil Mapping approaches for mapping GlobalSoilMap soil properties from legacy data in Languedoc-Roussillon (France)[END_REF] and not observed by [START_REF] Nussbaum | Evaluation of digital soil mapping approaches with large sets of environmental covariates[END_REF]. It may also explain why improving a covariate dataset often do not significantly improve the DSM products when overly sparse test datasets are used to calibrate the DSM model [START_REF] Samuel-Rosa | Do more detailed environmental covariates deliver more accurate soil maps?[END_REF][START_REF] Loiseau | Satellite data integration for soil clay content modelling at a national scale[END_REF].

Apart from its influence on the prediction error, the average spacing of sites also had an impact on the estimations of the associated uncertainty, which, to our knowledge, has not been observed before. Decreasing the average spacing reduced the underestimation of the overall uncertainty provided either by the bagging procedure of the Random Forest (Figure 4c) or the model-free evaluation process (Figure 4d), although these reductions were not as clear and regular as the prediction error because larger variabilities of results were observed within each tested spacing (see the error bars in Figures 4c and4d). A more complex behaviour was observed for the local estimation of the confidence interval by the model tested by PICP90 (Figure 4b). For the largest spacings, PICP90 estimates converged towards the nominal value of the confidence interval (90%) as the spacing decreased, whereas for the smallest spacings, PICP90 moved away from this nominal value. This latter unexpected result could be interpreted as the inclusion of outliers as the spacing decreased, which could perturb the estimates of the confidence interval bounds. It must be noticed that, as for the estimations of the overall uncertainty evoked before, a large variability in estimating PICP90

was observed within each tested spacing.

Finally, a final effect of the spatial sampling size is that it changes the amount and the drivers of the variations of performances observed within each sampling size. This will be developed in the next section.

The impact of the distribution of sites over the study region

The bars on Figures 4a through 4d show that the average spacing is not the only driver of DSM performance, especially with regard to the ability of DSM approaches to estimate overall (Figures 4 c and 4 d) and local uncertainties (Figure 4b). The matrices of Figure 5 confirmed many of the underlying hypotheses of the sampling strategies that have been proposed in the literature while providing new insights on the relation between spatial sampling and uncertainty estimation and nuancing the importance of some sampling characteristics according to the size of the spatial sampling.

From the matrices of Figure 5, it clearly appeared that the regularity of sampling and the representativeness in the geographical space improved the DSM results whatever the size of the spatial sampling and the considered indicators of performances. Therefore, the legacy soil data that are often characterized by both under-sampled and over-sampled sub-regions should be ideally completed using sampling strategies that could mitigate this irregularity of sampling in the geographical space [START_REF] Brus | Sampling for validation of digital soil maps[END_REF][START_REF] Adamchuk | Using targeted sampling to process multivariate soil sensing data[END_REF]. This would require harmonizing the legacy and the new dataset techniques for removing biases caused by differences of dates, field protocols and laboratory methods [START_REF] Baume | A geostatistical approach to data harmonization -application to radioactivity exposure data[END_REF][START_REF] Ciampalini | Detecting, correcting and interpreting the biases of measured soil profile data: A case study in the cap bon region (Tunisia)[END_REF]. An alternative to adding samples should be to better take into consideration in the DSM modelling the perturbing effects of the clusters of sites. This could be done by assigning different weights to the input sites according to their degree of remoteness [START_REF] Bel | CART algorithm for spatial data: Application to environmental and ecological data[END_REF], applying resampling techniques (Richer- [START_REF] Richer-De-Forges | Probability mapping of iron pan presence in sandy podzols in South-West France, using digital soil mapping[END_REF]Taghizadeh-Merjadhi et al., 2020) or restricting the predictions inferred from each cluster of sites to representative areas corresponding to well-identified and well-mapped soil systems [START_REF] Lagacherie | Mapping of reference area representativity using a mathematical soilscape distance[END_REF].

The spread of the spatial sampling with regard to the values of the target soil property (%out-of-range and variance of Clay in Figure 5) also seemed crucial for improving both the predictions of the soil property and the predictions of the associated overall and local uncertainties. However, correcting the existing legacy sample with regard to this characteristic is an uneasy task because it requires additional knowledge to anticipate the locations of the extreme values of the targeted soil property that should be preferentially sampled. A local pedological knowledge or a proxy of the soil property of interest [START_REF] Adamchuk | Using targeted sampling to process multivariate soil sensing data[END_REF] could be mobilized for that.

It is interesting to note that the sampling characteristics that involved the soil covariates only had an impact on the results for the smallest spacing of spatial sampling. In these cases, this impact was increased by strong relationships established between the covariates and the predicted soil property, whereas for the largest spacing, these relationships were too weak, even if the most related covariates were selected (results not shown in this paper).

This means that the strategies of sampling based on the regularity of coverage of the spatial sampling in the covariate spaces (Minasny and Mc Bratney, 2006;[START_REF] Carré | Estimation and potential improvement of the quality of legacy soil samples for digital soil mapping[END_REF]Zhang and Zhu, 2019) could not be effective for correcting overly sparse legacy datasets.

Conversely, for large datasets for which more covariates were involved in the model, a fair distribution of the covariates values would be required. However, the coverages and KLD indices in the covariate space and the geographical one were found to be highly correlated for the smallest average spacing (r > 0.90 for average sampling >= 1225 m), which means that taking into account the covariate space would be of little interest if the regularity of sampling in the geographical space is already ensured. However, this result would not hold in particular pedological contexts characterized by small inclusions of land with contrasted values of covariates that could be missed by regular samplings in the geographical space.

Finally, the stepwise regression coefficients of determination given in the last lines of the matrices of Figure 4 clearly showed that the selected sampling characteristics could not alone explain the variations of performances that were observed across the 12,000 trials.

This was particularly true for the sparsest spatial sampling when we considered the mean square error on predicted value (Figure 5a) or the biases of estimation of the latter (Figure 5c or 5d). The reverse tendency was observed for the error on PICP, for which the smallest average spacing of sites obtained the lowest coefficient of determination.

Uncertainty estimation biases

The two tested procedures of overall uncertainty estimation -the model-based and the model-free ones -exhibited non-negligible biases on uncertainty of clay content predictions (Figures 4c and4d). Although the two procedures could be considered as intrinsically unbiased, some specific characteristics of the legacy spatial sampling to which they were applied were responsible of these biases. Although it could be observed (Figure 4) that the sparse samplings were more prone to bias than the dense ones, the average spacing did not seem to be a first-order driver because much more variability occurs within a given sampling size (see the bars of Figures 4c and4d).

The correlations matrices (Figures 5c and5d) provided some insights on the causes of such biases. As far as the largest spacing were concerned, biases were all the more great that the sampling underestimated the real variations of clay content and thus left aside their extreme values. This observation can be related with the general difficulty of the inference models, such as Random Forest, to predict values that are out of the range of their learning sample (Conn et al., 2014). Alternately, overly clustered datasets (see examples in Figure 3, left column) resulted in selecting evaluation sites that could be too close from the calibration sites for satisfying the condition of independence , which may induce underestimations of the prediction errors. Therefore, to estimate the overall uncertainty as well, it is important to mitigate the perturbing effects of the clusters of sites techniques cited above.

Finally, it is worth noting that the estimations of the local uncertainty through a confidence interval calculated by the QRF algorithm seemed to be more robust than the estimations of the overall uncertainty (Figure 3b). Only the largest spacing (1732 m) and, to a lesser extent, smallest ones gave unsatisfactory results. This confirmed the results obtained by Vaysse and Lagacherie (2017) using the same algorithm.

Limitations and open questions.

Although some clear and coherent tendencies could be retrieved from the results of this numerical experiment (see above), some open questions remain. First, the characteristics of the spatial sampling that were considered in this paper did not explain the entire variability of DSM performances. The weak statistical relations observed in Figure 5 for the largest spacing suggest that some hidden factors should be evoked. Among others, we hypothesize that the random process used for optimizing the hyperparameters of the Random Forest generated a noise on DSM performances, the best possible combination of parameters not always being reached because of local optimal solutions, especially when the size of the learning sample is small. This hypothesis is supported by the fact that the average variabilities of the optimal QRF parameters provided by the optimization process decreased as the size of the sampling increased (average Coefficient of Variation from 35% to 21%).

Second, although a large range of soil sampling spacings were explored in this case study, the size of the study areas limited the testing of the sparser soil datasets that fed the DSM applications conducted at national (e.g. [START_REF] Mulder | GlobalSoilMap France: High-resolution spatial modelling the soils of France up to two meter depth[END_REF], continental [START_REF] Ballabio | Mapping topsoil physical properties at European scale using the LUCAS database[END_REF] or global scale [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF]. Whether or not the trends exhibited in [START_REF] Stuffler | The EnMAP hyperspectral imager -An advanced optical payload for future applications in Earth observation programmes[END_REF][START_REF] Steinberg | Prediction of Common Surface Soil Properties Based on Vis-NIR Airborne and Simulated EnMAP Imaging Spectroscopy Data: Prediction Accuracy and Influence of Spatial Resolution[END_REF]the Hyperspectral Infrared Imager -HyspIRI-, Lee et al., 2015), it could be envisaged to reproduce the same numerical experiment at a wider extent with, however, a loose spatial resolution. In the absence of such an experiment, the results obtained for the largest average spacing considered in this numerical experiment (1732 m) should help in orienting the future design of DSM approaches at these largest extents.

Conclusions

The main lessons of the numerical experiment are as follows

• User and producers of DSM products should be aware that the current methods of evaluation tend to underestimate the overall uncertainty, especially for sparse and unevenly distributed soil sampling

• Although decreasing the average spacing of soil inputs always brings improvements of DSM performances, one should be aware that, beyond a given threshold of average spacing, the improvement would need also to collect better soil covariates.

• The spatial distributions of the legacy data and the sampling strategies for correcting these distributions play a key role in reaching the best DSM performances. Sampling strategies that provide complete and even distributions in the geographical space and have as great a spread of the target soil property as possible should be privileged.

• Some hidden sampling characteristics that were not considered in this experiment seem to play a significant role, especially for sparse sampling. More research is required for identifying these characteristics.

Wright, M.N., Ziegler, A., 2017. ranger : A Fast Implementation of Random Forests for High Dimensional Data in C ++ and R 77.

Zhang, G., Zhu, A.X., 2019. A representativeness heuristic for mitigating spatial bias in existing soil samples for digital soil.pdf. Geoderma 351, 140-163.

A

  30-m ASTER digital elevation model (DEM) with specific ortho-rectification and mosaicking was produced for this area. The classical geomorphometric indicators found in the DSM literature were calculated. These include Elevation, Slope, Aspect, plan Curvature, Profile Curvature, Multi-Resolution Valley Bottom Flatness (MRVBF) and four variables describing the aspect: northness, easterness, northwesterness, and northeasterness. All of these indicators served as covariate candidates for representing the relationships between clay content and the relief. 2.2.3. Tunisian agriculture map A set of layers of the Tunisian agriculture map (STUDI-SCOT-SODETEG, 2001) at 1:20,000 scale was considered as covariates. This includes two soil properties, soil colour classes and textural classes, which were mapped by considering existing detailed soil maps completed by manual interpretation of remote sensing images. The other covariates extracted from the Tunisian agriculture map were the parent material (14 classes) and the land use (7 classes). 2.2.4. Indices from Sentinel-2 images Three spectral indices were used in this work, based on a Sentinel-2 image acquired on the 2 nd of November 2016 over the study area. The acquisition date of Sentinel-2 data was chosen to fit in the period of more extensive bare soils. Sentinel-2 data are composed by multispectral data in 13 bands covering the Visible, Near InfraRed and Short Wave InfraRed spectral domain with spatial resolutions ranging from 10 to 60 m. The Sentinel-2 data acquired on the 2 nd of November 2016 was downloaded from the Muscate platform of the French land data centre, called Theia (https://www.theia-land.fr/) in Level 2A, i.e., corrected from atmospheric effects, thanks to the three bands acquired at 60-m spatial resolution (coastal at 443 nm, water vapour at 1375-nm atmospheric correction, coastal at 443 nm, and cirrus at 1376 nm). The three spectral indices used in this work were the normalized difference vegetation index (NDVI), Redness index (RI) and Colour index (CI), calculated following Pouget et al. (1990) and Ghodalizeh et al. (2016):

  m, 1225 m, 1000 m, 866 m, 775 m, 707 m, 612 m, 548 m, 387 m, 316 m, 245 m and 173 m) using the following equation: /012/31 56/7893 = : $;$<= <>?< @ A?

  indicators quantifying the characteristics of the spatial distributions of the tested spatial sampling were calculated as candidates to explain the performances of the model summarized by the performance indicator presented above. The first three indicators (Coverage Index, Kullback-Leibler Divergence and percentage of out of range)

  show, for each average spacing (the columns of the matrices), the correlations between the spatial sampling indicators (lines of the matrices) and the indicators of DSM performance (one matrix per indicator). The last lines of the matrices show the coefficients of determination of the stepwise regressions between the indicator of performance of interest and the set of indicators of spatial sampling, which allowed the strength of the relation between these two types of indicators to be appreciated.

  be extrapolated to these applications remains an open question. With the next availability of hyperspectral VIS-NIR-SWIR satellite data (such as the French HYPerspectral X Imagery -HYPXIM-, Briottet et al., 2013; the Spaceborne Hyperspectral Applicative Land and Ocean Mission -SHALOM-, Bussoletti, 2012; the German Environmental Mapping and Analysis Program -EnMAP-,

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Location of the study area (a) and the spatial pattern of pseudo values of topsoil clay content (b) (afterLagacherie et al, 2009)

Figure 4 :Figure 5 :

 45 Figure4: Evaluation of DSM models (quantile Random Forests) using different size of soil input data (size is expressed by spacing): a) Mean Square error on predicted Clay value (MSEref), (in g 2 /kg 2 ).The green line is the total variance of Clay content over the study area b) 90% Prediction Interval Coverage Index (PICP90) (in %). The green line is the expected value of 90%. c) error on the QRF based estimation of MSEref (in % MSEref). The green line is 0 (no error) d) error on the model-free estimation of MSEref (removing 25% of the soil inputs for validation) (in % MSEref), ). The green line is 0 (no error). Red dots are averaged values per spacing and bars are +-the standard deviations (1000 models per spacing)

Figure 6 :

 6 Figure6: Evolution of the average performances of predictions (mean MSEref) with the spatially-structured variance ratio (SSVR). The red line is the linear regression using the seven spatial sampling (out of 12) having the smallest values of SSVR.

  b

Acknowledgements

This research was conducted within the "Centre d'Expertise Scientifique Cartographie Numérique des sols" granted by the CNES-TOSCA programme.

Dominique

Arrouays is coordinator and Philippe Lagacherie and Hocine Bourennane are collaborators of the GLADSOILMAP Consortium supported by the LE STUDIUM Loire Valley Institute for Advanced Research Studies