Thomas Loiseau

Dominique Arrouays
email: dominique.arrouays@inrae.fr

Anne Richer-De-Forges

Philippe Lagacherie²

Christophe Ducommun

Budiman Minasny

Density of soil observations in digital soil mapping: A study in the Mayenne region, France

Keywords: Digital Soil Mapping, topsoil particle-size distribution, sampling strategy, sampling density, prediction performance, Multiple soil classes, France

The density of soil observations is a major determinant of digital soil mapping (DSM) prediction accuracy. In this study, we investigated the effect of soil sampling density on the performance of DSM to predict topsoil particle-size distribution in the Mayenne region of France. We tested two prediction algorithms, namely ordinary kriging (OK) and quantile random forest (QRF). The study area is a region of ~5000 km 2 with the highest density of field soil observations in France (1 profile per 0.64 km 2 ). The number of training sites was progressively reduced (from n = 7500 to n = 400, corresponding to 1 profile per 0.7 km 2 to 1 profile per 13 km 2 ) to simulate the different density of observations. For OK and QRF, we tested random subsampling for splitting the data into training and testing datasets using k-fold cross validation. For QRF we also tested conditioned Latin hypercube sampling based on the point coordinates or the covariates. The results indicated that, with increasing density of observations, OK performed as well or even better than QRF, depending on the particle-size fraction. For silt prediction, OK was systematically better than QRF. However, the prediction intervals were much larger for OK than for QRF, and OK did not seem to estimate uncertainty correctly. Overall, the performance indicators increased with the density of observations with a threshold at about 1 profile per 2 km 2 which suggests that the main limitation of DSM prediction accuracy using QRF is the amount of data collected in the field, not the type of calibration sampling strategy. Future DSM activities should focus on gathering more field observations.

Introduction

It is generally accepted that a major limitation to Digital Soil Mapping (DSM) prediction performance is the density of soil observations that carry up-to-date information (e.g., Arrouays et al., 2014a[START_REF] Arrouays | Soil legacy data rescue via GlobalSoilMap and other international and national initiatives[END_REF][START_REF] Samuel-Rosa | Open legacy soil survey data in Brazil: geospatial data quality and how to improve it[END_REF]. [START_REF] Samuel-Rosa | Do more detailed environmental covariates deliver more accurate soil maps?[END_REF] evaluated whether investing in more spatially detailed environmental covariates improves the accuracy of digital soil maps. Their conclusions showed that more detailed covariates only result in a modest increase in prediction performance and it may be more useful to spend extra resources on collecting more soil observations. Indeed, substantial discrepancies in soil observations' density exist among, and even within, national soil databases [START_REF] Morvan | Soil monitoring in Europe: a review of existing systems and requirements for harmonisation[END_REF][START_REF] Arrouays | Soil legacy data rescue via GlobalSoilMap and other international and national initiatives[END_REF]. [START_REF] Somarathna | More data or a better model? Figuring out what matters most for the spatial prediction of soil carbon[END_REF] tested various machine learning algorithms to predict soil carbon and found that the density of observations was more important than the type of machine learning models. While many DSM studies focused on comparing various machine learning models [START_REF] Padarian | Machine learning and soil sciences: A review aided by machine learning tools[END_REF], less work focused on the impact of the density of observations. By using a large dataset of pseudo values of clay content obtained from hyperspectral data, [START_REF] Lagacherie | 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[END_REF] confirmed the importance of sampling density on DSM performances but showed that this importance diminished as sampling densities increased and that other spatial characteristics of the soil sampling (completeness and evenness in space) also had a strong effect.

Even within harmonized soil conventional mapping programs, the methodology and spatial coverage are far from homogeneous, even at the national scale. This is the case in France, where a national program for mapping soil-scapes at 1:250,000 is finalised [START_REF] Laroche | Le programme Inventaire Gestion et conservation des Sols. Volet Référentiel Régional[END_REF][START_REF] Richer-De-Forges | Mapping of Soils and Land-Related Environmental attributes in France: analysis of end-users' needs[END_REF]. In their prospective analysis on the future of the French national soil-mapping program, [START_REF] Voltz | Possible futures of soil-mapping in France[END_REF] proposed to move from conventional soil mapping to DSM and to improve the density of observed soil profiles to homogenize and enhance the accuracy of 1:250,000 soil maps for France. In practice, the French program of 1:250,000 soil mapping has been done using conventional mapping by administrative units i.e., the French "Départements" the mean area of which being about 6,000 km². Within these "Départements", the density of soil observations is highly variable, from 1 observation per 0.64 km² to 1 observation per 123 km². Three thresholds were set to indicate the quality of the maps (InfoSol, 2005): (i) at least 1 profile per 60 km² as a minimum standard designated as "operational level", (ii) more than 1 profile per 40 km² for a higher quality "advanced level" and (iii) more than 1 profile per 20 km² as "optimal level". These thresholds were set by experts and used to deliver labels of quality by the Ministry in charge of Agriculture, but their relevance was never assessed quantitatively. In particular, we never tested if increasing the density of points substantially would increase the performance of the predictions when using DSM techniques to map a "Département".

In this study, we chose the French "Département" with the highest soil profile density observations. We tested how the density of training profiles affects the prediction of topsoil particle size distribution, namely clay, silt, and sand contents. This paper has two objectives: i) to assess if we could identify thresholds of the density of soil profiles for prediction performances when applying DSM techniques to predict some basic soil properties and, ii) to compare the performances of a machine learning algorithm (quantile random forest) to pure spatial interpolation, ordinary kriging, as a function of the density of observations.

Material and methods

2.1.Study area and general framework of the study

We chose the "Mayenne" French "Département". (~8,100 points for a total area of 5208 km², corresponding to a density of one profile per 0.64 km²), which has the highest density of observations for a French department in the framework of the systematic conventional survey of soils at the scale of 1:250,000. Figure 1 shows the study area with the original sampling points from the soil surveyor and the elevation. Except for some urban areas, sampling points are rather homogeneously distributed. Of course, the density of points is not regular, because conventional survey allows the surveyor to place the points for maximum information on the soil-landscape, but we can reasonably expect that most of the soil/landscape situations are represented in this sampling. The region is located at the boundary of sedimentary deposits (mainly calcareous and loess) and crystalline and metamorphic rocks. Thus, we expect that we would be able to see contrasting spatial structures of soil texture within this region. In addition, the area is covered by interesting spatial covariates, namely, (i) a harmonized and revised lithological map (1:50,000 scale; [START_REF] Vernhet | Carte géologique harmonisée du département de la Mayenne[END_REF], recoded by Bialkowski et al., 2019), and (ii) an airborne gamma-ray survey.

Soil sampling and analysis

Soil sampling was done in the framework of a conventional soil survey; its density was adapted to the soil variations observed in the field by the soil surveyors. At each field sampling location, the soils were sampled either by digging a soil pit or by augering. The topsoil organo-mineral horizon (excluding O horizons when present) was sampled. For most of the agricultural plots, the corresponding thickness was the ploughed layer (from 0-20 to 0-30 cm, with an average thickness of 24 cm) whereas the thickness ranged from 0-5 to 0-30 cm (with an average of 20 cm) and from 0-5 to 0-25 cm (with an average of 18 cm) under permanent pastures and forest respectively. Samples were air-dried at 30°C to a constant mass and then gently crushed to a 2 mm sieve. The particle-size analysis was performed using the Robinson pipette method [START_REF] Robinson | The dispersion of soils in mechanical analysis[END_REF].

As the fractions clay (< 2 µm), silt (2-50 µm), and sand (50-2000 µm) sum to 1000 g kg -1 , the data were transformed with the additive log-ratio (alr-transform, [START_REF] Aitchison | The statistical analysis of compositional data[END_REF] for the spatial prediction of compositional data [START_REF] Lark | Cokriging particle size fractions of the soil[END_REF]. The alr-transformed variables were calculated as follows (Eqs (1 and 2)):

ln 1 ln 2
The alr predicted variables were then back-transformed to sand, silt, and clay through the alr inverse transformation (Eqs (3-5)):

. ! "#$% &'( 1 ) ! "#$% &'( ) ! *+#, &'( -1000 3 . ! *+#, &'( 1 ) ! "#$% &'( ) ! *+#, &'( -1000 4 12 . 1 1 ) ! "#$% &'( ) ! *+#, &'(
-1000 5

To assess the confidence interval of our predictions, we calculated the variance of the backtransformed clay, sand, and silt fractions. Román Dobarco et al. (2019) showed that since the prediction and variance are independent, the Taylor analysis method is more efficient than using Monte Carlo simulation to estimate the prediction variance of each particle-size.

Utilising the first-order Taylor analysis, the variance of prediction can be derived as follow (Eqs (6-8)):

• 4 5 . 6 7 "#$% &'( 8 9:::;< ='&> &'( ;9?< ='&> &'( ?< @A'B &'( C < D='&> &'( C ;9?< ='&> &'( ?< @A'B &'( C² 8 ) 7 *+#, &'( 8 F 9:::< ='&> &'( < @A'B &'( ;9?< ='&> &'( ?< @A'B &'( C D G 8 ) 2HIJ $#K , $#K 9: M < D='&> &'( < @A'B &'( ;9?< @A'B &'( C ;9?< ='&> &'( ?< @A'B &'( C D (6)
• 4 5 . 6 7 *+#, &'( 8 9:::;< @A'B &'( ;9?< ='&> &'( ?< @A'B &'( C < D@A'B &'( C ;9?< ='&> &'( ?< @A'B &'( C² 8

)

