PhD) a , Cédric Sautier Isabelle Rivals
email: isabelle.rivals@espci.fr

MS Guillaume Cognon

MS Vincent Cuzuel

EVALUATION OF DISTANCE-BASED APPROACHES FOR FORENSIC COMPARISON: APPLICATION TO HAND ODOR EVIDENCE

Keywords: Bayesian inference, dissimilarity measure, forensic science, human hand odor, likelihood ratio, logistic regression

The issue of distinguishing between the same-source and different-source hypotheses based on various types of traces is a generic problem in forensic science. This problem is often tackled with Bayesian approaches, which are able to provide a likelihood ratio that quantifies the relative strengths of evidence supporting each of the two competing hypotheses. Here, we focus on distance-based approaches, whose robustness and specifically whose capacity to deal with high-dimensional evidence are very different, and need to be evaluated and optimized.

A unified framework for direct methods based on estimating the likelihoods of the distance between traces under each of the two competing hypotheses, and indirect methods using logistic regression to discriminate between same-source and different-source distance distributions, is presented. Whilst direct methods are more flexible, indirect methods are more robust and quite natural in machine learning. Moreover, indirect methods also enable the use of a vectorial distance, thus preventing the severe information loss suffered by scalar distance approaches.

Direct and indirect methods are compared in terms of sensitivity, specificity and robustness, with and without dimensionality reduction, with and without feature selection, on the example of hand odor profiles, a novel and challenging type of evidence in the field of forensics. Empirical evaluations on a large panel of 534 subjects and their 1690 odor traces show the significant superiority of the indirect methods, especially without dimensionality reduction, be it with or without feature selection.

Introduction

A generic problem in courts of law is to decide whether a trace of an unknown origin, often drawn from a crime scene, and a specimen from a known source, stem from the same source, for example a person or a firearm. If the source is a person, the traces might be biometric such as a DNA profile [START_REF] Aitken | Statistics and the Evaluation of Evidence for Forensic Scientists[END_REF][START_REF] Puch-Solis | Practitioner Guide No 2: Assessing the Probative Value of DNA Evidence, Guidance for Judges, Lawyers, Forensic Scientists & Expert Witnesses[END_REF], fingerprints [START_REF] Neumann | Quantifying the weight of evidence from a forensic fingerprint comparison: A new paradigm[END_REF], a voice [START_REF] Morrison | A comparison of procedures for the calculation of forensic likelihood ratios from acoustic-phonetic data: Multivariate kernel density (MVKD) versus Gaussian mixture model-universal background model (GMM-UBM)[END_REF], an olfactory profile [START_REF] Cuzuel | Origin, analytical characterization and use of human odor in forensics[END_REF], or they might consist of footwear impressions [START_REF] Tang | Computational methods for the analysis of footwear impression evidence[END_REF], handwriting [START_REF] Tang | Likelihood ratio estimation in forensic identification using similarity and rarity[END_REF], etc. If the source is a firearm, the traces may be features such as striations and impressions on a bullet or on a cartridge case [START_REF] Mattijssen | Validity and Reliability of Forensic Firearm Examiners[END_REF].

The most common approach of forensic science to this problem is to estimate a likelihood ratio (LR), i.e. the ratio of the joint probability of occurrence of the two traces under the hypothesis that they arose from the same source and under the hypothesis that they arose from different sources. A convenient solution is to replace the joint probability of the traces by the probability of a distance between the two traces quantifying their dissimilarity [START_REF] Neumann | Computation of likelihood ratios in fingerprint identification for configurations of any number of minutiae[END_REF][START_REF] Riva | Automatic comparison and evaluation of impressions left by a firearm on fired cartridge cases[END_REF][START_REF] Tang | Computational methods for the analysis of footwear impression evidence[END_REF][START_REF] Ali | Sampling variability in forensic likelihood-ratio computation: A simulation study[END_REF][START_REF] Muehlethaler | Evaluation of infrared spectra analyses using a likelihood ratio approach: A practical example of spray paint examination[END_REF][START_REF] Riva | Comparison and interpretation of impressed marks left by a firearm on cartridge cases -Towards an operational implementation of a likelihood ratio based technique[END_REF][START_REF] Vergeer | A method for forensic gasoline comparison in fire debris samples: A numerical likelihood ratio system[END_REF]. If, as is most often the case, the distance is scalar, there is an important loss of information. Thus, we choose to focus on distance-based methods, but with the possibility to use a vectorial distance between traces. Furthermore, the distance-based LR estimate can be obtained either directly, by estimating the distance likelihoods under the two hypotheses, or indirectly, by first using logistic regression to discriminate between same-source and different source distance distributions, and then Bayes' formula to infer the LR. Whilst the direct method is more flexible, the indirect method is more robust and quite natural in machine learning [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. It is sometimes advocated for in the forensic context, for the same reasons and also because it enables score calibration and fusion with minimal mathematical complexity [START_REF] Enzinger | A demonstration of the application of the new paradigm for the evaluation of forensic evidence under conditions reflecting those of a real forensicvoice-comparison case[END_REF][START_REF] Morrison | Tutorial on logistic-regression calibration and fusion: converting a score to a likelihood ratio[END_REF]. Here, we show that the indirect method also enables the use of a vectorial distance, thus preventing the severe information loss suffered by scalar distance approaches. We discuss the direct and indirect methods in terms of robustness and ability to handle high dimensional evidence, with or without dimensionality reduction, and with or without feature selection. We evaluate them in terms of sensitivity, specificity and robustness on the example of traces consisting of a hand odor profile.

Materials and methods

Problem statement

The aim is, given the evidence consisting of a pair of traces (e.g. two olfactory profiles), to decide whether these traces have the same source (e.g. the same person) or not. In the following, Hss refers to the hypothesis that the two traces stem from the same source, and Hds to the alternative hypothesis that they stem from different sources. Given the a priori probabilities P(Hss) and P(Hds), the Bayesian formula yields the posterior probability of Hss given the evidence E:

(1) where f(E|Hss) and f(E|Hds) are the distributions of the evidence under Hss and Hds, or likelihoods. Jeffreys developed an absolute scale to evaluate the degree of confidence in the same-source hypothesis outside a decisional framework based on the posterior probability of Hss using the LR [START_REF] Jeffreys | Theory of probability[END_REF][START_REF] Robert | The Bayesian Choice: from Decision-Theoretic Motivations to Computational Implementation[END_REF] defined as:

(2) which is independent of the a priori probabilities. In fact, the observed evidence E consists of the two traces which are represented by n-dimensional vectors (whose components are the amounts of each odor compound). Since we focus on distance-based methods, the information contained in these two vectors is transformed into a distance or dissimilarity measure, which can be either a scalar or an n-vector (a distance for each feature of the trace,

P H E f E H P H f E H P H f E H P H

= + ( | ) ( | ) ( ) ( ) ( ) ( ) ( ) LR E f E H f E H ( ) ( | ) ( | ) ss ds =
here for each odor compound). In the following sections, this distance between the two traces will be denoted by d.

Candidate methods

The LR can be obtained either directly, i.e. by estimating the likelihoods of the distance under the two competing hypotheses, or indirectly, i.e. using first logistic regression to discriminate between same-source and different source distance distributions, and then formulas (1) and

(2) to infer the LR.

Direct methods

For direct methods, we need to be able to evaluate f(d|Hss) and f(d|Hds) for any value of d. For this purpose, part of the available dataset can be used to build pairs of traces of the two types:

same-source and different-source pairs. If d is scalar, or of dimension 2 or 3 at most, the two empirical distributions can be fitted, with Gaussian mixtures, for example, leading to parametric estimates of f(d|Hss) and f(d|Hds) [START_REF] Mattijssen | Validity and Reliability of Forensic Firearm Examiners[END_REF][START_REF] Riva | Automatic comparison and evaluation of impressions left by a firearm on fired cartridge cases[END_REF][START_REF] Riva | Comparison and interpretation of impressed marks left by a firearm on cartridge cases -Towards an operational implementation of a likelihood ratio based technique[END_REF]. In the case of many features, a fit of each component of d can be performed in the same way, and the overall likelihoods can be approximated through the product of the likelihoods in each dimension, leading to the naïve Bayes classifier. To be successful, the latter approach requires however that the features are not overly correlated.

Indirect methods

The aim of these method is to build a discriminative model of the boundary between the two categories of pairs (same-source and different-source pairs) rather than a generative model explicitly parameterizing the distributions in the two categories. It is well known that, under the hypothesis of single-Gaussian distributions with the same variance under Hss and Hds in the scalar case, or same covariance matrix in the multidimensional case, the posterior probability of Hss takes the form of a sigmoidal curve [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF][START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF], hence the motivation for a logistic regression approach. Despite its result being a discriminative model, it enables to calculate the posterior probability of Equation (1) as well as the LR of Equation (2). As a matter of fact, the logistic regression model with parameters q = [a T b] T has output:

(3)

r d a d b ( , ) 1
1 exp ( ) T ( ) θ = + - +
where b is a scalar, and a is either a scalar in the case of a scalar distance, or otherwise an nvector (of the dimension of the evidence). If the proportions of the same-source and differentsource categories in the calibration set are denoted by fss and fds, r(d, q) approximates:

(4) Thus, the posterior probability for a priori probabilities P(Hss) and P(Hds) can be retrieved with:

(5) and the likelihood ratio with: (6)