7 "#$% &'( 8 F 9:::< ='&> &'( < @A'B &'( ;9?< ='&> &'( ?< @A'B &'( C D G 8 ) 2HIJ $#K , $#K 9: M < ='&> &'( < D@A'B &'( ;9?< ='&> &'( C ;9?< ='&> &'( ?< @A'B &'( C D (7) • 4 5 12 . 6 7 "#$% &'( 8 F 9:::< ='&> &'( ;9?< ='&> &'( ?< @A'B &'( C D G 8 ) 7 *+#, &'( 8 F 9:::< @A'B &'( ;9?< ='&> &'( ?< @A'B &'( C D G 8 ) 2HIJ $#K , $#K 9: M < ='&> &'( < @A'B &'( ;9?< ='&> &'( ?< @A'B &'( C D (8)
For more detailed explanations on the use of first-order Taylor analysis in spatial models, we refer to [START_REF] Heuvelink | Propagation of errors in spatial modelling with GIS[END_REF].

2.3.Covariates

The list of covariates used in this study is provided in Table 1. All covariates were transformed to a 90 m resolution grid according to the GlobalSoilMap specifications (Arrouays et al., 2014b), either through downscaling for coarser-resolution layers, or upscaling for the finer ones. The upscaling was done in two ways: 1) by the average aggregation for quantitative data, 2) by a majority vote estimation for the qualitative covariates. In addition, we transformed the classes from the qualitative covariates to binary information (0: absence and 1: presence) to assess each class's importance in our model.

2.4.Digital soil mapping modelling

In order to assess the effect of point density on the performance of the predictions, we reduced the number of training samples progressively and adopted three different strategies to predict particle-size distribution over the "Département". To keep comparable testing sizes for large and small datasets, we selected 8100 to 600 sample points with a 1000 to 200 increment (Table 2). We used 200 testing points for cross-validation when the number of training points was larger than 2000, and 500 testing points when the number of training points was smaller or equal to 2000. This strategy aimed at keeping enough testing points for each fold, and as many folds as possible. Each splitting procedure was repeated 100 times. We tested four different prediction methods, summarized as follows:

1) Ordinary Kriging (OK), to assess the prediction of soil texture, only using soil sample location and measured particle-size values as information. Samples were randomly selected to form various sizes of training data.

2) Quantile Regression Forest (QRF) using the R package quantregForest [START_REF] Meinshausen | Quantile Regression Forests[END_REF] with covariates listed in Table 1. We chose this technique as it has been shown to be robust, and it enabled us to derive quantile distributions instead of a single mean value [START_REF] Vaysse | Using quantile regression forest to estimate uncertainty of digital soil mapping products[END_REF]. Samples were randomly selected as in OK.

3) QRF as previously, but we chose to sample the points using the conditional Latin Hypercube Sampling (cLHS, [START_REF] Minasny | A conditioned Latin hypercube method for sampling in the presence of ancillary information[END_REF]) method based on their coordinates to reduce the clustering of samples.

4) QRF as previously, but we chose to sample the points using cLHS applied on the covariates to ensure even coverage of the covariates.

2.5.Assessing the performance of the predictions

The performance of the predictions was evaluated through a k-fold cross validation by keeping a minimum of 200 to 500 samples apart for testing, depending on the strategy of partitioning the data. We assessed the model performance through a 100 times replication of the cross-validation method by reporting the medianvalue of the following indicators:

The coefficient of determination (R²), which measures the adequacy between our model and our observed value.

N 8 1 O ;1 O N PQP $RSTU,<R 8 C O 1 O V O 1 9
where

N PQP $RSTU,<R 8 ∑ Y A * Y̅ \ A]^² ∑ \ A]^YA Y̅ ² (10)
Here, _ + and _̅ are the value of the observation for point i and the mean of all observations, respectively; _ + * is the value of the prediction for point i.

The root-mean-square error (RMSE), which provides information on the statistical dispersion of our predictions in relation to our observations.

N` a b 1 c _ + O _ + * 8 P +d9 11
The concordance coefficient (CC; Lin, 1989), which is given as

2e7 Y * 7 Y 7 Y 8 ) 7 Y * ² ) _̅ O _̅ * 8 -100 12
where 7 Y 8 and 7 Y * ² are the observation and prediction variances, respectively, e is their correlation coefficient, and _̅ * is the mean of the prediction

The bias of the predictions, which is the difference between the mean predicted and mean observed values.

`f 1 f g 1 c P +d9 _ + * O _ + 13
The PICP (prediction interval coverage probability). This indicator assesses the uncertainty between the observed and predicted distribution. The PICP is the probability that the target of an input pattern lies within the prediction limits [START_REF] Shrestha | Machine learning approaches for estimation of prediction interval for the model output[END_REF]. The PICP was calculated from the confidence interval of 90% with the extreme quantile 5% and 95% of the model prediction.

hi h ∑ _ j% * l _ + l _ mj% * P +d9 -100 14

Results and discussion

3.1.Summary statistics

The distribution of the particle-size fractions in g.kg -1 , over the Mayenne department, is shown in Table 3. We observed a large dominance of the silt fraction over the samples.

Moreover, testing the skewness and kurtosis indices, silt and sand seemed to follow a rather symmetric distribution with fewer and less extreme outliers than the normal distribution. The clay distribution was positively skewed and had a high kurtosis value, which shows its asymmetry and more outliers in the right part of its distribution.

3.2.Performance of the predictions

In order to assess the performances of each sampling strategy, we compared the performance of the predictions for each sampling density and algorithms used. The R 2 for the QRF algorithms for silt and sand showed a sharp increase with an increasing number of training points until about 1800 and 2500 respectively, followed by a small increase (Fig 2). This shape is similar to those obtained by [START_REF] Lagacherie | 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[END_REF]. These threshold were less evident for clay. The three QRF algorithms gave quite similar results.

Interestingly the R 2 for silt was always higher when using OK than using QRF methods. This may be due to the fact that the covariates we used did not capture silt distribution effectively.

Moreover, when the number of training samples was very large, OK R² gave similar results compared to QRF for clay and slightly better results for sand. Note, however, that OK gave poorer results than QRF for clay and especially for sand at the lowest densities of sampling and that for these fractions, the decrease of R 2 was more pronounced for low sampling densities. until reaching a threshold of about 2500, followed by a slight decrease thereafter. No real trend was observed for silt and sand. No substantial difference was observed between the algorithms (Fig. 3) except for silt, in which OK always performed slightly better than the QRF methods, especially at high sampling densities, which is consistent with the results of R². Note that the RMSE values are always lower than the interquartile range (see Table 3), suggesting that the mapping is effective. The red line separates the strategies using 500 testing points from those using 200 testing points.

Q3

The mean error for clay seemed to be constant for all models at a very low value for clay (around -2 g.kg -1 ). Silt and sand exhibited much larger ME, even if the OK presented better results for these fractions. For all QRF models, the ME median value stayed somewhat similar. Interestingly, silt and sand showed opposite trends, increasing ME for silt until about The CCC indicator increased with increasing density for all fractions (Fig. 5). Interestingly, OK worked better than QRF for large sampling densities, and even for all densities for silt.

For sand, the CCC trend was similar to R 2 , with a sharp increase with the increasing density of training points and a lower increase when the number of training points became higher (> 2500). For the highest training densities, it was higher for OK than for QRF, especially for silt. On the contrary, for the lowest densities, it became lower for OK than for QRF. The red line separates the strategies using 500 testing points from those using 200 testing points.

For clay (Fig. 6), the PICP was always lower than 90%, which shows that the 90% prediction intervals were slightly under-estimated. This was more evident for OK. From a practical point of view, this shows that kriging smoothed out high/low values, thus under-estimated the prediction intervals. The PICP is mostly over 90% for the prediction of silt by OK, although it

shows acceptable values (close to 90%) for the QRF models. This shows that contrary to clay, OK over-estimated PIs for silt prediction, which counter-balanced the better prediction results obtained by OK. For sand prediction, all PICPs were larger than 90%. The contrast between the trends observed for clay and sand could be due to the negative correlation between these two fractions. We observed a rather erratic behavior of the PICP for the lowest calibration sample values. This is likely due to the low number of training and testing samples and the small number of k-fold cross-validations. This is consistent with the findings of Lagacherie et al. (2019), who showed that a small number of validation samples is not robust enough to assess the performance of the model predictions. However, this erratic behavior remains in a narrow range and we should recall here that we replicated 100 times the splitting between training and testing samples. There were no substantial differences between the 3 different sampling strategies for the QRF algorithms. This may be because all subsampling strategies that we tested were denser than the sampling strategies commonly used in France. Moreover, when the number of selected samples is large, cLHS has no advantage over other sampling designs [START_REF] Wadoux | Sampling design optimization for soil mapping with random forest[END_REF].