Pros and cons

The indirect method offers several advantages:

-it spares the necessity to fit the likelihoods, -in the multi-dimensional case, contrary to the naïve Bayes classifier, the independence assumption is not necessary, because the logistic regression automatically takes care of the correlation between features, -by construction, the log LR is defined by a hyperplane, and thus robust with respect to the equal variance assumption, and to outliers or sparse data far from the boundary, -in the forensic context, since the log LR is directly proportional to a T d + b, the logistic allows a convenient and interpretable calibration of the dissimilarity score d, and a fusion of scores in the multidimensional case [START_REF] Morrison | Tutorial on logistic-regression calibration and fusion: converting a score to a likelihood ratio[END_REF].

On the other hand, the indirect method might suffer from:

-a reduced flexibility since it amounts to assuming single-Gaussian distributions, -a possibly important computation time in the case of high-dimensional evidence and of a distance of the same dimension.

These advantages and disadvantages will be examined and discussed on the example of hand odor evidence using a large panel of subjects.

f d H f f d H f f d
( ) = + - + LR d a d b f f ( ) exp T ds ss ( ) = +

Dataset description

A panel of 534 volunteers was set up which gathers 218 men and 316 women aged 7 to 94 years (median 28, interquartile interval [22 ; 48]), see Table 1 for the detailed composition in terms of gender and age. Note that this composition does not aim at reflecting that of a precise target population, such as one which is more likely to commit a crime, but to be as representative of the diversity of odors as possible. As a matter of fact, criminal investigations also often necessitate to look for victims, or to discriminate between traces from different people present at a crime scene, including those of victims or witnesses, who might be women as well as men, children or seniors as well as middle-aged adults. All data were completely anonymized prior to analysis, and no personal information was stored.

The goal was here to identify the subjects by their hand odor, whose volatile profile was shown to display a between-subject variability which is sufficient for differentiation [START_REF] Curran | The differentiation of the volatile organic signatures of individuals through SPME-GC/MS of characteristic human scent compounds[END_REF][START_REF] Cuzuel | Origin, analytical characterization and use of human odor in forensics[END_REF]. Also, in the forensic context, the hands have the advantage to be more likely to be directly in contact with objects at a crime scene, and to be easier to sample during a police interrogation.

The volatile profiles were obtained by a direct sampling procedure using identical sample collection kits of 4 small polymer strands that the subjects were asked to rub together in their hands for 15 minutes. The polymer strands were thermodesorbed, and the concentrated substances were separated by comprehensive bidimensional gas chromatography (GCxGC) coupled with mass spectrometry (MS). The sampling method and the optimization of the GCxGC-MS analysis were extensively described in [START_REF] Cuzuel | Sampling method development and optimization in view of human hand odor analysis by thermal desorption coupled with gas chromatography and mass spectrometry[END_REF][START_REF] Cuzuel | Human odor and forensics. Optimization of a comprehensive gas chromatography method based on orthogonality: how not to choose between criteria[END_REF]. Data were acquired, converted to .mzXML files with GC Real Time Analysis 4.20 (Shimadzu software), and then processed with MatlabTM (Natick, MA, USA) version 9.6.0.1150989 (R2019a), its Statistics and Machine

Learning Toolbox version 11.5 and its Bioinformatics Toolbox version 4.12.

Using a "home-made" Matlab script [START_REF] Cuzuel | Développement d'une stratégie de caractérisation chimique de la signature odorante d'individus par l'analyse chimiométrique de données issues de méthodes séparatives bidimensionnelles[END_REF], the preliminary manual processing of 25 chromatograms obtained on 3 subjects between 23 and 26 years old of both genders sampled several times at different time instants enabled us to draw up a first list of several hundreds of peaks. A library was built to store their retention times, their linear retention index, their mass spectrum, and the name of the corresponding compound when it could be identified using the NIST library. Indeed, if the availability of its mass spectrum is compulsory, a compound does not need to be formally identified for the comparison of chromatograms. We also checked whether compounds described in the literature as constituents of the human hand odor were present in this library, otherwise they were included. The library was then continuously enriched as the panel was increased with compounds potentially relevant to human hand odor because of their empirical frequency in new samples. This work led us to a customized library of 741 compounds, which were looked for in each chromatogram. As a result, each sample was characterized by the peak areas of 741 compounds.

In order to compensate for uncontrolled variations of the total area of the chromatograms, the sum of these areas was normalized in logarithmic scale to unit value, see Table 2 for a comparison of the reproducibility of the data without normalization, and normalization in scalar and logarithmic scales. Not knowing whether all 741 compounds are really relevant for identification, such a normalization across all compounds might be questionable. Thus, the possibility to avoid the problem by working on the dichotomized areas, i.e. 1 if the compound is present, or 0 if it is absent from the sample, was also investigated. Also, this approach might be of interest for forensic identification problems dealing with intrinsically binary features, such as gradient, structural and concavity (GSC) binary features in handwriting identification [START_REF] Srihari | On the Discriminability of the Handwriting of Twins[END_REF]. In the following, we refer to these traces of 741 features, continuous or dichotomized, as "odor traces".

As stated above, the subjects were sampled in quadruplicate, but due to unavoidable mishaps with some samples (like accidently dropping a polymer on the floor during sampling) and to chromatographic problems (such as failures of the cryogenic modulator), 1690 odor traces were obtained for the 534 subjects (44 were sampled once, 77 twice, 160 three times, and the remaining 253 subjects four times, leading to an average of 3.2 odor traces per subject).

This data set was split into a calibration set for training and validation, and an independent test set for performance estimation. Since gender [START_REF] Chang | Comparison and combination of ear and face images in appearance-based biometrics[END_REF] and age [START_REF] Haze | 2Nonenal newly found in human body odor tends to increase with aging[END_REF] are known to impact odor traces, the split was made so as to respect the gender proportions, with subject of all ages in the two sets, and odor traces of the same subject being put in the same set. As a result, the calibration set comprises 412 subjects and their 1 299 odor traces (corresponding to 1 594 Hss and 841 457 Hds pairs), and the test set comprises the remaining 122 subjects and their 391 odor traces (leading to 481 Hss and 75 764 Hds pairs). The way the odor traces distribute between calibration and test set can be grasped through the Principal Component Analysis (PCA) of Figure 1.

Implementation

Three methods are implemented:

1) the direct method using a scalar distance between odor traces, 2) the indirect method using a scalar distance between odor traces, 3) the indirect method using a vectorial distance, i.e. a distance on each odor compound.

Note that, given the large dimension of the problem (n=741) and the known correlations between features, we did not attempt to implement the direct method using a vectorial distance (i.e. a scalar distance on each feature and the naïve Bayes classifier).

Distances between two odor traces

Let xi denote the n-vector representing odor trace i. A standard choice of scalar dissimilarity measure between odor traces i and j is the Euclidian distance between vectors xi and xj, i.e.:

But this distance is not robust with respect to shifts and linear transformations of the features.

Thus, a more appropriate distance would be Pearson's linear correlation based distance:

However, Pearson's correlation is sensitive to non-linearities, whereas Spearman's nonparametric correlation coefficient on the ranks is able to capture monotonic nonlinear associations as well as linear ones [START_REF] Daniels | The relation between measures of correlation in the universe of sample permutation[END_REF]. Thus, the Spearman correlation based distance (same as Equation ( 8) with the xi k and xj k replaced by their ranks) is expected to be more robust with respect to nonlinear variations of the peak areas. The three distances were compared in a previous study where Spearman's correlation based distance indeed clearly outperformed the two other distances [START_REF] Cuzuel | Human odor and forensics: towards Bayesian suspect identification using GCxGC-MS characterization of hand odor[END_REF]. Therefore, the direct and indirect method using a scalar distance have been implemented with Spearman's correlation based distance.

The chosen vectorial distance for the indirect method using a vectorial distance is simply the vector of the absolute differences between feature values:

(9) ∑ = - = d x x x x ( , ) ( ) Euclid i j i k j k k n 2 1 ∑ ∑ ∑ = - - - - - = = = d x x x x x x x x x x ( , ) 1 ( )( ) ( ) ( )
Pearson i j i k i j k j k n i k i k n j k j k n 1 2 1 2 1 = - - d x x x x x x ( , ) | |...| | vectorial i j i j i n j n T 1 1

Estimation of the likelihoods for the direct method

The calibration set was used to build pairs of odor traces of same and different sources, and to compute their distances. The empirical densities being essentially uni-or bimodal (see Figure 2), they were best fitted with a two-Gaussian mixture distribution, using Matlab's function "fitgmdist", leading to estimates of the likelihoods f(d|Hss) and f(d|Hds).

When performing feature selection (see section 2.5.5 below), the same-source and differentsource distributions vary, depending on the subset of features. However, checking the optimality of the two-Gaussian mixture for all subsets of features during the feature selection process would be too time consuming. Nevertheless, the convergence of the fits was checked for, each fit being repeated three times with a different random set of initial parameters in order to retain the fit with the largest likelihood, and the optimality of the fit obtained with the selected number of features was tested. Of course, as stated in section 2.3, the necessity to fit the likelihoods and to optimize these fits, or to suffer suboptimal fits, is the first disadvantage of the direct method.

Estimation of the logistic model for the indirect methods

The logistic regression model of Equation ( 3) with parameters q = [a T b] T was fitted through maximum likelihood, by minimizing the cross-entropy cost function with Matlab's function "glmfit".

Likelihood ratio and performance estimation

For the direct method, the LR was evaluated using the estimates of the likelihoods f(d|Hss) and f(d|Hds) and Equation ( 2), as a function of the distance d. For the indirect methods, the LR was obtained from the fitted logistic regression and Equation [START_REF] Tang | Computational methods for the analysis of footwear impression evidence[END_REF], and plotted as a function of the distance d or of the score a T d + b, depending on d being scalar of vectorial.

Since there is no true reference for the LR, the performance of the different methods was evaluated by estimating the posterior probability P(Hss|d) according to Equations ( 1) and ( 5)

for the direct and indirect methods respectively, and by performing a binary classification using equal prior probabilities (P(Hss) = P(Hds) = 0.5). Varying the decision threshold on P(Hss|d), the sensitivity and the specificity were estimated on the calibration and test sets, and used to compute the corresponding areas under the receiver operating characteristic "ROC" curve (AUC) [START_REF] Hanley | The Meaning and Use of the Area under a Receiver Operating Characteristic (ROC) Curve[END_REF]. The performance was further characterized by the sensitivity Sn and specificity Sp maximizing Youden's index [START_REF] Youden | Index for rating diagnostic tests[END_REF], i.e. Sn + Sp -1.

Feature selection

In a previous study [START_REF] Cuzuel | Human odor and forensics: towards Bayesian suspect identification using GCxGC-MS characterization of hand odor[END_REF], improved results were obtained using feature selection. The idea is to retain the features that contribute the most to the difference between distance densities under Hss and Hds. Given the large number of features (odor compounds) and the large size of the data set, a filter approach to this selection was chosen. Filter approaches are based on a statistical measure of the difference between the two densities for each feature. They are hence independent from the main algorithm (direct or indirect method), as opposed to the so-called wrapper approaches, which sequentially evaluate the relevance of each feature subset based on the performance of the whole procedure, i.e. feature selection together with main algorithm [START_REF] Venkatesh | A Review of Feature Selection and Its Methods[END_REF]. Since filter approaches consider the features independently, they might retain redundant features, but in turn this gives them more robustness and, most importantly, they require less computation time. For each feature, considering the absolute values of the differences for the Hss and Hds pairs, we chose Wilcoxon's non-parametric test statistic as statistical measure in the case of continuous features, and Fisher's exact test statistic in the case of dichotomized features. Then, the features were ranked in decreasing order of the onesided p-value of the test, which is a one-sided test since smaller distances between features under Hss than under Hds are sought for. The number of features maximizing the AUC was estimated on the calibration set using 3-fold cross-validation. The cross-validation partitions were randomly chosen with the constraint that the odor traces of the same subject were put in the same partition. Note that cross-validation also enabled us to estimate the uncertainty on the AUCs through the mean standard deviation on the three partitions.

Results and discussion

The three methods are first evaluated using all the features of the odor traces (baseline comparison) and then, the possibility to further improve their performance using feature selection is investigated.

Baseline comparison of the three methods (without feature selection)

The results obtained with the three methods on the calibration and test sets are summarized in Tables 3 and4 for dichotomized and continuous features respectively. As a first remark, the performance of the direct and indirect methods using a scalar distance in terms of AUC and of sensitivity and specificity are almost identical, for both dichotomized and continuous features.

Thus, the higher flexibility of the direct method does not increase the performance. On the contrary, its lack of robustness can be visualized on Figure 2 depicting the posterior probability and the LR obtained with the dichotomized features: due to the larger variance of the likelihood under Hss, the posterior probability and the LR, instead of being monotonous, start to increase with the distance at some point (d » 0.7). Whereas whatever the situation with the indirect method, posterior probability and LR always decrease with d, see Figure 3 depicting the posterior probability and the LR obtained with indirect method, this time on the continuous features.

Also noteworthy, the performance obtained with the indirect method using a vectorial distance is significantly better than those of the methods working with a scalar distance: the AUC on the calibration and test sets jumps from 91-92% to 97-98%, the standard deviation of the AUC being estimated at 0.7% using 3-fold cross-validation on the calibration set. The distributions of the score (a T d + b) resulting from the logistic regression, the regression itself, the posterior probability and the LR are shown in Figure 4. The only drawback lies in the increased, but perfectly tractable computational cost (10 minutes instead of a few seconds, on a 4,2 GHz Intel Core i7).

Finally, the binarization of the features decreases the performance, but only marginally (the AUC is decreased by » 1%). Note that, in this precise case where the features quantify the amount of odor compounds, this could be due to the fact that the normalization of the compound proportion uses all these compounds whereas it is not known whether they are all relevant. Note also that the normalization was improved by performing it in the logarithmic scale rather than in the linear scale (the former improving the reproducibility, see Table 2), with which continuous features did not outperform dichotomized features, as shown in a

previous study [START_REF] Cuzuel | Human odor and forensics: towards Bayesian suspect identification using GCxGC-MS characterization of hand odor[END_REF]. Finally, other normalization methods specific to GCxCG-MS data might advantageously be investigated [START_REF] Chen | Influences of Normalization Method on Biomarker Discovery in Gas Chromatography-Mass Spectrometry-Based Untargeted Metabolomics: What Should Be Considered?[END_REF]], but are outside the scope of this paper.

Comparison of three methods with feature selection

The number of selected features using the filter approach is reported in Tables 5 and6, together with the corresponding results on the calibration and test sets, for dichotomized and continuous feature respectively.

Again, there is almost no difference in performance between the direct and indirect methods with a scalar feature, be it on dichotomized or continuous features. In terms of AUC, the selection is more efficient on binary features than on continuous ones (94.4% with selection instead of 91.5% without for dichotomized features, 93.5% instead of 93.0% for continuous features, on the test set), with an important reduction of the number of the dichotomized features (267 instead of 741), and a moderate one for continuous features (535 instead of 741). Note also that in both cases, this increased performance benefits the specificity, which is highly desirable in a forensic application (it is crucial in this context not to reject the different-source hypothesis, i.e. the defense hypothesis, when in fact it is true).

For both dichotomized and continuous features, the indirect method using a vectorial distance is again significantly better than the previous ones (AUCs around 97-98% instead of 94-95% on both calibration and test sets), with a similar number of selected features (440 for dichotomized features, 500 for continuous ones). However, in both cases, the parsimony due to feature selection does not increase the performance as compared to the baseline method,

it is quasi-identical with and without selection. In return, this testifies to a robustness of the indirect method with respect to possibly irrelevant features. And of course, not to have to perform the selection spares computation time.

Discussion of the choice of equal priors

In this manuscript, the methods are compared in terms of AUC, sensitivity and specificity. In the case of the indirect methods, the regression being obtained by fitting a logistic function to the data, whatever the prior probabilities P(Hss) and P(Hds), the posterior probability P(Hss|d)

given by Equation ( 5) is also a logistic function. Thus, when the threshold on P(Hss|d) is varied from 1 to 0, the same ROC curve is described, whose AUC only depends on the distance distributions under Hss and Hds: only the threshold maximizing Youden's index changes.

With the direct method, the choice of the prior has an influence on the shape of the posterior probability, so that the threshold on P(Hss|d) can possibly be varied in a different interval (see Figure 3 where P(Hss|d) never reaches 0 for example). However, in practice, there is no influence on AUC, sensitivity and specificity because, again, the AUC depends essentially on the distance distributions under Hss and Hds. The only palpable change is on the threshold yielding the best compromise between sensitivity and specificity, threshold which adjusts to P(Hss) by roughly following it.

Thus, the assumption of equal prior probabilities has practically no impact on the LR estimate.

Limitations

From a practical point of view, our work suffers several limitations for a real-world forensic application. First, for practical reasons, the subjects were sampled at a single time point, so that the variability of the data is essentially due to the analytical variability. Second, the chromatograms being compared are of the same nature, i.e. obtained on samples provided by directly sampling the subjects (with contact with the subjects' whereas in real life, the unknown source sample will be obtained indirectly from an object on crime scene (without contact with the subject). Third, the odor collected on the crime scene might be contaminated by other odors, from the environment or from other people. A study focused on mixtures of odors, contaminations, and weathered traces has not been carried out yet but is considered. However, despite these controlled conditions, the PCA of Figure 1 and the statistics of Table 2 show that the data is already of limited reproducibility, so that the good results we have obtained are encouraging concerning the robustness of the best method to more realistic sampling conditions.

From a methodological point of view, the proposed methods are based on a common source scenario, where it is asked whether the two traces originate from the same source or from different sources without specifying which sources are considered, and not on a specific source scenario, where the question is whether the two traces stem specifically from the known source [START_REF] Neumann | Defence Against the Modern Arts: the Curse of Statistics -Part II: "Score-based likelihood ratios[END_REF]. The problem of the common source scenario is that it does not take account of the typicality of the source, contrary to recommendations for a better estimation of the strength of evidence through the LR [START_REF] Morrison | Tutorial on logistic-regression calibration and fusion: converting a score to a likelihood ratio[END_REF][START_REF] Tang | Likelihood ratio estimation in forensic identification using similarity and rarity[END_REF][START_REF] Morrison | Score based procedures for the calculation of forensic likelihood ratios -Scores should take account of both similarity and typicality[END_REF]. But to implement a specific source scenario, a number of traces from the known source are needed in order to be able to estimate the distribution under Hss (for the direct method) or to discriminate between the Hss and Hds populations (for the indirect methods), which is quite unpractical when dealing with human hand odor, and not feasible at this stage of the study (at most four usable odor traces were obtained, for only 253 subjects among the 534).

Conclusion

To summarize, the advantages expected from an indirect method are fully obtained, in particular the dispensation to parameterize the likelihoods, and the robustness with respect to differences in their variance and/or to possible outliers. Moreover, an increase in performance of the indirect method as compared with the direct one is not obtained with a scalar distance between odor traces, but when using the vector of the distances between each feature of the odor traces. This improvement was not really expected, because, especially in the forensic context, it is often advocated to convert multivariate data to a univariate datum summarizing the relationship between features. Finally, the indirect method with a vectorial distance proves also robust with respect to potentially irrelevant features since removing them does not modify the performance, an appealing quality for dealing with traces which are not yet solidly characterized, such as odor traces.

FIGURE CAPTIONS Figure 1 .

 CAPTIONS1 FIGURE CAPTIONS

Figure 2 .

 2 Figure 2. Results of the direct method with the scalar distance d on dichotomized feature, as functions of d, on the calibration set. a) Empirical distribution of d under Hss (1 594 pairs), and estimated density (mixture of 2 Gaussians). b) Empirical distribution of d under Hds (841 457pairs), and estimated density (mixture of 2 Gaussians). c) Posterior probability of Hss obtained using Bayes' formula with equal priors. d) Likelihood ratio.