When the number of training samples was equal or larger than 2500, some QRF performance indicators were slightly improved with increasing sampling density. This suggests that collecting new data or rescuing more legacy data should further improve the prediction accuracy, as indicated by [START_REF] Arrouays | Soil legacy data rescue via GlobalSoilMap and other international and national initiatives[END_REF] and [START_REF] Samuel-Rosa | Do more detailed environmental covariates deliver more accurate soil maps?[END_REF]2020). When the sampling density was very high, OK performed as well, or even better than machine learning methods, which validates the general framework proposed by [START_REF] Minasny | Methodologies for global soil mapping[END_REF] for global digital soil mapping, that was thereafter adopted by the GlobalSoilMap initiative (Arrouays et al., 2014b). Interestingly, the silt fraction had a particular behavior compared to other fractions, i.e., except for PICP, the performance was better when using OK for most sampling densities. This may be due to the fact that the covariates that we selected were less responsive for silt than for the other fractions. A threshold of the number of points (around 2500) was identified, under which most of the performance indicators performed worse. Over this threshold, the performance indicators remained stable or increased more slowly with the density of points, suggesting that the main limitation of the goodness of DSM predictions is the amount of data collected in the field. We note that the number of testing points was changedfrom 200 to 500 at 2500 training points, which may affect the results. To keep a minimum number of k-fold the same time, we needed to change the number of testing points. But we assessed the model performance through a 100 times replication of the cross-validation method. This avoids a random chance of under-or over-estimation of error. Moreover, our results are consistent with those obtained recently by [START_REF] Lagacherie | 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[END_REF].

It is worth noting that this threshold of 2500 points (about 1 profile per 2 km 2 ) corresponds to a nearly ten times larger density that the density recommended in the traditional 1:250,000 soil mapping programs, even for those considered as optimal (1 profile per 20 km 2 ). The conventional survey recommendation at 1:250,000 may be adequate if applied by a skilled soil surveyor, but certainly not for fine resolution DSM applications. The large majority of the French departments have a number of profile observations much less than this threshold of 1 profile per 2 km². One can, therefore, assess the enormous efforts needed if the objective is to produce fine resolution DSM products with an acceptable level of quality at the department level using only training points and available covariates. This is even more critical if we consider that the most important covariates in our case (airborne gamma-ray, [START_REF] Loiseau | Could airborne gamma-spectrometric data replace lithological maps as co-variates for digital soil mapping of topsoil particle-size distribution? A case study in Western France[END_REF] are not available for the entire French territory.

3.3.Maps of predictions

To assess the influence of increasing the number of calibration points on mapping predictions, we mapped several examples for the four trials that we tested for each particle-size fraction (Figures 7 to 9).

The trend showed in these figures is a stabilization of the spatial prediction with the increasing number of training points. This trend is also more pronounced when we focus on OK. The increasing number of calibration points allowed the models to describe more detailed spatial structures, and decreased the smoothing effect of the kriging.

When comparing QRF methodologies sampling, few differences were observed between sampling densities. However, for the lowest training densities, the predictions were more pixelated and some spatial patterns were less visible on the maps. At thresholds ranging from 1800 to 2500 training points, the general patterns seemed to remain stable for clay, silt, and sand.

For clay (figure 7), the spatial patterns with high values (over 200 g.kg -1 ) seemed to remain stable from 1800 training points in the middle part of the department. The predictions of sand (figure 9) presented also similar variations, with a quite similar stabilization of the predictions, for all QRF models, over a threshold close to 1800 points. Several thresholds, ranging from 1800 to 2500 training points, can be observed for the prediction maps of particle-size over the Mayenne department. Indeed, the sampling strategies did not impact the prediction maps compared to the number of training points. However, as shown in the precedent observations, the particle size variability was not entirely explained by our selected covariates. In complex regions, more observations could explain the spatial variation of the soil in more detail. Moreover, using good quality covariates increased a lot in the details of the spatial prediction. In comparison to OK, QRF seemed to indicate the variability of the predictions closer to the observed soil properties variability and produced less smoothed predictions. There are limitations in this study; we still missed some crucial covariates and that our selection of training samples was based on observations collected by a traditional soil survey, which may not be statistically optimal.

3.4.Interval of predictions and mapping

To assess the influence of the training dataset size over our textural predictions, we estimated and mapped their 90% prediction interval (PI). Figures 10 and11 show the PI of each sampling strategy for training sets of 400 and 2500 points. We observed for each sampling, a decrease of the PI with the increase of the number of training points, except for coordinates sampling where the PI seemed relatively stable or better, depending on the particle size, with a smaller dataset (Figure 9). For clay, PIs were lower and showed the most extensive variations over the metamorphic area, whereas the uncertainty was much larger with a few points for each sampling strategy.

Silt showed large variations in the North, but also in the South-Est. This can be due to its variability and the fact that complex deposition processes could not be captured by our model and covariates, even with a large set of training points (Figure 11). The global variation of sand over the department resulted from the cumulative error from the back transformation and presented the largest variation compared to clay and silt. Note that the PIs were always much larger for OK than for QRF models and that, as expected, the PI intervals strongly decreased when going from a small number to a large number of training points. The spatial patterns of PIs were very different from OK to QRF as already observed by [START_REF] Vaysse | Using quantile regression forest to estimate uncertainty of digital soil mapping products[END_REF]. Note also that the PI maps were very similar for all the QRF models, suggesting that the sampling strategy is not important at this high number of sampling points. These results were expected and supported environmental information integration in DSM models for mapping against pure geostatistical interpolation. With a large set of data, the QRF sampling strategies were indifferent, and the variability of prediction was more affected by the number of calibration samples. However, for small size of training points, cLHS appeared to be useful (i.e., for QRF, sampling capturing the maximum combinations of the values of relevant covariates). We suggest that for planning soil sampling for DSM, the use of sampling design that aims to capture the maximum of information is required. There are still discussions about the best sampling strategy for a small number of samples: it is beneficial to use methods such as cLHS or k-means clustering based on covariates when covariates exhibit contrasted values, which is the case in the present study. However, if the contrast of covariates is smoother, then it may be more important to get a rather regular coverage of the geographical space. Overall, increasing the number of training points to more than 2500 led to a relatively small increase in prediction performances, which suggests that acquiring more points may not be worth considering the gain in accuracy that can be achieved. This observation might be different if we consider more complex models such as deep learning [START_REF] Padarian | Using deep learning for digital soil mapping[END_REF][START_REF] Ng | The influence of training sample size on the accuracy of deep learning models for the prediction of soil properties with near-infrared spectroscopy data[END_REF] or new (or finer resolution) covariates that can explain soil texture variability better. This proposal will be tested in future work.

Conclusions

We evaluated the effect of soil observations density on the digital soil mapping model's performance for topsoil particle-size distribution with four sampling strategies and with a decreasing number of training points. We demonstrated a threshold of sampling density, and that sampling with a density less than about 1 profile per 2 km 2 could lead to a substantial decrease in the performance accuracy. This result is significant, considering that the density recommended for conventional soil mapping at 1:250,000 is 1 profile per 20 km 2 , which is 10 times smaller. We showed that sampling strategies based on covariates or coordinates with cLHS did not produce different results from random sampling when the number of training samples was large. We showed that increasing the number of training samples produced better improvements in the predictions maps when using QRF models compared to ordinary kriging. We conclude that the main limitations of DSM prediction accuracy are the amount of data collected in the field and high-quality covariates. Future DSM activities should focus on gathering more field observations. However, increasing the number of training points may in some cases not be worth considering the gain in accuracy that can be achieved.

Fig. 1 .

 1 Fig. 1. Map of the study area with soil sampling points used for the conventional mapping and

  Figures 2 to 6 show the evolution of selected statistical indicators (R², RMSE, ME, CCC, and PICP) as a function of training sample densities.

Fig 2 .

 2 Fig 2. Median R² for testing: comparison for a) Clay, b) Silt and c) Sand fraction for each

Fig 3 .

 3 Fig 3. Median RMSE for testing: comparison for a) Clay, b) Silt and c) Sand fractions for

Fig 4 .

 4 Fig 4. Median error for testing: comparison for a) Clay, b) Silt and c) Sand fraction for each

Fig 5 .

 5 Fig 5. Median Concordance for testing: comparison for a) Clay, b) Silt and c) Sand fractions

Fig 6 .

 6 Fig 6. Median PICP for testing: comparison for a) Clay, b) Silt and c) Sand fraction for each

Fig 7 .Fig 8 .

 78 Fig 7. Clay predictions comparison for the four maps a) OK, b) QRF random, c) QRF 352

Fig 9 .

 9 Fig 9. Sand predictions comparison for the four maps a) OK, b) QRF random, c) QRF