Figure 3 .

 3 Figure 3. Results of the indirect method with the scalar distance d on continuous features, as functions of d, on the calibration set. a) Empirical distribution of d under Hss (1 594 pairs). b) Empirical distribution of d under Hds (841 457pairs). c) Logistic regression (dotted line), and deduced posterior probability of Hss with equal priors (continuous line). d) "Likelihood ratio" exp(a T d +b) corresponding to the logistic regression (dotted line), and likelihood ratio (continuous line).

Figure 4 .

 4 Figure 4. Results of the indirect method with the vectorial distance d on continuous features, as functions of the score a T d +b, on the calibration set. a) Empirical distribution of the score a T d + b under Hss (1 594 pairs). b) Empirical distribution of the score a T d + b under Hds (841 457pairs). c) Logistic regression (dotted line), and deduced posterior probability of Hss with equal priors (continuous line). d) "Likelihood ratio" exp(a T d +b) corresponding to the logistic regression (dotted line), and likelihood ratio (continuous line).

Figure 1 .

 1 Figure 1. PCA of the data set showing the distribution of the odor traces between calibration and test sets. The PCA was performed on the covariance matrix of the continuous features, i.e. normalized in logarithmic scale.

Figure 2 .

 2 Figure 2. Results of the direct method with the scalar distance d on dichotomized feature, as functions of d, on the calibration set. a) Empirical distribution of d under Hss (1 594 pairs), and estimated density (mixture of 2 Gaussians). b) Empirical distribution of d under Hds (841 457pairs), and estimated density (mixture of 2 Gaussians). c) Posterior probability of Hss obtained using Bayes' formula with equal priors. d) Likelihood ratio.

Figure 3 .

 3 Figure 3. Results of the indirect method with the scalar distance d on continuous features, as functions of d, on the calibration set. a) Empirical distribution of d under Hss (1 594 pairs). b) Empirical distribution of d under Hds (841 457pairs). c) Logistic regression (dotted line), and deduced posterior probability of Hss with equal priors (continuous line). d) "Likelihood ratio" exp(a T d +b) corresponding to the logistic regression (dotted line), and likelihood ratio (continuous line).

Figure 4 .

 4 Figure 4. Results of the indirect method with the vectorial distance d on continuous features, as functions of the score a T d +b, on the calibration set. a) Empirical distribution of the score a T d + b under Hss (1 594 pairs). b) Empirical distribution of the score a T d + b under Hds (841 457pairs). c) Logistic regression (dotted line), and deduced posterior probability of Hss with equal priors (continuous line). d) "Likelihood ratio" exp(a T d +b) corresponding to the logistic regression (dotted line), and likelihood ratio (continuous line).

Table 1 .

 1 Panel composition in terms of gender and age.

	age	7-17	18-64	65-94	total
	gender
	man	18	182	18	218
	woman	28	268	20	316
	total	46	450	38	534

Table 2 .

 2 Repeatability of the odor trace features (peak areas), without and with two normalizations, estimated on the 490 subjects sampled at least twice (IQI stands for interquartile interval).

	Normalization	Median relative standard deviation [IQI] in %
	None	61.0 [32.1 ; 78.0]
	In scalar scale	56.1 [31.9 ; 74.2]
	In logarithmic scale	33.2 [17.9 ; 48.8]

Table 5 .

 5 Comparison between the three methods on the calibration and test sets, using dichotomized features, in terms of AUC, sensitivity (Sn) and specificity (Sp) maximizing Youden's index, with selection of the number of features among the 741 by cross-validation of the calibration set (AUC-CV3 the mean 3-fold cross-validation AUC on the calibration set, and #feat. denotes the number of features selected by cross-validation).

	Calibration	Test

Table 6 .

 6 Comparison between the three methods on the calibration and test sets, using continuous features, in terms of AUC, sensitivity (Sn) and specificity (Sp) maximizing Youden's index, with selection of the number of features among the 741 by cross-validation of the calibration set (AUC-CV3 the mean 3-fold cross-validation AUC on the calibration set, and #feat. denotes the number of features selected by cross-validation).

	Calibration	Test


Sophie Cornu
email: sophie.cornu@inrae.fr

Jean-Dominique Meunier

Céline Ratie

Fréderic Ouedraogo

Yves Lucas

Patricia Merdy

Doris Barboni

Camille Delvigne

Daniel Borschneck

Olivier Grauby

Catherine Keller

Allophanes, a significant soil pool of silicon for plants

Keywords:

Introduction

The variability of Si content in soils is large, from <1 % wt in Histosols to greater than 45 % wt in some Podzols [START_REF] Sommer | Silicon pools and fluxes in soils and landscapes-a review[END_REF]. Si is included in a large number of soil solid phases of various origins: (1) minerals derived from the parent material as quartz, feldspar, phyllosilicates (e.g. muscovite, biotite) etc.; (2) solid phases resulting from the transformation by weathering and pedogenesis of the primary minerals into clay minerals (kaolinite, smectite, illite among others), crystalline and amorphous material, including poorly crystalline Fe oxides and short-range ordered (SRO) aluminosilicates (e.g. allophanes and imogolites); (3) amorphous silica particles (BSi) formed by Si biocycling and synthesized by plants (phytoliths), animals (sponge spicules) or protists (testate amoeba shells, diatom frustules) [START_REF] Sommer | Silicon pools and fluxes in soils and landscapes-a review[END_REF][START_REF] Puppe | The protozoic Si pool in temperate forest ecosystems-quantification, abiotic controls and interactions with earthworms[END_REF]. While soils offer a wide range of Si sources, plants absorb Si only from the soil solution as Si(OH)4. Therefore, deciphering the source of soluble Si (DSi) is challenging because of the complex biogeochemical interactions involving Si in soils [START_REF] Schaller | Silicon Cycling in Soils Revisited[END_REF]. The conversion of forest to cropland was shown to lead to several Si pool changes: (1) a decrease in the phytolith pool [START_REF] Conley | Deforestation causes increased dissolved silicate losses in Hubbard Brook Experimental Forest[END_REF][START_REF] Struyf | Historical land use change has lowered terrestrial silica mobilization[END_REF][START_REF] Guntzer | Long-term removal of wheat straw decreases soil amorphous silica at Broadbalk, Rothamsted[END_REF][START_REF] Vandevenne | Silicon pools in human impacted soils of temperate zones[END_REF], in the easily soluble Si pool [START_REF] Clymans | Anthropogenic impact on amorphous silica pools in temperate soils[END_REF] and the weathering of silicate minerals [START_REF] Struyf | Historical land use change has lowered terrestrial silica mobilization[END_REF] and (2) an increase in the pedogenic Si fraction i.e. clays [START_REF] Vandevenne | Silicon pools in human impacted soils of temperate zones[END_REF]. In addition, the growing crops such as rice, sugar cane, or wheat [START_REF] Hodson | Phylogenetic variation in the silicon composition of plants[END_REF][START_REF] Guntzer | Long-term removal of wheat straw decreases soil amorphous silica at Broadbalk, Rothamsted[END_REF][START_REF] Katz | Silica phytoliths in angiosperms: Phylogeny and early evolutionary history[END_REF] that are among the most important Si accumulators (1 % or more Si by dry weight) are predicted to change the global Si biogeochemistry in the future by increasing the reservoir of BSi on earth [START_REF] Carey | Human appropriation of biogenic silicon-the increasing role of agriculture[END_REF].