Fig 10 .

 10 Fig 10. Particle-size prediction interval (90%) for the four maps using 600 points (400 for

Fig 11 .

 11 Fig 11. Particle-size prediction interval (90%) for the four maps using 3000 sampling points

Table 1 .

 1 Spatially exhaustive co-variates used in the Mayenne department.

		Covariates		Resolution Data type	Reference
	Climate
	Climate type				250 m	Qualitative Joly et al., (2010)
	Vegetation
	Forest type				1:25 K	Qualitative	Inventaire Forestier National (2006)
	Land use 2016				10 m	Qualitative CESBIO (2016)
	Principal (N=3)	Component		NDVI	500 m	Quantitative Loiseau et al., (2019)
	Soil maps
	Clay content France (%)		500 m	Quantitative Ballabio et al., (2016)
	Silt Content France (%)		500 m	Quantitative	Ballabio et al., (2016)
	Sand content France (%)		500 m	Quantitative Ballabio et al., (2016)
	Rate	of	river	network
	development and persistence	1:50 K	Quantitative BRGM (2014)
	(IDPR)
	Topography
	Elevation				25 m	Quantitative IGN (2014)
	Compound topographic index (CTI)	25 m	Quantitative	IGN (2014)
	Multi-Resolution Valley Bottom Flatness (MRVBF)	25 m	Quantitative IGN (2014)
	Multi-Resolution Flatness (MRRTF)	Ridge	Top	25 m	Quantitative IGN (2014)
	Roughness				25 m	Quantitative IGN (2014)
	Curvature				25 m	Quantitative IGN (2014)
	Exposition				25 m	Quantitative IGN (2014)
	Slope position				25 m	Quantitative IGN (2014)
	Slope, cosines(slope)		25 m	Quantitative IGN (2014)

Table 2 .

 2 Training and testing datasets used in this study.

	Total samples (n) 600 800 1000 1200 1400 1600 1800 2000 3000 4000 5000 6000 7000 8100
	Training samples 400 600 800 1000 1200 1400 1600 1800 2500 3500 4500 5500 6500 7600

Table 3 .

 3 Summary statistics of particle-size distribution of the samples over the Mayenne department.

		Clay	Silt	Sand
	Min	18,83	20	19,06
	Q1	141	448	189
	Mean	177,22	529,26	293,33
	Median	164	540	267

coordinates and d) QRF covariates, when the number of training data increases.

Acknowledgments

The airborne gamma-spectrometric data used in this study were made available by BRGM to INRAE, under the framework of this study and the license agreement n°2019/04. The sampling and most of the soil analyses were funded by a French Scientific Group of Interest on Soils, the "GIS Sol", involving the French Ministry for Ecology and Sustainable Development, the French Ministry of Agriculture, the French Agency for Energy and Environment (ADEME), the National Institute for Agronomic Research (INRAE), the Institute for Research and Development (IRD), the National Forest Inventory (IFN) and the French Agency for Biodiversity, and by local department funds. We thank all of the soil surveyors and technical assistants involved in sampling the sites, preparing soil samples, and producing soil analyses. This research was conducted within the "Centre d'Expertise Scientifique Cartographie Numérique des Sols" granted by the CNES TOSCA program. We are grateful to the INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020. Migale bioinformatics Facility, doi: 10.15454/1.5572390655343293E12) for providing computing and storage resources. We thank Mercedes Román Dobarco for her help with the implementation of the Taylor first-order analysis. Thanks to Nicolas Saby for fruitful discussions about this work. The work of Thomas Loiseau was funded by the French Ministry in charge of agriculture. We thank David G. Rossiter and an anonymous reviewer for their helpful and constructive comments on the paper. D.A. is the coordinator, B.M. is a member, and A.R.d.F. and P.L. are collaborators of the research consortium GLADSOILMAP supported by LE STUDIUM Loire Valley Institute for Advanced Studies through its LE STUDIUM Research Consortium Programme.

Authors' contribution


Mathieu Vangrevelynghe

Lucie Le Nouvel

Cécile Pesenti

Rodolphe Sonnier
email: rodolphe.sonnier@mines-ales.fr

Laurent Ferry

Emmanuelle Gesta

Christian Lagrève

A method to quantitatively assess the modes-of-action of flame-retardants

Keywords:

The assessment of the modes of action of flame retardants (FRs) and their quantification is a challenging issue. This study highlights that the condensed phase efficiency θ can be assessed by combining different flammability tests (pyrolysis-combustion flow calorimeter PCFC, Limiting Oxygen Index LOI, cone calorimeter). Relations were proposed to calculate LOI and pHRR in cone calorimeter from several parameters including θ . A satisfying agreement between experimental values and calculated ones were found for a set of several dozens of FR-free polymers assuming that θ = 1 (i.e. no action in condensed phase).Then, θ was evaluated for EVA/PE blends containing Aluminium TriHydroxyde (ATH) and/or Magnesium DiHydroxyde (MDH) and for ABS containing Ammonium PolyPhosphate (APP) and/or DecaBro-moDiPhenylOxide (DBDPO) assuming that LOI and pHRR of flame-retardant compositions follow the same relationships as those of unfilled polymers. The evolution of θ (as well as χ, i.e. the combustion efficiency) highlights the modes of action of FRs and can contribute to the discussion about possible synergy or antagonism between FRs. Especially, an apparent synergism between APP/DBDPO is observed in cone calorimeter from a performance point of view. Nevertheless, the calculation of χ and θ challenges this conclusion.

Introduction

Polymers are easily flammable, forcing to add flame-retardants to formulations. Many flame retardant systems improve the flame retardancy through different modes-of-action [START_REF] Laoutid | New prospects in flame retardant polymer materials: from fundamentals to nanocomposites[END_REF] in condensed phase (through char promotion or endothermic decomposition), or in gas phase (through radical trapping or dilution of fuels). Synergism between FRs is also often desired to improve the fire performances. For example, Ethylene Vinyl Acetate compounds are commonly used in the wire and cable industry, due to their outstanding properties (flexibility and elongation at break). EVA is flame retarded with high amount of fillers (up to 65%), especially with aluminium and/or magnesium hydroxide (ATH & MDH) [ 2 , 3 ]. In order to improve the fire retardancy of EVA compounds, many fillers were used in combination with ATH or MDH like nanofillers such as montmorillonite [START_REF] Laoutid | Intumescent mineral fire retardant systems in ethylene-vinyl acetate copolymer: effect of silica particles on char cohesion[END_REF][START_REF] Beyer | Flame retardant properties of EVA-nanocomposites and improvements by combination of nanofillers with aluminium trihydrate[END_REF][START_REF] Cárdenas | Mechanical and fire retardant properties of EVA/clay/ATH nanocomposites -Effect of particle size and surface treatment of ATH filler[END_REF][START_REF] Wang | The utility of nanocomposites in fire retardancy[END_REF] , or phosphorous based additives [START_REF] Lv | Flame retardancy and mechanical properties of EVA nanocomposites based on magnesium hydroxide nanoparticles/microcapsulated red phosphorus[END_REF] . Other well-established synergistic systems can be cited, like halogen-antimony synergy [START_REF] Khanna | Flame -Retardant Polymeric Materials[END_REF] which has an effect in both condensed and gas phases thanks to Sb 2 O 3 reaction with hydrogen halides. Ternary intumescent systems are another synergistic systems which promote the formation of a charred layer by a combined reaction between an acid source, a char forming agent and a blowing agent [ 10 , 11 ].

Assessing the relative efficiency of each mode-of-action would be useful to design better FR systems and to identify synergism between FRs. While actions in gas phase can be estimated through the measurement of the combustion efficiency χ, the efficiency of the action in condensed phase, i.e. the relative efficiency of heat and mass transfer at the solid surface (often called θ ) cannot be easily evaluated. Moreover, the importance of modes-of-action depends on the fire scenario (i.e. on the flame test). Therefore, the efficiency of a FR may be significantly different from one test to another because test conditions and pass/fail criteria greatly differ. In LOI, the main criterion governing rating is the flame extinction after removal of the ignition source. In cone calorimeter, the burning is forced by exposing the material to a radiative source during the whole test.

As an example, Cogen et al. [START_REF] Cogen | Correlations between pyrolysis combustion flow calorimetry and conventional flammability tests with halogen-free flame retardant polyolefin compounds[END_REF] compared performances of halogen-free flame retardants in EVA/PE blends for wire and cable industry using different tests. They found some correlations between PCFC and cone calorimeter, but no correlation was no-ticed between LOI and PCFC. Schartel et al. [START_REF] Schartel | Pyrolysis combustion flow calorimeter: a tool to assess flame retarded PC/ABS materials?[END_REF] studied seventeen PC/ABS blends retarded with various additives. They found correlation between total heat release (THR) and char residue with the limiting oxygen index (LOI). But they showed a poor correlation with the heat release capacity (HRC or sumHRC) with a high dispersion of LOI from 15 to 45 for HRC close to 400 J/g.K.

This lack of correlation may be also used to assess the modesof-action of FRs. Sonnier et al. [START_REF] Sonnier | Combining cone calorimeter and PCFC to determine the mode of action of flame-retardant additives[END_REF] have evaluated the barrier effect of different fillers combining two methods (namely PCFC and cone calorimeter), based on the fact that some effects like flame inhibition or barrier effect cannot be observed in PCFC [START_REF] Sonnier | Pyrolysis-Combustion Flow Calorimetry: a Powerful Tool To Evaluate the Flame Retardancy of Polymers[END_REF] but are efficient in cone calorimeter.

Different studies proposed equations to calculate LOI using various parameters. For example, Van Krevelen [START_REF] Van Krevelen | Some basic aspects of flame resistance of polymeric materials[END_REF] proposed a relation between oxygen index and char yield measured in thermogravimetry (anaerobic pyrolysis). Lyon et al. [START_REF] Lyon | Screening flame retardants for plastics using microscale combustion calorimetry[END_REF][START_REF] Lyon | A molecular-level fire growth parameter[END_REF] suggested an equation where LOI is predicted from flammability parameters measured by PCFC [START_REF] Lyon | Pyrolysis combustion flow calorimetry[END_REF]. Nevertheless, these approaches are suitable to a limited range of polymers and/or lead to meaningless values for some materials. More recently, Lyon et al. [START_REF] Lyon | A molecular-level fire growth parameter[END_REF] develop a Fire Growth Capacity (FGC) from PCFC data. This index involves the specific heat of complete combustion and the temperatures at 5% and 95% conversion of solid to gaseous fuel, that can be assimilated respectively to the ignition temperature and the burning temperature at macroscale. FGC was correlated (R 2 = 0,7) with LOI for a wide range of polymers. Other equations have also been proposed to calculate the pHRR in cone calorimeter from parameters measured at microscale [START_REF] Lyon | Screening flame retardants for plastics using microscale combustion calorimetry[END_REF].

Moreover, a better quantification of FR modes-of-action is needed to assess properly synergism. Indeed, synergism is often claimed but rarely proved. Synergism [START_REF] Schartel | Recommendations on the scientific approach to polymer flame retardancy: part 1 -Scientific terms and methods[END_REF] can be postulated based on performances only or by establishing specific mechanisms. Camino et al. [START_REF] Camino | The chemical kinetics of the polymer combustion allows for inherent fire retardant synergism[END_REF] focused on the cyclic chemical process that controls the polymer combustion, and showed that the maximum synergistic effect should be obtained considering a quantitative ratio of the condensed to gas phase action.

Lewin [START_REF] Lewin | Synergism and catalysis in flame retardancy of polymers[END_REF] and Holdsworth [START_REF] Holdsworth | Novel metal complexes as potential synergists with phosphorus based flame retardants in polyamide 6.6[END_REF] have defined several parameters to evaluate the synergy according to the results obtained in LOI.

This study principally focuses on two approaches proposed to assess the actions of different FR systems in two different fire tests (Limiting Oxygen Index and Cone calorimeter) using flammability data from Pyrolysis-Combustion Flow Calorimetry (PCFC). The goal of the study is to characterize quantitatively the efficiency in condensed phase of FRs, with the aim of clarifying the means of quantitatively evaluating the synergy between FRs. Several examples based on EVA/PE and ABS compounds are discussed. Different modes of action are studied using a halogen FR (DBDPO) which acts in the gas phase action and a phosphorus based compound (APP) which acts in condensed phase in order to evaluate if synergistic effects were obtained [START_REF] Laoutid | New prospects in flame retardant polymer materials: from fundamentals to nanocomposites[END_REF] .

Experimental part

Description of the method

The determination of the limiting oxygen index involves the sustained ignition and the downward propagation of a flame along a vertically orientated sample. Therefore, LOI should be related to the laws which govern these two phenomena and particularly the flame propagation. Quintiere proposed a simplified theory for flame spread rate based on the hypothesis that spread is due to the ignition of a small volume heated by the flame. For thermally thin solids, which can be considered as a reasonable assumption for LOI specimen, the flame spread velocity V f can thus be described by Eq. ( 1) [START_REF] Quintiere | A simplified theory for generalizing results from a radiant panel rate of flame spread apparatus[END_REF] :

V f = ˙ q f δ f e ρC p T ig -T 0 (1)
With ˙ q f the heat flux from the flame, δ f the length heated by the flame, e the sample thickness, ρ the sample density, C p the heat capacity, T ig the ignition temperature, T 0 the ambient temperature. In Eq. ( 1) , V f is supposed to be the flame velocity in free atmosphere, i.e. the oxygen concentration being 21 vol%. Additionally, V f is affected by the variation of oxygen concentration X O 2 . Nasr et al. in their work on confined fire showed that the heat flux from the flame ˙ q f is composed of a convective part which is proportional to X O 2 and of a radiative part which is proportional to X 4 O 2 [START_REF] Nasr | Fuel Mass-Loss Rate Determination in a Confined and Mechanically Ventilated Compartment Fire Using a Global Approach[END_REF] . In the case of LOI with downward flame spread, the radiative part of ˙ q f should be preponderant. Hence, the flame spread velocity should vary as:

V f = k X 4 O 2 (2)
LOI can be considered as the value of oxygen concentration for which the flame spread velocity reaches a critical value V f c , this value being independent of oxygen concentration. From these considerations, it can be proposed that:

LOI = X f ree atm O 2 V f c V f ree at f 1 4 = 21 × ⎛ ⎝ V f c ˙ q f δ f e ρC p ( T ig -T 0 ) ⎞ ⎠ 1 4
(3)

In the present method, it is proposed that some parameters of Eq. ( 3) are directly related to parameters measured in PCFC experiments. Thus, T ig can be related to T max the temperature at maximum pyrolysis rate in PCFC. Moreover, the heat flux from the flame ˙ q f can be written as:

˙

q f = α × HRC × χ × θ LOI (4)
With α an ad hoc parameter (g.K.s -1 .m -2 ), HRC the heat release capacity (J.g -1 .K -1 ) measured in PCFC and assessing the decomposition rate and the dilution of the polymer in condensed phase, χ the combustion efficiency assessing the gas phase action and θ the relative efficiency of FR in condensed phase (barrier effect or endothermic effect, and more generally all phenomena which are not effective in PCFC, i.e. which do affect neither HRC nor T max ).

The combustion efficiency χ was calculated as follows using PCFC and cone calorimeter parameters:

χ = EH C cone HC C PCF C ( 5
)
with

EH C = T H R 1 -f cone residue and H CC = T H R 1 -f PCF C residue ( 6
)
f residue is the mass fraction of residue and THR is the total heat release. χ was calculated by dividing the effective heat of combustion (EHC) of fuel gases in the cone calorimeter test by HCC, Heat of complete combustion, measured at PCFC [ 26 , 27 ]. χ characterizes especially the flame inhibition. It cannot be calculated using PCFC in standard conditions because combustion is assumed to be complete. χ may be calculated using PCFC only if a combustion temperature of isoconversion is known (i.e. the temperature in combustor for which heat of combustion in PCFC is equal to EHC in cone calorimeter). But the data are too scarce to establish a priori this temperature of isoconversion. Then, in this study, χ was calculated using PCFC and cone calorimeter and is assumed to be similar in LOI. From Eqs. ( 3) and (4) , it comes: 9) and [START_REF] Bourbigot | Intumescent ethylene-vinyl acetate copolymer: reaction to fire and mechanistic aspects[END_REF] .

LOI = 21 × β HRC× χ ×θ LOI ( T max -T 0 ) 1 4 (7)
The coefficient β includes all parameters supposed to be constant whatever the materials ( V f c , δ f , e , α) but also the thermophysical properties ( ρ, C p ) because they were shown to have minor influence on LOI.

The method proposed in this study assumes that LOI can be predicted by a so-called Calculated Index (CI) identified in Eq. ( 7) and mainly based on flammability parameters measured in PCFC as shown in equation 8.

CI

= HRC * χ * θ LOI ( T max -T 0 ) (8)
To check the validity of Eq. ( 7) relating LOI and CI, a set of 26 pure polymers (see Table S1 in supporting information) from [START_REF] Lyon | Polymer Flammability[END_REF] was used assuming that θ = 1 since no flame retardant effect in condensed phase was expected. This reference provides the LOI values as well as the heat release capacity, the temperature of ignition that will be assimilated to the temperature at maximum HRR in PCFC and the combustion efficiency. Note that this assumption ( θ = 1) appears reasonable for most of polymers. In case of charring polymers, this assumption may be more questionable.

The experimental data were fitted using Eq. ( 7) . The higher coefficient of determination (R ² = 0.88) was obtained with β = 2 , meaning that

LOI = 25 × ( CI ) -1 4 ( 9
)
The data were also fitted using a power law letting free the exponent. The highest coefficient of determination (R ² = 0,89) was obtained for the following relation: [START_REF] Bourbigot | Intumescent ethylene-vinyl acetate copolymer: reaction to fire and mechanistic aspects[END_REF] Since the exponent in Eq. ( 10) was close to that in Eq. ( 9) deduced from the physical theory, it was chosen to keep Eq. ( 10) that gives the highest coefficient of determination. Fig. 1 highlights that Eq. ( 10) with the new calculated index CI enables to obtain a good prediction of the LOI for pure polymers considering θ = 1.

LOI = 22 . 784 × ( CI ) -0 . 341
As mentioned before, considering thermophysical properties in the denominator of CI:

CI therm = HRC * χ * θ LOI ( T max -T 0 ) * ( ρ * cp )
did not improve the prediction of LOI (R 2 = 0.86), while HRC, χ and T have to be taken into account to maximize the correlation between experimental LOI and CI (with θ = 1). It can be noticed that if HRC is not used in CI the correlation collapses to R 2 = 0.65 whereas without χ or T, R 2 decreases moderately respectively to 0.83 and 0.84. This can be explained by a larger variation range of HRC or combustion efficiency compared to that of thermophysical properties.

Then for flame retardant polymers (i.e. when θ LOI is not a priori equal to 1), θ LOI can be extracted as follows by replacing CI in Eq. [START_REF] Bourbigot | Recent Advances for Intumescent Polymers[END_REF] :