Various methods have been used to estimate DSi and bioavailable Si [START_REF] Liang | Silicon in Agriculture[END_REF]. Their results suggest that bioavailable Si concentration is both negatively correlated with total Si content, reflecting the predominance of low solubility minerals such as quartz in the soils [START_REF] Yanai | Evaluation of available silicon content and its determining factors of agricultural soils in Japan[END_REF][START_REF] Landr´e | Do climate and land use affect the pool of total silicon concentration? A digital soil mapping approach of French topsoils[END_REF] and positively correlated with soil properties such as weatherable aluminosilicates (Li et al., 2020a), phytolith content (Li et al., 2020b;[START_REF] De Tombeur | Plants sustain the terrestrial silicon cycle during ecosystem retrogression[END_REF], pH [START_REF] Klotzbücher | Impact of agriculture practices on plant-available silicon[END_REF][START_REF] Meunier | pH as a proxy for estimating plant-available Si? A case study in rice fields in Karnataka (South India)[END_REF][START_REF] Miles | Extractable Silicon in Soils of the South African Sugar Industry and Relationships with Crop Uptake[END_REF][START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF], soil organic matter, iron oxides [START_REF] Yanai | Evaluation of available silicon content and its determining factors of agricultural soils in Japan[END_REF] and the clay (<2 μm) fraction contents [START_REF] Liang | Silicon in Agriculture[END_REF][START_REF] Yanai | Evaluation of available silicon content and its determining factors of agricultural soils in Japan[END_REF][START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF].

In temperate soils, despite the observed decrease of soil phytolith content under agriculture, higher DSi (estimated by the CaCl2 method referred to herein as SiCaCl2) concentration was observed under agriculture compared to forest developed on sedimentary rocks [START_REF] Vandevenne | Silicon pools in human impacted soils of temperate zones[END_REF][START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF]. This was related to a higher pH observed in agricultural soils compared to forest soils [START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF] pointing out two mechanisms: greater phytoliths dissolution, or greater adsorption onto clay minerals at higher pH. [START_REF] Cornu | Evidence of short-term clay evolution in soils under human impact[END_REF] also observed that the liming of agricultural soils in Europe induced a change in the nature of the clay minerals, and hence, potentially in the relative proportion of the different Si pools. Consequently, despite their lower solubility compared to phytoliths [START_REF] Fraysse | Surface chemistry and reactivity of plant phytoliths in aqueous solutions[END_REF], clays may play an important role in maintaining high values of SiCaCl2, somehow cancelling out the phytolith pool decrease [START_REF] Vandevenne | Silicon pools in human impacted soils of temperate zones[END_REF][START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF].

To identify the mechanisms controlling SiCaCl2 in temperate soils and to understand the impact of land use on these mechanisms in different soil types, we designed a paired-site approach in which we sampled adjacent soils differing only in their land use while other pedological factors (topography, parent material, climate, and age) were the same. Seven paired sites were selected with calcaric and non-calcaric soils, developed on sedimentary parent materials (limestone and loess). Beside the SiCaCl2 pool, we measured the pools of Si associated to crystalline and amorphous material, including poorly crystalline Fe oxides and SRO-aluminosilicates, and biogenic Si particles. Only the topsoil horizons (0-0.3 m) were analysed as they were considered the most impacted by land use change.

Material and methods

Identification of the paired sites and sampling strategy

Paired sites are not easy to locate as they have to differ only by their land use, all other pedological factors (parent material, climate, topography, and time) being kept constant between the two sites [START_REF] Filippi | Chapter Four -Monitoring and Modeling Soil Change: The Influence of Human Activity and Climatic Shifts on Aspects of Soil Spatiotemporally[END_REF]. We chose soils on sedimentary materials as land use impacts SiCaCl2 on these soils, and calcaric and non-calcaric soils as their behaviour with respect to SiCaCl2 differs [START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF]. In France, most cultivated soils with these characteristics are developed either on limestone or on decarbonised loess. We also chose to study soils with different degrees of pedogenetic evolution: calcaric and hypereutric Cambisols for limestone; Luvisols and Albeluvisol for loess [START_REF] Wrb | World reference base for soil resources[END_REF]. To identify the paired sites, we used the French Soil Quality Network (RMQS) of 2111 georeferenced sites in France spread over a regular grid (16*16 km) for which soil type, parent material and land use are provided [START_REF] Jolivet | Manuel du R´eseau de Mesures de la Qualit´e des Sols[END_REF]. All the RMQS sites with the appropriate parent material and soil type either cultivated or under forest, were extracted from the database. For each selected RMQS site, an adjacent site (i.e. with the same climate), with the same topographical position (therefore expected to be on the same type of soil and parent material) and under the complementary land use (either forest or cultivation) was identified on the basis of: (i) the maps of the French National Institute of Geographic and Forest Information (IGN, 1:25000); (ii) Cassini's maps, the oldest available maps with records of land use [START_REF] Cassini | Composite: Carte de France[END_REF], both available on Geoportail (https://www.geoportail.gouv.fr/); (iii) aerial photographs taken at different times during the 20th century. Only paired sites for which the land use change occurred more than 100 years ago were kept. In most cases, the forested areas have been preserved from cultivation because they belonged to the state or to castle parks. Nineteen pairs of adjacent plots were identified. Another paired site of Albeluvisol [START_REF] Wrb | World reference base for soil resources[END_REF] on already well-characterised loess was added [START_REF] Montagne | Changes in the pathway and the intensity of albic material genesis: Role of agricultural practices[END_REF]. The 20 sites were in the northern half of France (Table SI1), where the cultivation of wheat is common. Most of them included cereals in their rotation (Table SI1). Forests were mixed deciduous forests except for site 62 that had a conifer forest (Table SI1). All the selected paired sites were auger-sampled and described in the spring of 2015. In the cultivated area, one sample was taken in the upper horizon which is homogenized by ploughing (usually 0-25 cm). In the nearby forested area two samples were taken (usually 0-5 and 5-25 cm). On the basis of the soil morphology, eighteen out of the 20 paired sites were kept and analysed for rare earth element (REE) concentrations and particle size distribution to detect potential changes in parent materials. Analyses were performed at the Laboratory of Soil Analysis (Arras, France). Seven paired sites were finally kept (Table SI1), as their REE and grain size skeleton patterns did not differ between land uses considering the analytical uncertainties (see supplementary information, SI, part 3 for more details) and thus could be considered as having the same parent material. Bulk density measurements were performed using the volumetric cylinder method in cultivated areas and the excavation method in forested areas because the root density prevented the use of the cylinder. When bulk density could not be measured (e.g., at site 1166), it was estimated using the pedotransfer function proposed by [START_REF] Alexander | Bulk Densities of California Soils in Relation to Other Soil Properties[END_REF], which relates bulk density to soil organic carbon (SOC) content. To compare the forested and cultivated sites, composite values (0-25 cm) were calculated for the forested sites considering the thickness of the two sampled horizons (0-5 cm and 5-25 cm) and their respective bulk densities. 2.2. Samples characterisation 2.2.1. Pedological analysis Air dried samples were analysed for pedological characteristics at the Soil Analysis Laboratory of INRAE (Arras, France). The particle-size distribution was obtained by wet sieving and the pipette method (NF X 31-107). Total SOC concentration was measured by dry combustion (NF ISO 10694), CEC (cation exchange capacity) by the cobaltihexamine chlorine method (NF X 31-130), and pH in water with a 1:5 soil:water ratio (NF X 31-107). The particle size fraction contents are hereafter denominated < 2 μm, 2-20 μm, 20-50 μm, 50-200 μm, 200-2000 μm. The < 2 μm fraction CEC (CEC<2μm) was calculated following equation ( 1).

CEC<2μm = ((CEC CECSOC* SOC) / < 2μm)*1000 (1)
where SOC and < 2 μm are the contents in SOC and in the < 2 μm fraction, respectively. We used 200 cmol+ kg 1 as an estimate of CECSOC [START_REF] Baize | Guide des analyses en p´edologie[END_REF] of both the ploughed horizon of the cultivated site and for the 5-25 cm horizon of the forested plots and 100 cmol+ kg 1 as an estimate of CECSOC of the 0-5 cm horizon of the forested plots that contains a greater fraction of particulate organic matter with a lower CEC (Balesdent, pers. comm., 2018).

Total major elements concentrations in soil, including Si (Sitot), were also analysed at the Soil Analysis Laboratory of INRAE (Arras, France) by ICP-MS analysis after HF tri-acid dissolution (Altot, Fetot, Natot, Catot, Mgtot, Ktot, Mntot) or alkaline fusion (Sitot, NF ISO 14869-2).

Chemical extraction of different soil Si pools

We used three types of soft chemical extraction in parallel to access different soil Si pools:

(1) Soluble and bioavailable Si was extracted using the CaCl2 method (SiCaCl2; [START_REF] Haysom | Some aspects of the calcium silicate trials at Mackay[END_REF][START_REF] Narayanaswamy | Evaluation of selected extractants for plant--Available Silicon in Rice Soils of Southern India[END_REF]. The soils were mixed with a 0.01 mol L-1 CaCl2 solution with a solid/solution ratio of 10, in two replicates. The suspension was stirred for 16 h and filtered at 0.45 μm.

(2) Poorly crystalline Fe oxides and short-range ordered (SRO) aluminosilicates such as allophanes and imogolites were extracted by the oxalate method (Sio, [START_REF] Tamm | Eine Method zur Bestimmung der anorganishen Komponenten des Golkomplex in Boden[END_REF][START_REF] Wada | Allophane and imogolite[END_REF][START_REF] Harsh | Encyclopedia of Soils in the Environment[END_REF]. It is assumed that this method does not dissolve amorphous silica particles (i.e., phytoliths) [START_REF] Kodama | Tiron dissolution method used to remove and characterize inorganic components in soils[END_REF][START_REF] Saccone | Assessing the extraction and quantification of amorphous silica in soils of forest and grassland ecosystems[END_REF]. The Sio pool thus represents Si associated with amorphous co-precipitated Fe and Al compounds. Allophane content was estimated using [START_REF] Harsh | Encyclopedia of Soils in the Environment[END_REF] equation ( 2):

% allophane = Sio/( 0.067 Alo/Sio + 0.27) (2)
(3) Iron oxides were extracted by the dithionite-citrate-bicarbonate (DCB) method (Sid, [START_REF] Mehra | Iron Oxide Removal from Soils and Clays by a Dithionite-Citrate System Buffered with Sodium Bicarbonate[END_REF], which also dissolves the poorly crystalline fraction and may dissolve phytoliths [START_REF] Georgiadis | Development of a method for sequential Si extraction from soils[END_REF]. Si extracted by these three methods and Fe and Al extracted by the latter two were measured using ICP-OES (Inductively Coupled Plasma Optical Emission Spectroscopy). Both Sio and Sid were extracted and measured at the Soil Analysis Laboratory of INRAE (Arras, France).

The different soil Si forms were then estimated considering that (i) the DSi equals the SiCaCl2; (ii) Si in SRO minerals (allophanes, imogolites and poorly crystalline Fe oxides) equals (Sio -SiCaCl2); (iii) Si associated to crystalline Fe oxides or to the part of phytoliths dissolved by the extractant equals (Sid -Sio); and (iv) Si in crystalline silicates plus part of the phytoliths equals (Sitot -Sid). We have to keep in mind that all these extractions are operationally defined and may extract part of soil solid phases other than the targeted ones as already extensively discussed in the literature (e.g. [START_REF] Georgiadis | Development of a method for sequential Si extraction from soils[END_REF].

Error bars were estimated considering 5 % uncertainty on the chemical and bulk density measurements.

Phytolith extraction and characterisation

Soil phytoliths were quantified, identified and their weathering state was evaluated.

BSi particles were extracted from the soils using wet digestion and heavy liquid separation [START_REF] Jouquet | The distribution of silicon in soils influenced by termite bioturbation in South Indian forest soils[END_REF]; see SI part 2 for more details). No replicates have been done. Extracts, which happened to include phytoliths as well as low density particles such as micro-charcoals and amorphous silica skeletons of chrysophytes, diatoms, sponges, etc., were dried and weighed before being mounted on microscope slides to evaluate the real proportion of phytoliths versus non-phytolith particles in the extracts and, therefore, provide an estimation of the phytolith contents.

We considered that phytoliths contained 42 % of Si [START_REF] Alexandre | Plant impact on the biogeochemical cycle of silicon and related weathering processes[END_REF] and could thus estimate the soil phytolith Si pool based on the phytolith abundances in the soils.

Since cultivated plants produce morphotypes of phytoliths different from forest species (trees, shrubs, ferns, forbs), phytolith taxonomic identification was performed according to the International Committee for Phytolith Taxonomy (2019) and counting was carried out.

Lastly, the analysis of the phytolith surface is a good indicator of their weathering state [START_REF] Sommer | Silicon pools and fluxes in soils and landscapes-a review[END_REF][START_REF] Jouquet | The distribution of silicon in soils influenced by termite bioturbation in South Indian forest soils[END_REF], we measured the proportion of corroded vs non-corroded phytoliths (for more details, see SI part 2).

Mineralogical characterisation of the < 2 μm fraction

Dried soils were dispersed in 200 to 300 mL of ultrapure water by rotative shaking for 22 h using plastic bottles and 20 agate balls (Ø 5 mm) [START_REF] Balesdent | Effet des ultrasons sur la distribution granulom´etrique des mati`eres organiques des sols[END_REF]. The mass of dried soil was chosen to obtain ca. 4 g of < 2 μm fraction. The suspension thus obtained was further dispersed by addition of 2 mL of Na-hexametaphosphate at 150 g L-1, mechanical shaking, and 5-min ultrasonic dispersion. The < 2 μm fraction was then sampled according to the Stokes law and the settled soil particles resuspended for a new extraction step. The process was repeated until more than 80 % of the total fraction was recovered (that is more than 3.2 g of < 2 μm fraction out of 4 g). An aliquot of the < 2 μm suspension was dried to determine the < 2 μm suspension concentration. Oriented slides were made by deposition on glass slides of 10 mg of < 2 μm particles spiked with 10 % of boehmite added as a reference for X-ray diffraction (XRD) studies. The glass slides were then analysed with a Panalytical X'Pert Pro X-ray diffractometer equipped with a cobalt anti-cathode (λ = 1.79 Å) at 40 kV and 40 mA, with three treatments following [START_REF] Robert | M´ethode de pr´eparation des argiles des sols pour des ´etudes min´eralogiques[END_REF]: air-dried (AD) and ethylene glycol (EG) solvation and heating at 490 •C for 4 h.

Peak areas of the main mineral phases and boehmite were measured using the Fityk software [START_REF] Wojdyr | Fityk: a general-purpose peak fitting program[END_REF]. The peak area ratio with respect to boehmite was used to compare samples. Transmission electron microscope (TEM, JEOL JEM 2011) coupled to an energy-dispersive X-ray spectrometer (EDX, X-Flash Silicon Drift Detector 5030, Bruker) analysis was performed on the < 2μ m fraction of sample BS3F1 (Table SI1) for checking the presence of allophanes.

Data analysis

Pedological characteristics and major element concentrations of the seven studied paired sites were analysed by principal component analysis (PCA). To identify different groups of soils from a pedological and chemical point of view, a hierarchical classification was performed on the first two factors of the PCA. Soil characteristics of the paired sites were considered as different when the difference for a given soil characteristic exceeded the analytical uncertainties. As only one measurement by site was performed, no statistical tests could be performed to further test the statistical significance of this difference.

Significance of Pearson correlations between the SiCaCl2 concentrations and different soil parameters was estimated based on the p-values considering a threshold level of 5 %. All the statistical analyses were performed using the XLSTAT software for windows (Addinsoft).

Results

Seven paired sites with contrasted pedological characteristics and total concentrations of major elements

Pedological characteristics and total contents of major elements in the topsoils of the seven paired sites are reported in Tables SI1 andSI2. They were computed by PCA. The first two components explained 72 % of the dataset variance with 60 % for the first component (Fig. 1). The first component was correlated positively to Catot, the < 2 μm fraction, Fetot, Altot, Mgtot and the SOC contents, as well to pH and the CEC, and negatively to Sitot, Natot, silt (2-20 μm and 20-50 μm fractions) and fine sand (50-200 μm fraction) contents. The second component (12 % of the variance) was correlated positively to the 2-20 μm fraction, and negatively to the 50-200 μm fraction (Fig. 1a). The PCA discriminated three groups of soils confirmed by a hierarchical clustering (Fig. 1b). The first group corresponds to calcaric Cambisols and is positively correlated with the first PCA component, which means high pH, < 2 μm fraction, Fetot and SOC contents. The second group corresponds to soils developed on loess (Luvisols and Albeluvisols) and is negatively correlated to the first PCA component, which means high 20-50 μm fraction and Sitot contents, low pH, Fetot and SOC contents. The third group corresponds to hypereutric Cambisols and is intermediary between the other two groups, showing more similarities with the soils developed on loess than with the calcaric Cambisols.

Although the PCA did not highlight the impact of land use, pH and Si content in SRO were higher and SOC values lower under cultivation than under forest (Fig. 2). The CEC of the < 2 μm fraction did not differ with land use in Cambisols developed on limestone but was higher under cultivation in Luvisols and Albeluvisols developed on loess (Fig. 2).

Phytoliths of grasses were dominant in the two land uses

The relative abundances of the main phytolith types grouped according to plant types and confidence intervals are given in Table 1. On average, the biogenic silica extracts contained 89 % of phytoliths (range 51-99 %), the remaining fraction consisting of 3 % of other biogenic silica particles (range 0-8 %) and of 9 % charcoal (range 0-48 %); all particles had roughly the same size. For each soil, the mass of phytoliths corrected by phytolith relative abundance (in %) (Table 1) indicates that phytolith pools in soils ranged from < 0.01 to 0.616 g 100 g -1 corresponding to 90 and up to 3978 kg Si ha -1 . For comparison, 800 kg Si ha -1 was found in a rice field from Camargue, France [START_REF] Desplanques | Silicon transfers in a rice field in Camargue (France)[END_REF], and values ranging from 450 to 1400 kg Si ha -1 in French forest soils [START_REF] Bartoli | The biogeochemical cycle of silicon in two temperate forest ecosystems[END_REF], using a similar procedure (heavy liquid extraction).

The phytolith contents were very low in calcaric Cambisols and Albeluvisols (Fig. 2i), for which there was no clear differences due to land use. For Luvisols, we observed clearly higher phytolith contents under forest (Fig. 2i). For the hypereutric Cambisols, a clear difference in the amount of phytoliths between land uses was also observed, although not always in the same direction: for the BS2 site more phytoliths were recorded under forest, while for the 1166 site more phytoliths were observed under cultivation (Fig. 2i).

Grass phytoliths (silica short cells and various silicified epidermal tissues) represented >64 % in all forest and cultivated soil samples, except in forest samples 755-F1, 755-F2, and 1166-F1 (25 % only) (Fig. 3a). Phytoliths assigned to forest species (see Tables in SI part 2 for details) only represented 1-20 % of the total phytoliths sum, except in the forest soil sample 1166-F1 (45 %) (Fig. 3b).

The proportions of ELONGATE DENDRITIC phytoliths, mainly produced by grass inflorescences, tended to be higher in the cultivated soils of the hypereutric Cambisols (Fig. 3c). Again, in the other soils their abundance was too low to discuss any differences between land uses (Fig. 3c). These results showed that in the forest soils, except for the forest site 1166, the forest-indicator phytoliths exhibited a very low abundance relative to grass phytoliths (Fig. 3), suggesting that wild grasses were the main producers of phytoliths in forest understorey. The proportion of corroded phytoliths (ELONGATE and BLOCKY) was not different between forest and cultivated soils, except for samples from sites 1166 and BS2 where it was higher in the forest soil and in the cultivated soil respectively (Table 1).

Mineralogy of the < 2 μm fraction in the topsoils of the seven paired sites

The mineralogical composition of the < 2 μm fractions was rather similar for all forest soils (Fig. 4). All the fractions contained quartz (peak at 3.34 Å in all treatments), kaolinite (peak at 7.13 Å in both AD and EG treatments), illite (peak at 10 Å in all treatments), chlorite (peak at 14 Å in all treatments), poorly crystalline smectites and/or interlayered illite-smectite minerals (large bulgearound 12-11 Å for Na-exchanged samples shifting around 14 Å with glycol with no well-marked peak, more or less marked depending on the considered sample). Soils on loess clearly contained more quartz in their < 2 μm fraction than soils on limestone as shown by their higher quartz:boehmite peak ratio (Table SI3). Cultivation appeared to slightly increase the amount of smectite and interlayered minerals and the amount of illite and decrease the amount of chlorite in soils on loess and in hypereutric Cambisols (Fig. 4 b to d, and Table SI3). In calcaric Cambisols, however, cultivation seemed to result in a decrease of all clay minerals, possibly due to the nondissolution of the carbonate fraction that acts as clay cement, with no clear change in proportion between the different types of clay minerals (Fig. 4a). Transmission electron microscope (TEM) images showed the presence of discrete spherical particles with external diameter of 6 nm (Fig. 5). These spherical particles appear in the form of a disc with a white central part interpreted as being the projection of a hollow sphere.

Si pools in the topsoils of the seven paired sites assessed by chemical extractions

Whatever the soil type, more than 99 % of Sitot was contained in crystalline silicates (Sitot -Sid) (Fig. 2j). The other Si pools (Si in SRO minerals (Sio -SiCaCl2), Si associated to crystalline oxides (Sid -Sio) and Si in phytoliths) represented<0.3 % of the total soil Si. Si in SRO minerals was the most abundant of these pools in calcaric Cambisols, especially when cropped (Fig. 2g), while Si in phytoliths was the most abundant in Luvisols (Fig. 2i). Si associated to crystalline Fe oxides (Sid -Sio) represented also a non-negligible part of the nonsilicate Si pools, notably under forest (Fig. 2h).

Part of this pool may originate from amorphous silica partly dissolved by the DCB reagent [START_REF] Georgiadis | Development of a method for sequential Si extraction from soils[END_REF]. However, the pool associated to crystalline oxides (Sid-Sio) was clearly lower than the Si content in phytoliths, indicating that the phytoliths were not entirely dissolved by DCB.

The SiCaCl2 pool represented a negligible proportion of the Sitot (Table 1) and was larger in soils developed on limestone than in soils developed on loess (Fig. 2f). Cultivation increased both the SiCaCl2 pool, especially in soils on loess, and the SRO-Si pool, especially in soils on limestone parent material (Fig. 2g). Cultivation, however, did not seem to modify the amount of Si associated to wellcrystalline iron oxides and included in crystalline silicates (Fig. 2h andj).

Discussion

Allophanes are responsible for the SiCaCl2 concentration in forest soils

When considering the 21 soil samples of the seven paired sites, the correlations between SiCaCl2 and other soil parameters were in a good agreement with previous studies for the following parameters:

1) positive correlations between SiCaCl2 and the < 2 μ m fraction and its composition (CEC<2μm, Altot, Fetot), pH, CEC (Table 2), as already observed by [START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF], but also between SiCaCl2 and Sio and Sid (that contains Sio see section 2.2.2). Such positive correlations of SiCaCl2 with Sio and Sid were already observed for many soil types by [START_REF] Yanai | Evaluation of available silicon content and its determining factors of agricultural soils in Japan[END_REF].

2) a negative correlation of SiCaCl2 with Sitot (Table 2), which can be explained by the predominance of low solubility minerals such as quartz in the total Si pool [START_REF] Yanai | Evaluation of available silicon content and its determining factors of agricultural soils in Japan[END_REF][START_REF] Landr´e | Do climate and land use affect the pool of total silicon concentration? A digital soil mapping approach of French topsoils[END_REF]. However, contrary to previous studies made for very different soils, SiCaCl2 was not positively correlated with iron oxide (estimated by Fed) contents (a sandy soil with very low content of SiCaCl2; de [START_REF] De Tombeur | Plants sustain the terrestrial silicon cycle during ecosystem retrogression[END_REF] or the amount of phytoliths (estimated by Na2CO3 for Podzol and Chernozem, Li et al., 2020a) (Table 2). Li et al. (2020b) demonstrated that the extraction by Na2CO3 gave very different results from the extraction of phytoliths through heavy liquid separation. In the studied soils, the phytolith abundances estimated through heavy liquid separation were negatively correlated to the pH of forest soils (r--0.826, p-value = 0.022). This result agrees with an increasing solubility of phytoliths with increasing pH [START_REF] Fraysse | Surface chemistry and reactivity of plant phytoliths in aqueous solutions[END_REF]. As a result, the amount of phytoliths in the calcaric Cambisol and in the hypereutric Cambisol 1166, that had a pH value close to 7 (Table SI1), was low. It was also the case for the Albeluvisol despite a soil pH value of 5.2 (Table SI1). In this soil, reducing conditions occurring in winter and resulting in a transient high pH [START_REF] Cornu | Soluble and colloidal translocation of Al, Fe, Si and Mn in an artificially drained French Retisol[END_REF] may be responsible for an increase in the phytolith dissolution over that period. For the Luvisols and the BS2 hypereutric Cambisol, the amount of phytoliths in forest increased as the pH decreased. The absence of correlation between SiCaCl2 and phytoliths suggests that this pool is however not the main pool controlling the SiCaCl2 in temperate soils and that another pool of pedogenic origin must exist as already evoked by [START_REF] Vandevenne | Silicon pools in human impacted soils of temperate zones[END_REF] and [START_REF] Keller | Are clay minerals a significant source of Si for crops? A comparison of amorphous silica and the roles of the minerals type and pH[END_REF].

Sio was very strongly correlated with Alo (Fig. 6a), excluding the cultivated calcaric Cambisol 701 that was very rich in SOC compared to the other cultivated soils (Fig. 2c and discussion below). Contrarily Sio was not significantly correlated with Feo in cultivated plots (r = 0.6, p-value = 0.2) or in forests (r = 0.5, p-value = 0.3). In addition, the Alo/Sio ratio observed for both forest and cultivated soils was close to 2, a typical value for allophanes [START_REF] Harsh | Encyclopedia of Soils in the Environment[END_REF][START_REF] Wada | Allophane and imogolite[END_REF]. Hence, it is likely that the Sio pool consisted in SRO aluminosilicates. SRO aluminosilicates are found in most soils, including those developed on loess and sedimentary rocks as in this study, insofar as their DSi are between 0.1 and 4 mmol L-1 [START_REF] Harsh | Encyclopedia of Soils in the Environment[END_REF]. DSi values close to 0.1 mmol L-1 were indeed obtained for the soil solution of the BS3 soil [START_REF] Cornu | Soluble and colloidal translocation of Al, Fe, Si and Mn in an artificially drained French Retisol[END_REF]. In addition, the spherical particles observed by TEM (Fig. 5) are commonly interpreted as allophanes [START_REF] Kitagawa | The "unit particle" of allophane[END_REF][START_REF] Henmi | Morphology and composition of allophane[END_REF][START_REF] Parfitt | Allophane in new Zeland -a review[END_REF].

Our dataset shows also a strong correlation between the < 2 μm fraction and the allophane contents (Fig. 6b), suggesting that allophanes may either be a product of clay mineral weathering or be part of the clay mineral formation, notably of the finest clays, the smectites and/or interlayered clay minerals that were shown to be present in all the soils considered (Fig. 4). Li et al. (2020a) also observed the formation of allophanes in a pot experiment as the result of weathering. When the Alo/ Sio ratio is above 2, which is the case in most forest soils (Table 1 andSI2), the Alo content in excess is present as Al-OC complexes [START_REF] Harsh | Encyclopedia of Soils in the Environment[END_REF]. These Al-OC complexes partially preclude the formation of short range ordered aluminosilicates [START_REF] Harsh | Encyclopedia of Soils in the Environment[END_REF]. SiCaCl2 was well correlated to allophane content (Fig. 7a), excluding the BS2 site for which part of the Sio extracted was not contained in allophanes contrarily to all the other sites (as calculated by eq. ( 2)).

Therefore, allophanes are assumed to control a fraction of SiCaCl2 in soils. This agrees well with the high SiCaCl2 content observed in Andosols, which contain a large amount of imogolites and allophanes [START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF]. [START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF] concluded that SiCaCl2 was controlled by both sorption and dissolution processes, which are strongly controlled by pH for soils with < 2 μm fractions both larger than 50 g kg 1 and lower than 325 g kg 1. For Indian soils this relation is not linear but polynomial, with a stabilisation of the SiCaCl2 concentrations for pH higher than 7 [START_REF] Meunier | pH as a proxy for estimating plant-available Si? A case study in rice fields in Karnataka (South India)[END_REF]. We used a comparable relationship for the soils considered here excluding the BS2 and 62 forest soils (Fig. 7b). The pH value of the 62 forest plot lower than 4.5 suggests the beginning of podzolisation [START_REF] Caner | Shorttime clay-mineral evolution in a soil chronosequence in Ol´eron Island (France)[END_REF][START_REF] Sauer | Podzol development with time in sandy beach deposits in southern Norway[END_REF][START_REF] Fekiacova | Can Fe isotope fractionations trace the pedogenetic mechanisms involved in podzolization?[END_REF]. The conifer vegetation observed on this plot (Table SI1) is compatible with the onset of such a process. This low pH also corresponds to the upper limit of the Al/Fe oxide buffer range [START_REF] Stützer | Early stages of podzolisation in young aeolian sediments, western Jutland[END_REF].

Podzolisation induces dissolution of the clay minerals, which may be responsible for the higher SiCaCl2 concentration observed with respect to the regression model presented in Fig. 7b. The above findings have the following consequences on the terrestrial Si cycle described by [START_REF] Struyf | Historical land use change has lowered terrestrial silica mobilization[END_REF]. In all soil types, 99 % of Si was contained in crystalline silicates. However, in soils developed on loess almost 60 % of the crystalline silicates consisted of quartz (e.g., Montagne et al., 2008 for BS3), a poorly weatherable mineral. In the calcaric Cambisols, the total Si content in crystalline silicate was much lower than in soils developed on loess and was mainly contained in clay minerals and some quartz and feldspar. Despite this mineralogical difference, we concluded that the main weatherable mineral compartment contributing to SiCaCl2 consists in the clay minerals (Fig. 8). This compartment was almost twofold larger in the calcaric Cambisols than in the soils on loess (Fig. 2a). The clay compartment controlled allophanes (Fig. 6b) which in turns controlled the SiCaCl2 as discussed above.

Allophanes are an easily available source of DSi as shown by dissolution rates much higher than those of phytoliths at pH lower than 9 [START_REF] Ralston | Dissolution rates of allophane with variable Fe contents: implications for aqueous alteration and the preservation of X-Ray amorphous materials on Mars[END_REF][START_REF] Fraysse | Surface chemistry and reactivity of plant phytoliths in aqueous solutions[END_REF]. Consequently, phytoliths did not appear to be the first source of SiCaCl in the considered soils, in agreement with [START_REF] Keller | Are clay minerals a significant source of Si for crops? A comparison of amorphous silica and the roles of the minerals type and pH[END_REF], who found that the clay minerals can compete with phytoliths as a main source of plant Si when the proportion of clay is high.

Reflecting the pH-dependent dissolution rate of phytoliths [START_REF] Fraysse | Surface chemistry and reactivity of plant phytoliths in aqueous solutions[END_REF], the phytolith pool was lower in calcaric Cambisols (pH around 7, Table SI1) than in soils developed on loess (pH around 5, Table SI1). Consequently, the higher SiCaCl2 in calcaric Cambisols compared to soils on loess could be explained by the higher amount of allophanes in calcaric Cambisol, but also, to a lesser extent, by a higher dissolution of the phytolith pool.

Impact of land use on the Si biogeochemical cycle and its various mineralogical pools

As expected, cultivation increased the soil pH (Fig. 2b) and decreased the SOC content (Fig. 2c) in all soils. For soils developed on loess, the pH was two units higher under cultivation than under forest because of liming. For calcaric Cambisols, the pH was one unit higher under cultivation than under forest because of the acidifying capacity of the forest (Fig. 2b, Table SI1). We confirm here that cultivation increases the amount of SiCaCl2 for soils developed on sediments [START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF][START_REF] Vandevenne | Silicon pools in human impacted soils of temperate zones[END_REF][START_REF] Klotzbücher | Impact of agriculture practices on plant-available silicon[END_REF].

The change in pH and SOC induced by cultivation may act on the allophanes and phytoliths pools and therefore modify the terrestrial Si cycle. Based on the literature, a lower phytolith content was expected in cultivated soils compared to topsoils under forest [START_REF] Guntzer | Long-term removal of wheat straw decreases soil amorphous silica at Broadbalk, Rothamsted[END_REF][START_REF] Vandevenne | Silicon pools in human impacted soils of temperate zones[END_REF]. This was clearly observed for Luvisols, and to a lesser extent for the BS2 hypereutric Cambisol (Fig. 2i). This decrease in phytoliths may be explained by both a larger dissolution of phytoliths under cultivation due to the pH increase and phytolith exportation due to harvesting. The relative abundance of phytoliths (Fig. 3) suggests that wild grasses were the main producers of phytoliths in forest undergrowth. This is in agreement with temperate trees being poor Si accumulators [START_REF] Bremond | Advantages and disadvantages of phytolith analysis for the reconstruction of Mediterranean vegetation: an assessment based on modern phytolith, pollen and botanical data (Luberon, France)[END_REF].

Cultivation had no clear impact on the phytolith pool in the calcaric Cambisol or in the Albeluvisols, where the amounts of extracted phytoliths were very low. The proportion of corroded phytoliths (ELONGATE and BLOCKY) was higher in the forest 1166 hypereutric Cambisol (Table 1), that also exhibits a lower amount of phytoliths than under cultivation. Differences in surface areas between grass and forest phytoliths might explain the observed proportions with a higher rate of dissolution of forest phytoliths [START_REF] Wilding | Contributions of Forest Opal and Associated Crystalline Phases to Fine Silt and Clay Fractions of Soils[END_REF].

However, this hypothesis is not supported by the laboratory measurements [START_REF] Fraysse | Surface chemistry and reactivity of plant phytoliths in aqueous solutions[END_REF]. Crop management through the restitution strategy of inedible plant parts and the crop rotation system may also lead to site fluctuation in exportation, maybe explaining the exception of the 1166 site. Indeed, crop straw recycling in agricultural soils have proven to be an efficient strategy for limiting desilication [START_REF] Puppe | Crop straw recycling prevents anthropogenic desilication of agricultural soil-plant systems in the temperate zone-Results from a long-term field experiment in NE Germany[END_REF].

By decreasing SOC (Fig. 2c), cultivation promotes the formation of allophanes that are twice as abundant under cultivation as under forest in the calcaric Cambisols and soils developed on loess (Fig. 2c and5b). We observed Alo/Sio ratios higher than 2 in forest but close to 2 in cultivated soils (Table 1), evidencing a decrease in the Al-OC complexes that prevent allophane formation under forest as discussed above. For the cultivated soils with a pH close to 8, pH may be less favourable to Si adsorption on allophanes that are at zero point of charges [START_REF] Harsh | Encyclopedia of Soils in the Environment[END_REF]. The high SiCaCl2 values may then be explained by Si adsorbed on carbonate surfaces [START_REF] Caubet | Agriculture increases the bioavailability of silicon, a beneficial element for crop, in temperate soils[END_REF]. [START_REF] Struyf | Historical land use change has lowered terrestrial silica mobilization[END_REF] observe that cultivation decreases the exportation of DSi at the watershed outlet. They interpret this decrease in exportation by a reduction in the silicate mineral weathering. On the contrary, our results suggest that this decrease might be due to a higher retention of Si in the soils by the increase of allophane formation under cultivation (Fig. 8).

Conclusion

In temperate soils, phytoliths may not be the main Si pool contributing to bioavailable Si in contrast to tropical soils that are strongly depleted in Si and have kaolinite as dominant clay mineral. Temperate soils indeed exhibit a more complex clay mineral assemblage (smectite, illite, kaolinite and interlayered minerals). Allophanes -either precursors or dissolution products of these clay minerals -appear as the main Si pool contributing to DSi, estimated by CaCl2 extraction. This explains the correlation between SiCaCl2 and the clay fraction reported in the literature. Cultivation seems to impact soil Si pools and DSi in two ways. First, cultivation increases pH by liming in the case of non-calcaric soils and by suppressing forest acidification. The increase in pH enhances Si sorption and phytolith dissolution. Secondly, cultivation decreases soil organic carbon, promoting the formation of allophanes. At pH lower than 4.5, other processes occur, such as podzolisation, modifying the relationship between DSi and pH. So far, there has been a lot of emphasis put on phytoliths as a source of DSi. Our work highlights the importance of allophanes as the missing link between the soil solution and the clay minerals that could explain why clay minerals may be so important, despite lower solubility rates than phytoliths.

Acknowledgments

This work was performed in the framework of the French ANR BioSiSol (ANR-14-CE01-0002) project. The authors are grateful to the CEES team (10.15454/VCZ6BS) (INFOSOL, INRAE, Orl´eans) for their support in site selection, sampling and bulk density measurements and to B. Angeletti for ICP measurements at CEREGE. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.geoderma.2022.115722.

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.