θ LOI = T HRC × χ × LOI 22 . 784 -1 0 , 341 (11)
According to Lyon et al. [START_REF] Lyon | Screening flame retardants for plastics using microscale combustion calorimetry[END_REF] , peak of HRR measured at cone calorimeter can be related to the heat release capacity measured with PCFC by the following equations:

pH RR = H RC × q ext n g ( 12
)
with

n g = h g T p X × θ CC (13)
θ CC can be extracted from Eqs. ( 12) and ( 13) :

θ CC = pHRR × h g × q ext HRC × T p × X ( 14
)
q ext is the external heat flux, h g is the enthalpy of gasification with an average value considered to be 2 kJ/g and T p is the pyrolysis temperature interval (measured as the half width of HRR peak and considered equal to 50 K for pure polymers). Hence, assuming that θ = 1 for pure polymers (see table S2 in supporting information), Fig. 2 shows experimental pHRR versus pHRR calculated from Eq. [START_REF] Cogen | Correlations between pyrolysis combustion flow calorimetry and conventional flammability tests with halogen-free flame retardant polyolefin compounds[END_REF] . Considering that the coefficient of correlation is acceptable (R 2 = 0.84), this equation can be used to extract θ CC which quantifies the condensed phase efficiency in the cone calorimeter test.

In the following, these methods will be applied to flame retardant polymers in order to determine the value of θ in the case it would be different from 1. Hence condensed phase actions can be quantified and possible synergism between FR can be highlighted.

For that purpose, the fire behaviour of several formulations was investigated.

Materials and methods

Polymers used in this study were a mixture of ethylene vinyl acetate with 26wt% of vinyl acetate (EVA Escorene UL00226CC from Exxon Mobil Chemicals) and polyethylene (LDPE PE-019 from Repsol), and acrylonitrile butadiene styrene (ABS Terluran GP22 natural from Ineos). The FRs used were Magnesium hydroxide (MDH -FR20 from ICL-IP), Aluminium hydroxide (ATH -Martinal OL104 from Huber-Martinswerk), Ammonium polyphosphate (APP Exolit AP422 from Clariant) and Decabromodiphenyloxide (DBDPO FR1210 from ICL-IP).

Compounds were extruded using a Clextral BC21 twin-screw extruder (length 90 mm, speed 200 rpm). The obtained pellets were injection moulded (Krauss Maffei 180-CX 50 t) to prepare square specimens (100 × 100 × 4 mm 3 ). The formulations carried out for this study are detailed in the Table 1 .

Each formulation was characterized using three methods. The cone calorimeter experiments were carried out using a Fire Testing Technology (FTT) apparatus. A horizontal sample sheet of 100 × 100 × 4 mm 3 was placed at 25 mm below a conical heater. Samples were tested at a heat flux of 50 kW/m 2 with piloted ignition. All samples were tested twice. The PCFC analyses were performed with a FTT apparatus. The samples were heated from 150 to 900 °C at 1 °C/s in a pyrolyzer under nitrogen flow. Gases produced during pyrolysis were sent into a combustor and burnt at a temperature of 900 °C in presence of 20% of oxygen. In both methods, heat release rate was calculated according to the consumption of oxygen thanks to the Huggett's relation which considers that 1 kg of consumed oxygen during the combustion corresponds to 13.1 MJ of released energy. LOI was measured according to ISO 4589 specifications with a precision of 1% imposed by the accuracy of the machine using a Stanton Redcroft instrument. Tests were performed on 100 × 10 × 4 mm 3 samples that were cut from the injected sheets.

Results and discussion

Fire performances of fr systems

Cone calorimeter

Some HRR curves in cone calorimeter are shown in Fig. 3 . As expected, pHRR is clearly reduced when the content of fillers in the polymer increases. Particularly, in EVA/PE, the incorporation of metallic hydroxides allows increasing time-to-ignition due to the release of water through endothermic decomposition, which cools the surface of the polymer and dilutes the fuel in gas phase [START_REF] Camino | Effect of hydroxides and hydroxycarbonate structure on fire retardant effectiveness and mechanical properties in ethylene-vinyl acetate copolymer[END_REF] . Moreover, MDH seems to be better than ATH to increase time-to ignition, at 60wt% TTI = 133 s for MDH and 98 s for ATH. This phenomenon is maybe due to the release of water, which occurs at low temperature for ATH when the pyrolysis has not started yet. MDH releases water at higher temperature during the first step of degradation of EVA (i.e. release of acetic acid). Percentage of ATH and MDH must be at least 60wt% in order to reach an important improvement in terms of pHRR which becomes lower than 200 kW/m 2 . Nevertheless, it may be noticed that at 60wt% the protective layer formed by dehydrated ATH is effective longer than that formed by dehydrated MDH. This may be related to the thickness of the mineral layer. Indeed, residue from alumina (after ATH decomposition) is around 4mm-thick, while magnesium oxide forms a very thin layer. Combining the two fillers allows to combine the effects: indeed, the presence of 30wt% MDH leads to a high TTI of 120 s while the presence of 30wt% of ATH is enough to reach a significant decrease of HRR and a flaming period higher than 700 s. The second peak at 700 s may be caused by cracking char or when the pyrolysis front reaches the bottom of the sample.

Concerning flame retarded ABS compounds, using APP or DB-DPO at 20wt% and 30wt% content decreases TTI. ABS containing 30wt% and 40wt% of APP shows a curve typical of a material exhibiting a barrier effect with an increase in HRR until an efficient char layer is formed and then a progressive decrease of HRR up to the flame-out.

On the contrary, DBDPO-containing ABS curves are characteristic for a non-charring material with no residue and an increase up to a quasi-static HRR that goes ahead until the pHRR occurs just before the flame out. DBDPO seems to act more efficiently than APP. The combination of the two FRs (20/20wt%) decreases the pHRR to a constant value under 200 kW/m 2 and leads to the disappearance of both the first peak of HRR observed after ignition with APP, and the second peak observed with DBDPO at the end of the test. This behaviour may indicate that DBDPO by its gas phase action hinders the peak of HRR just after ignition while APP promotes the formation of a char that leads to controlled decomposition rate up to the flame-out.

Results obtained with cone calorimeter are summarized in Figs. 4 and5 . These Figures also highlight that combining ATH and MDH do not provide a remarkable effect of synergy. From performances point of view, the combination of 20wt% APP and 20wt% DBDPO may be considered as synergistic, i.e. at the same filler content, the combination of both APP and DBDPO leads to lower pHRR than the incorporation of APP or DBDPO alone.

The thermal and thermo-oxidative decomposition of flame retarded EVA has been widely reported in the literature [START_REF] Holland | The thermal degradation of poly(vinyl acetate) measured by thermal analysis-Fourier transform infrared spectroscopy[END_REF] . It takes place in two steps: the first one corresponds to a deacetylation, with the elimination of acid acetic between 300 °C and 400 °C, the second step is attributed to the degradation of unsaturated products obtained during the first step. The PCFC curves show the two pHRR corresponding respectively to the first and the second steps of the EVA decomposition. PCFC results show a decrease of HRR for each filler ( Fig. 6 ). The reduction in pHRR for all compounds is less distinct than in cone calorimeter test. Moreover, for EVA/PE compounds, the decrease of HRR is proportional to the filler content: ATH and MDH act in a very similar way, with a decrease in HRC of 20%, 40% and more than 50% for respectively 20wt%, 40wt% and 60wt% of fillers, then this decrease is almost exclusively assigned to the dilution of the condensed phase. Results obtained with PCFC are summarized in Figs. 7 and8 . The mixture of ATH and MDH does not show any synergy since HRC of the blends is approximately equal to the mean value of HRCs of the components. PCFC does not consider physical effects like the barrier effect because of the sample size (2 mg). Flame inhibition (from DBDPO) is also not effective in PCFC in standard conditions. Nevertheless, for ABS compounds, the decrease of HRR is lower with APP than with DBDPO. As an example, with a 40wt% filler content, the HRC decreases by 56% with DBDPO while the decrease is only 37% with APP. Concerning the mixture of the two FRs, the decrease in HRC for 40wt%(APP + DBDPO) is close to the formulation with 40wt% of DBDPO.

In LOI, for both polymers and with the different flameretardants, the increase of FR content leads to the increase of LOI. This effect appears exponential for ATH or MDH in EVA/PE and DB-DPO in ABS. Fig. 9 shows that in EVA/PE containing ATH or MDH, LOI values are not significantly different and, once again, no effect of synergy was observed for this parameter when combining ATH and MDH. For ABS, DBDPO is much more efficient than APP. No effect of synergy was observed, i.e. the mixture of both additives does not have a greater effect than each one considered separately. On the contrary, the formulations containing both fillers have the same behaviour than those containing only APP. Hence, the fire spreads slowly and the test ended because of the time of combustion was higher than 3 min for oxygen index (OI) < LOI. On the contrary, for ABS filled with DBDPO only, the behaviour was completely different with a rapid extinction for OI < LOI and a faster flame spread for OI ≥ LOI. Then, the mode of action of APP seems to limit the flame inhibition of DBDPO in LOI test. In our case, it seems that mixing a brominated FR, which acts in gas phase by flame inhibition, and APP that acts in condensed phase as char promoter leads to an antagonistic effect in LOI test.

All main data from cone calorimeter, PCFC and LOI tests are summarized in Table 2 . The combustion efficiency χ is close to one for the majority of the formulations, except for DBDPO compounds with a combustion efficiency around 0.40, even in presence of APP. This low value illustrates the flame inhibition of this halogenated flame retardant. Horrocks et al. [START_REF] Horrocks | Quantification of Zinc Hydroxystannate and Stannate Synergies in Halogen-containing Flame-retardant Polymeric Formulations[END_REF] applied a linear methodology in order to define a synergistic effect, using the following equation:

SE = ( O I F R + S+ P -O I p ) ( ( O I F R + P -O I P ) + ( O I S+ P -O I P ) ) ( 15
)
Where OI FR + P + S is the OI for polymer containing a flame retardant with a synergist. In our case, for the formulations containing two flame retardants, one of the additives is considered as the synergist. Thus Eq. ( 15) can be applied in order to evaluate the potential synergistic effect between ATH and MDH in EVA/PE and between APP and DBDPO in ABS. In the case of EVA/PE filled with 60wt% of ATH and MDH, SE was equal to 2.6, supposing a good synergy between ATH and MDH, while SE was only 1 for 40wt% of ATH and MDH. This confirms the need to use high filler content in order to optimize the FR efficiency. Nevertheless, previously we showed that 60wt% of ATH or MDH have the same impact than 60wt% of [ATH + MDH]. Then the index proposed by Horrocks is not repre-sentative in our case. For flame retarded ABS with 40wt% of APP and DBDPO, SE was equal to 0.9 illustrating an antagonistic effect, as actually observed.

Contrary to LOI, some synergistic effects can be observed in other tests. The most obvious synergy is observed between APP and DBDPO in ABS in the cone calorimeter test since pHRR of the formulations containing the additive mixture is lower than that of compositions containing only one filler at the same total content. Since this synergy was not noticeable in the LOI measurements, it can be concluded that synergies depend on the test used. Note that synergism is only considered from a performance point of view, i.e. without considering the supposed underlying mechanism.

Assessment of FR efficiency in condensed phase

Assuming that LOI and pHRR in cone calorimeter of flameretardant polymers can be predicted according to Eqs. ( 11) and ( 14) , θ was determined for the EVA/PE and ABS flame retarded formulations. Although the Eqs. ( 10) and ( 12) allow assessing LOI and pHRR with an acceptable accuracy for pure polymers, θ is not strictly equal to 1. In order to compare the efficiency of different FR systems in a same polymer (ABS or EVA/PE), reduced values of relative efficiency of heat and mass transfer at the surface θ LOI r or θ cc r are introduced. θ r = θ F R polymer θ pure polymer [START_REF] Van Krevelen | Some basic aspects of flame resistance of polymeric materials[END_REF] θ values are summarized in Table 3 . Uncertainties of θ were evaluated using the variance formula, Eq. ( 17) , and considering the standard deviation of each parameter. Typically, standard deviations are ± 15 J/g.K, ± 2 K, ± 0.03, ± 500 J/g, ± 1% and ± 20 kW/m 2 , respectively for the parameters HRC, T and T p , χ, h g LOI and pHRR.

u 2 c ( y ) = n i =1 df dxi 2 u 2 ( xi ) (17)
With u(x i ) the standard deviation of the parameter x i . The uncertainty of θ CC r is much more important than that of θ LOI r due to the high standard deviation of the enthalpy of gasification (20 0 0 ± 500 J/g according to Lyon et al. [START_REF] Lyon | Screening flame retardants for plastics using microscale combustion calorimetry[END_REF] ).

Assessment of FR efficiency in condensed phase in cone calorimeter

Fig. 10 reports θ CC values for the flame-retarded compounds. It is noteworthy that, despite some large standard deviations, the discrepancy between formulations is much important, therefore conclusions can be drawn.

In EVA/PE blends, no effect in gas phase was observed for ATH or MDH according to combustion efficiency values (close to 1 in all cases). ATH exhibits a higher condensed phase efficiency compared to MDH in EVA/PE compounds. This result is related to the thickness of the mineral residue and is in agreement with the HRR curve shape in cone calorimeter. According to [START_REF] Schartel | Development of fire-retarded materials-Interpretation of cone calorimeter data[END_REF] , the curve shape is representative of an efficient barrier effect in the case of ATH-based formulations (at 60wt% of FR). Additionally, the results confirm that there is no synergism between ATH and MDH.

θ CC r for ABS formulations containing APP decreases (about 0.6 at 40wt% content). This result evidence the charring promotion by the phosphorus FR. On the contrary, θ CC r is close to 1 for DBDPObased formulations, highlighting that FR acts only in gas phase and does not promote any charring. The smaller is θ CC r , the better is the efficiency in condensed phase.

Table 4 shows the residue content r obtained in PCFC and cone calorimeter. On the whole, for EVA/PE formulations, residues are roughly the same in both tests. It may be noticed that the fraction of residue is correlated with θ . This confirms the condensed phase mode of action of those fillers even if not only the char content but also its structure has an effect on its insulating character as barrier layer. For ABS, the residue of APP containing formulations is significantly higher in cone calorimeter test: it means that barrier effect in cone calorimeter leads to incomplete pyrolysis (residue contains not only char but also a fraction of unpyrolyzed polymer). On the contrary, formulations with DPDPO are almost totally consumed whatever the test.

Assessment of FR efficiency in condensed phase in LOI

The influence of each FR system was assessed for all the flame retarded polymers using the calculated index. The values of θ for the LOI test can be deduced from Eq. ( 16) using all the parameters collected at the previous tests. The values are listed in Table 2 .

In EVA/PE compounds, θ LOI r decreases from 1 to 0.33 and 0.32 for 60wt% of ATH and MDH respectively. However, with only 20wt% and 40wt% of fillers, θ LOI r does not decrease. At these filler contents, LOI increases moderately from 21% to 25% (40wt% ATH) and 26% (40wt%MDH), that is mainly due to the decrease in HRC, resulting from the dilution of the condensed phase (i.e. the polymer is partly replaced by the mineral fillers). It is well known that ATH and MDH need to be used at a content equal or greater than 60wt% to be efficient, below this content, the mode of action in condensed phase is not significant in LOI. For the ABS-based formulations containing APP, θ LOI r decreases from 1 to 0.50 with 40wt% fillers. By comparison, DBDPO induces a limited decrease of θ r from 1 to 0.75 while LOI increases more than with APP. As expected, this shows that APP, as a phosphorus flame retardant, has an important condensed phase effect with the formation of a char. On the contrary, DBDPO mainly acts in the vapour phase and traps radicals propagating thermal oxidation in the flame. That is why θ r remains higher for DBDPO compounds while χ is assumed to be significantly reduced (according to cone calorimeter data) ( Fig. 11 ). Furthermore, APP and DBDPO lead to a strong antagonistic effect in condensed phase since Fig. 11 shows an increase of θ to values largely higher than 1. The combustion efficiency remaining low for the APP/DBDPO combination ( χ between 0.39 and 0.46 in cone calorimeter), the low values of LOI would not be related to an inefficient action of DBDPO in gas phase. It may be noticed that this result depends on the assumption that χ which was calculated using PCFC and cone calorimeter is similar in LOI. Further investigation is needed to check this assumption.

Respective contribution of combustion efficiency and phase condensed efficiency

Since DBDPO and APP act through different modes-of-action, the discussion is focused on the ABS formulations in cone calorimeter. The gas phase effect of DBDPO inhibits the combustion whereas the charring promotion from APP traps the fuel in the solid by charring and/or promotes barrier effect through the formation of protective residue layer.

In flame retarded ABS compounds, both modes of action can be active. Fig. 12 plots θ CC r versus χ for ABS formulations. The first diagonal represents the situation where the combustion efficiency is equal to the condensed phase efficiency. Thus, the contribution of the combustion efficiency and the condensed phase efficiency is revealed using such a graph. Hence, for ABS/APP-DBDPO [70/30],

the decrease in θ CC r is 37% while the decrease in χ is 46%. Thus, the respective contribution of condensed phase efficiency can be evaluated to be 37/(37 + 46)x100 = 45% and the part of combustion efficiency is equal to 46/(37 + 46)x100 = 55% .

The decrease in θ CC r corresponds to the condensed phase effect whereas the decrease in χ corresponds to the combustion efficiency illustrating the gas phase mode of action. Hence, "isoefficiency" curves for which χθ = constant can be plotted. χθ is considered as the burning efficiency and isoefficiency curves allow to illustrate the graphical position of all mixtures with the same efficiency. Thus, it appears that ABS filled with 20wt% APP and 20wt% DBDPO has a burning efficiency of 0.25.

For this latter composition, a synergy (i.e. a lower pHRR than corresponding formulations with the same content of only one FR) was highlighted in cone calorimeter. Isoefficiency curves reveal that the burning efficiency is actually lower for the formulations containing both FRs. However, the term "synergy" is debatable. Indeed, Fig. 13 shows the evolution of θ CC r and χ in ABS as function of respectively APP and DBDPO content. This clearly reveals that θ and χ do not decrease linearly with the FR content. In both cases, there is a strong decrease at low filler content and then the FR efficiency seems to reach a plateau from 20wt% for DBDPO and from 30wt% for APP. In these conditions, the rule of mixtures is probably not appropriate to evidence a synergy effect. Obviously, the efficiency of 20wt% DB-DPO in ABS/DBDPO/APP (60/20/20) will be higher than half of the efficiency of 40wt% DBDPO in ABS/DBDPO (60/40) without considering any interaction between DBDPO and APP.

Considering Fig. 13 and assuming an additive contribution of gas phase and condensed phase effects, it can be supposed that the best burning efficiency at the lowest filler content is reached for DBDPO/ABS ratio of 20/80 and APP/ABS ratio of 30/70. Solving this equation system would lead to an optimal ABS/DBDPO/APP composition with the ratio 60/15/25. θ CC r for ABS with 20wt% APP and χ for ABS with 20wt% DB-DPO are equal to 0.58 and 0.40 respectively. The product is equal to 0.23, i.e. very close to the burning efficiency obtained for the formulation containing 20wt% of APP and 20wt% of DBDPO (0.25) which exhibits values of 0.52 and 0.39 for θ CC r and χ respectively.

Based on such analysis, the combination APP/DBDPO seems to act according to a simple additive effect without interactions between both fillers. This analysis illustrates that the synergy is not so obvious between the two additives. This is particularly due to the fact that their separate effect does not vary linearly with their content as shown in Fig. 13 .

Conclusion

In this study, the efficiency of different FR systems in condensed phase was assessed using equations predicting the LOI and pHRR from a couple of parameters. The equation for pHRR in cone calorimeter has been proposed by Lyon while the equation for LOI was obtained from a set of 26 pure polymers. In both cases, the reliability of these equations was checked assuming that the efficiency in condensed phase is equal to 1 for pure polymers.

The values for the efficiency in condensed phase are in agreement with the expected modes-of-action of the FR systems stud- ied. ATH and MDH act in condensed phase but their action is efficient only for high contents (40wt% in cone calorimeter and 60wt% in LOI). APP acts only in condensed phase (low θ and χ close to 1) while DBDPO is an efficient flame inhibitor acting manly in gas phase ( θ close to 1 and low χ).

The method allows a quantitative assessment of the modesof-action by discriminating actions in gas and condensed phases. Isoefficiency curves can be drawn to identify the modes-of-action, from pure flame inhibitor ( θ close to 1 and low χ) to char promoter ( χ close to 1 and low θ ).

It also provides information about synergistic or antagonistic effect. Combining ATH and MDH does not seem to be a suitable strategy. On the contrary, strong negative or positive interactions have been highlighted between APP and DBDPO. Antagonism was pointed out in LOI, through a high θ value (considering that χ is the same in LOI and in cone calorimeter).

Nevertheless, synergism may be falsely postulated because the effect of a flame retardant is not proportional to its content. Indeed, synergism was apparently highlighted in the cone calorimeter test from a performance of view. The combination APP/DBDPO (20/20) shows a low χ (i.e. a high efficiency in gas phase) and an intermediate θ (between the values for DBDPO and APP only). Nevertheless, the efficiency burning is similar to the product of θ for ABS filled with 20wt% APP and χ for ABS filled with 20wt% DBDPO (i.e. considering that both FR act independently, one in condensed phase and the second in gas phase exclusively). Therefore, synergism between APP and DBDPO is questionable.

Fig. 1 .

 1 Fig. 1. LOI versus calculated index: experimental data and Eqs. (9) and (10) .

Fig. 2 .

 2 Fig. 2. Relationships between experimental pHRR measured at 50 kW/m 2 versus calculated pHRR.

Fig. 3 .

 3 Fig. 3. Cone calorimeter curves for flame retarded EVA/PE and ABS at 50 kW/m 2 .

Fig. 4 .

 4 Fig. 4. pHRR in cone calorimetry (50 kW/m 2 ) for EVA-PE filled with ATH (A), MDH (M) and ATH + MDH ( A + M ).

Fig. 5 .

 5 Fig. 5. pHRR in cone calorimetry (50 kW/m 2 ) for ABS filled with APP (AP), DBDPO (D) and APP + DBDPO (AP + D ).

Fig. 6 .

 6 Fig. 6. PCFC curves of the flame retarded EVA/PE and ABS compounds.

Fig. 7 .

 7 Fig. 7. HRC measured by PCFC for EVA-PE filled with ATH and MDH.

Fig. 8 .

 8 Fig. 8. HRC measured by PCFC for ABS filled with APP and DBDPO.

Fig. 9 .

 9 Fig. 9. LOI versus FR content for EVA-PE systems and ABS systems.

Fig. 10 .

 10 Fig. 10. θ CC r versus FR content for EVA-PE compounds and ABS compounds.

Fig. 12 .

 12 Fig. 12. θ CC r versus χ for flame retarded ABS compounds and several isoefficiency curves.

Fig. 13 .

 13 Fig. 13. Evolution of θ CC r and χ according to the content of APP and DBDPO respectively.

Table 1

 1 Detailed mass composition of the compounds studied (in wt%).

		1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
	polymers
	EVA/PE	80	60	40	80	60	40	60	40
	ABS									80	70	60	80	70	60	70	60
	FR
	ATH	20	40	60				20	30
	MDH				20	40	60	20	30
	APP									20	30	40				15	20
	DBDPO												20	30	40	15	20

Table 2

 2 Main results of cone calorimeter, PCFC and LOI tests.

			Cone calorimeter			PCFC				LOI
		FRwt%	pHRR (kW/m ²)	EHC (kJ/g)
	EVA/PE	0%	1407	35.8	0.92	910	38.9	452	38.5	21
	ATH	20%	934	33.2	0.98	668	33.9	457	37.3	23
		40%	397	29.4	0.93	537	31.8	457	34.4	25
		60%	156	23.5	0.93	355	25.1	454	32.4	43
	MDH	20%	1166	33.8	0.95	703	35.6	456	36.6	23
		40%	608	31.0	0.94	524	33.0	456	35.6	26
		60%	340	25.7	0.98	321	26.4	450	39.2	44
	ATH-MDH	40%	613	30.3	0.95	533	31.9	460	35.5	25
		60%	225	25.0	0.97	338	25.7	455	34.4	44
	ABS	0%	1231	30.0	0.86	638	34.8	414	53.2	19
	APP	20%	535	25.0	0.73	523	34.4	415	58.1	23
		30%	432	28.5	0.98	428	29.1	413	56.2	24.5
		40%	385	30.7	1.12	401	27.4	411	56.4	26
	DBDPO	20%	438	11.7	0.40	411	29.3	404	63	25
		30%	323	8.6	0.35	370	24.8	394	61.3	33
		40%	227	7.0	0.34	284	20.8	394	68.1	38
	APP-DBDPO	30%	251	14.1	0.46	347	30.3	397	69.4	25
		40%	179	11.4	0.39	295	29.3	397	70.8	28

χ HRC (J/g.K) HCC (kJ/g) T( °C) T p ( °C) LOI (%)

Table 3 θ

 3 and θ r for the formulations studied.

		Wt% FR	θ LOI	θ LOI r	θ CC	θ CC r
	EVA/PE	/	0.74	1 ± 0.06	1.74	1 ± 0.08
	ATH	20	0.74	1.00 ± 0.06	1.53	0.88 ± 0.23
		40	0.76	1.04 ± 0.1	0.93	0.53 ± 0.14
		60	0.24	0.33 ± 0.03	0.58	0.33 ± 0.09
	MDH	20	0.72	0.98 ± 0.08	1.91	1.09 ± 0.3
		40	0.69	0.93 ± 0.13	1.38	0.79 ± 0.22
		60	0.24	0.32 ± 0.06	1.11	0.64 ± 0.2
	ATH-MDH	40	0.75	1.02 ± 0.09	1.36	0.78 ± 0.2
		60	0.23	0.31 ± 0.03	0.79	0.46 ± 0.12
	ABS	/	1.37	1 ± 0.07	1.79	1 ± 0.08
	APP	20	1.50	1.09 ± 0.09	1.34	0.58 ± 0.16
		30	0.87	0.63 ± 0.04	1.13	0.44 ± 0.11
		40	0.68	0.50 ± 0.04	0.83	0.36 ± 0.1
	DBDPO	20	1.33	0.97 ± 0.05	2.14	1.01 ± 0.27
		30	1.15	0.84 ± 0.05	2.02	0.98 ± 0.26
		40	1.03	0.75 ± 0.08	1.90	0.83 ± 0.25
	APP-DBDPO	30	2.05	1.49 ± 0.09	1.24	0.53 ± 0.14
		40	2.08	1.52 ± 0.15	1.25	0.52 ± 0.15

Table 4

 4 Experimental percentage of residue for flame retarded EVA/PE and ABS in cone calorimeter and PCFC.Fig. 11. θ LOI r versus FR content for EVA-PE compounds and ABS compounds.

		EVA/PE	ATH			MDH			ATH-MDH	ABS	APP			DBDPO			APP-DBDPO
	Wt%	100%	20%	40%	60%	20%	40%	60%	40%	60%	0%	20%	30%	40%	20%	30%	40%	30%	40%
	r CC (%)	0	14	26	42	16	28	43	28	42	0	14	24	30	0.4	4	7	16	20
	r PCFC (%)	0	11	26	39	13	28	41	25	40	0	11.8	8.3	17.4	4.7	4	2.7	13	20.6

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to the work reported in this paper.

Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.polymdegradstab.2021. 109 767.

CRediT authorship contribution statement