Nicolas Bottinelli
email: nicolas.bottinelli@ird.fr

Maeght Jean Luc

R D Pham

Christian Valentin

C Rumpel

V Q Pham

T T Nguyen

D H Lam

A D Nguyen

T M Tran

; J L Maeght

; R D Pham

Q V Pham

; T T Nguyen

; A D Nguyen

R Zaiss

P Jouquet

Anecic earthworms generate more topsoil than they contribute to erosion -Evidence at catchment scale in northern Vietnam

Keywords: tropical soil, bioturbation, earthworm cast

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

Introduction

Earthworms ingest and transform huge amounts of organic and inorganic material, corresponding to up to 30 times their own weight per day [START_REF] Lavelle | Consommation annuelle d'une population naturelle de vers de terre (Millsonia anomala Omodeo, Acanthodrilidae: Oligochetes) dans la savane de Lamto (Côte d'lvoire)[END_REF]. As a result they deposit organic matter rich aggregates at the soil surface (casts), which impact numerous soil functions, such as those associated to the dynamic of carbon and nutrients [START_REF] Van Groenigen | How fertile are earthworm casts? A meta-analysis[END_REF]. High casting activity and stability of these structures are also considered to (1) shape the organization of soil structure [START_REF] Lavelle | Soil aggregation, ecosystem engineers and the C cycle[END_REF], (2) explain the granular organization of A horizons in many ecosystems [START_REF] Jongmans | Soil structure and characteristics of organic matter in two orchards differing in earthworm activity[END_REF], and (3) the burial of archaeological items [START_REF] Darwin | The formation of vegetable mould through the action of worms: With observations on their habits, The Formation of Vegetable Mould Through the Action of Worms: With Observations on their Habits[END_REF][START_REF] Stein | Earthworm activity: a source of potential disturbance of archaeological sediments[END_REF]. Despite their significant influence on main ecological functions in soil, a major obstacle to the quantification and modeling of the influence of earthworms on soil functioning lies in the fact that the amount and dynamics of casts produced at landscape scale remain poorly documented.

Surface earthworm casting activity is highly variable and reaches its record in the tropics with 250 t ha -1 year -1 measured in Nigeria [START_REF] Madge | Field and laboratory studies on activities of 2 species of tropical earthworms[END_REF]. Once deposited on the soil surface, casts are degraded by the rain and/or by the trampling of large animals depending on environmental and ecological factors. For instance, the production and degradation of casts is expected to be higher in tropical than in temperate environments, because of higher temperature and higher precipitations and greater biological activity. Vegetation may protect casts from destruction though raindrop impacts and extend their lifespan [START_REF] Decaëns | Degradation dynamics of surface earthworm casts in grasslands of the eastern plains, of Colombia[END_REF].

However, based on our visual assessment in the field, we hypothesized that the casts' resistance to erosion also depends on soil properties such as their carbon and clay content or on the cast type (granular < past-like < globular) formed in different situations. If not degraded on site, casts can also contribute to the sediment load following water erosion, especially in tropical sloping land [START_REF] Blanchart | Effects of tropical endogeic earthworms on soil erosion[END_REF][START_REF] Jouquet | Utilization of near infrared reflectance spectroscopy (NIRS) to quantify the impact of earthworms on soil and carbon erosion in steep slope ecosystem. A study case in Northern Vietnam[END_REF]. Generally, the impact of erosion on the fate of earthworm casts was estimated from indicators, as for example soil aggregate stability to water or measured directly in the field at small scale on erosion plots under natural or simulated rainfalls [START_REF] Hazelhoff | The exposure of forest soil to erosion by earthworms[END_REF]Jouquet et al., 2008b[START_REF] Jouquet | Influence of earthworms and termites on runoff and erosion in a tropical steep slope fallow in Vietnam: A rainfall simulation experiment[END_REF][START_REF] Jouquet | Surface casting of the tropical Metaphire posthuma increases soil erosion and nitrate leaching in a laboratory experiment[END_REF][START_REF] Le Bayon | Annual variations in earthworm surface-casting activity and soil transport by water runoff under a temperate maize agroecosystem[END_REF]Le Bayon andBinet, 1999, 2001;Podwojewski et al., 2008;[START_REF] Sharpley | Potential role of earthworm casts for the phosphorus enrichment of run-off waters[END_REF]. These approaches are not necessarily representative of the processes taking place at the landscape scale.

Dry and water-stable casts are usually considered to be resistant to raindrop impact.

They increase soil roughness and water infiltration (i.e. reduce water runoff), reduce soil detachment and protect soil from crusting [START_REF] Jouquet | Influence of earthworms and termites on runoff and erosion in a tropical steep slope fallow in Vietnam: A rainfall simulation experiment[END_REF]. In contrast, freshly emitted casts are unstable, prone to dislocation and can increase seal formation and soil loss by erosion [START_REF] Le Bayon | Annual variations in earthworm surface-casting activity and soil transport by water runoff under a temperate maize agroecosystem[END_REF][START_REF] Le Bayon | Earthworm surface casts affect soil erosion by runoff water and phosphorus transfer in a temperate maize crop[END_REF].Therefore, the casting activity could contribute to soil accrual or its erosion. The occurrence of both phenomena may vary in different landscape positions due to microclimatic conditions. This study investigated the balance of both processes to answer the question if earthworms contribute more to soil accrual or its erosion. We addressed this issue by quantifying the impact of earthworm casting activity at the landscape scale in a tropical sloping catchment (50 ha). To this end, we established small plots (0.25 m 2 ) under three vegetation units and quantified (1) casts produced by the anecic earthworm Amynthas adexilis and (2) their impact on water runoff and soil detachment during 2.5 years. Additionally, we assessed cast production in the entire catchment by grid sampling and measured soil transported out of the catchment. We hypothesized that (1) earthworm casting activity is influenced by climatic variables and shows highest activity during the rainy season, and (2) that cast production might counteract soil erosion. [START_REF] Podwojewski | Land-use impacts on surface runoff and soil detachment within agricultural sloping lands in Northern Vietnam[END_REF], derived from the weathering of volcano-sedimentary schists of Mesozoic age. Soils are over 1 m deep but with marked variation in thickness. Their clay content is higher than 50% and they are very porous with a bulk density of 1 g cm -3 . They have a homogenous brown color (10 YR 4/4 to 7.5 YR 4/6), and a weak differentiation. The clay fraction is almost exclusively composed of kaolinite with a low CEC. The catchment is covered by woodlands, secondary forests, meadows and shrublands. Surface earthworm activity in the catchment is high with the earthworm Amynthas adexilis, formerly identified as Amynthas khami (Jouquet et al., 2008a) as the only species. It belongs to the ecological category of anecic earthworms, which show vertical borrowing activity with deposition of earthworm casts at the soil surface. The casts of Amynthas adexilis reach 20 cm length and are produced through an accumulation of fecal pellets deposited the one on the other (Fig. 1).

Methods

Study site

They are always enriched in organic C and can be more water stable than the surrounding soil depending on their degradation stage [START_REF] Bottinelli | To what extent do ageing and soil properties influence Amynthas khami cast properties? Evidence from a small watershed in northern Vietnam[END_REF][START_REF] Mer | Age matters: Dynamics of earthworm casts and burrows produced by the anecic Amynthas khami and their effects on soil water infiltration[END_REF]. In addition, A. adexilis digs vertical burrows that are connected to the surface and can act as preferential pathways for water flow (Le [START_REF] Mer | Age matters: Dynamics of earthworm casts and burrows produced by the anecic Amynthas khami and their effects on soil water infiltration[END_REF].

Earthworm cast production and soil detachment at the plot scale

Eight plots of 0.25 m 2 were set up from July 2016 to December 2018 in each of the three representative vegetation units (woodland, shrubland and meadow) to measure the casting activity of A. adexilis and its impact on runoff and soil detachment. The plot location was chosen randomly in each vegetation unit to account for the field variability of earthworm casting activity. The mean slope was 24, 17 and 13 degrees in woodland, shrubland and meadow, respectively. In each vegetation unit, four plots were kept untouched, while four other plots were used to quantify the weekly production of casts. All earthworm casts were removed from plots before the beginning of the experiment.

Rigid metal frames (50 cm height) were inserted to a depth of 5 cm and used to delimit plots.

A 30 L bucket was put at the outlet of each plot to collect runoff water and sediment. Soil height was measured twice (July 2016 and December 2018) on 110 regular points using a laser distance meter. Soil detachment and water runoff were collected five times at the end of the rainy season in 2018 (18 of August, 1, 22, 29 of September and 27 of October). The total volume of runoff was calculated from the measurement of water height in the collecting bucket. The runoff coefficient corresponded to the ratio of the runoff divided by rainfall. The concentration of soil particles was measured from a 500 mL aliquot, which was filtered and weighted after air-drying. We calculated soil detachment as the product of the sediment concentration in the collected samples and the runoff volume. Cast production was measured through the sampling and weighing of casts after oven drying at 105°C during 24 h. At the end of the experiment (December 2018), all remaining casts were collected from the untouched plots and weighed in the field. Subsamples were taken to the laboratory to determine their moisture after drying them in an oven at 105°C during 24 h gravimetrically.

The sum of cast mass collected every week minus the mass of casts collected after 2.5 years in the untouched plots was used to assess the mass of casts transformed into soil or washed away out of the plots by water runoff.

Earthworm cast production and soil erosion at the catchment scale

The sampling took place in February 2017, during the dry season because of the difficulty of walking in the catchment during the rainy season. 195 points were sampled on a regular grid (50 x 50 m) covering the entire catchment. At each point, on three replicates of 1 m 2 , casts were classified into three groups (i) active (standing and wet); (ii) non-active (stand and dry) and (iii) degraded (dry and broken). The number of active and non-active casts was quantified and the weight of the three categories was determined in the field and corrected for moisture, which was determined on subsamples dried in the oven at 105°C during 24h.

Soil erosion at the catchment scale was measured at the outlet of the catchment using a Vnotch weir. Suspended sediments were collected using an automatic water sampler. The automatic sampler was triggered by the water level recorder to collect water after every 2-cm water level change during flood rising and every 5-cm water level change during flood recession. The concentration of suspended particulate matter in each sample was measured after filtration and evaporation at 105 °C for 48 h. Bed particulate matter was determined from the sediments that were retained in the stilling basins of the weir. Monthly, or if the basin was full, the volume of deposited sediments was measured considering a density of 1 g cm -3 for soft sediments and of 2.65 g cm -3 for stones. The annual erosion rate represented the cumulative suspended and bedload sediments divided by the catchment area.

An automatic weather station measuring on a one-minute basis the temperature and precipitation was located at the bottom of the catchment. Soil temperature sensors were buried at 10 cm depth and watermark sensors were buried at 5, 10, 30, 50, 100, 150 cm depth to continuously monitor soil temperature and soil water tension (0-250 cm) in the three vegetation units.

Statistical analyses

One-way ANOVA and repeated measures ANOVA were performed to test the effect of the vegetation unit and the sampling time on the production of casts in the 0.25 m 2 plots, respectively. T-test was used to test the effect of presence or absence of earthworm casts on water runoff and soil detachment in 0.25 m 2 plots for each vegetation unit. Prior to running ANOVA and T-test, data were tested for homogeneity of variances and normality.

Differences among treatments were declared significant at the 0.05 probability level.

Between-class analysis (BCA) and Monte Carlo permutation tests (1000 permutations) was carried out to discriminate the three vegetation units in 2017 and 2018 in function of soil and air temperature and water tension. BCA measures the amount of variance restricted to the grouping factor as a percentage of the total inertia [START_REF] Dray | Revisiting guerry's data: introducing spatial constraints in multivariate analysis[END_REF]. Random forest regression was used to identify climatic variables influencing cast production in each vegetation unit. Random forest is a nonparametric method, which consists of a large number of individual tree models trained from bootstrap samples of the data [START_REF] Breiman | Random forests[END_REF]. The results of all individual trees are aggregated to make a single prediction. The variables tested in were total precipitations, precipitations > 20 mm h -1 , air temperature, soil temperature and soil water tension from 10 to 150 cm depth. The dataset for each vegetation unit represented 115 measurements of cast production and was randomly split into calibration dataset (80%) and validation dataset (20%). The optimization of parameters for the final models was based on minimizing the root mean square error (RMSE) between the measured and the estimated value of the production of casts for the calibration dataset based on 10-fold repeated crossvalidation with five repetitions. The number of variables for each tree (mtry) and the number of trees in the forest (ntree) are two user-defined parameters which need to be optimized. The number of trees (ntree) was 1000 in this study. Mtry was identified as those returning the lowest RMSE by iterating mtry values from 1 to 11 representing the number of predictors.

The relative importance of variables was estimated from the mean decrease in predictive accuracy. The model performances for prediction of the production of casts on the test dataset were evaluated using coefficient of determination (R 2 ), root mean square error of prediction (RMSE) and residual predictive deviation (RPD). The RPD denotes the ratio of the standard deviation of the measured cast production to the RMSE calculated between the measured and predicted cast production.

The relationship between the number of casts collected in the three vegetation units after one week of production in February 2017 and the cumulative mass of casts produced for 2017 for each plot was calculated. The relationship was used to predict the mass of casts produced in 2017 at every grid sampling points using the number of active casts collected during the sampling of earthworm casts at catchment scale in February 2017. All statistical calculations and plots were carried out using R software using 'caret', 'car', 'agricolae' and 'ggplot2' packages.

Results

Earthworm cast production at the plot scale and consequence on soil detachment

The production of casts at the plot scale varied from 16 to 220 t ha -1 year -1 and decreased in the order meadow>woodland>shrubland (p < 0.01) (Fig. 2). The production of casts was higher in 2017 than 2018 in woodland and meadow (p < 0.01 for in both cases), whereas the opposite was found in shrubland (p < 0.05). In shrubland and meadow, the production followed a periodicity of 6 months with lowest values occurring in December and July. In woodland, cast production followed a periodicity of 12 months with lowest values in December. We used random forest models to investigate the climatic variables influencing the cast production at plot scale. These models did not perform well at predicting the mass of casts with ratio of prediction deviation (RPD) of 1.2, 1.1 and 1.2 for the woodland, shrubland and meadow, respectively (Fig. 3a). For the shrubland and meadow, the most important variables (Fig. 3b) explaining the production of casts were related to soil water tension (at 150 and 5 cm depth in shrubland and meadow, respectively). Conversely in the woodland, the most important variable was soil temperature.

The production of casts led to the formation of a new horizon (Fig. 4a) reaching up to 6.5 cm after 2.5 years (Fig. 4b). The degradation of casts measured for 2.5 years represented between 2 to 34 kg m -2 (Fig. 4c).

For the five rain events combined (152, 207, 80, 69 and 148 mm of precipitations), runoff coefficient (Fig. 5a) and soil detachment (Fig. 5b) were greatest in the absence of casts. However, results were only statistically significant for the soil detachment in meadow.

Earthworm cast production and soil erosion at the catchment scale

Cast production in the entire catchment was evaluated in February 2017. The mass of casts found in the catchment varied from 0 to 6.5 kg m -2 and represented on average 7.5 t ha -1 (Fig. 6). 80% of the casts showed a degraded aspect. The number of active and non-active casts varied between 0-13 and 0-14 per m 2 , respectively. The number of casts collected in 0.25 m 2 plots after one week of production in February 2017 varied between 0 to 20 per m -2 and was positively related (R 2 = 0.90) to the total mass of casts collected during the year (Fig. 7).

Using this relationship for assessing the annual mass of casts produced in the entire catchment with the number of active casts, we predicted that the cast production amounted on average to 36 t ha -1 of casts in 2017 (from 0.2 to 1.7 cm of soil using a bulk density of 1 g cm -3 ). This is 5 times higher than the cast mass observed on average at catchment scale (see above) and 12 times higher than the annual erosion rate measured at the outlet of the catchment, which represented only 3 t ha -1 .

Discussion

Earthworm cast production varies in space and time

The quantity of surface casts produced in the three vegetation units during the years 2017 and 2018 varied from 16 to 220 t ha -1 at plot scale (i.e., between 0.2 to 2.2 cm of soil considering the soil bulk density of 1 g cm -3 ), which is of the same order of magnitude as was found in other tropical regions [START_REF] Hauser | Effect of land use change, cropping systems and soil type on earthworm cast production in West and Central Africa[END_REF][START_REF] Madge | Field and laboratory studies on activities of 2 species of tropical earthworms[END_REF]. The most important factor explaining this variability is the vegetation. Differences are most likely linked to the abundance of earthworms and not their activity since soil properties (data not shown) and soil climate were relatively similar between vegetation units. Indeed, assuming that the number of fresh casts quantified every week is an estimation of the number of individuals per plot [START_REF] Jiménez | Earthworm communities in native savannas and man-made pastures of the Eastern Plains of Colombia[END_REF], the density of A. adexilis may be different between the three vegetation units (Fig. 7). Temperature and rainfall intensity changes throughout the year had strong impacts on the production of casts. In the three vegetation units, decreasing cast production was found in December, when the air temperature and precipitation were the lowest. This is in agreement with other studies, which found relationships between cast production and soil temperature [START_REF] Whalen | Quantifying surface and subsurface cast[END_REF] and moisture (Zaller and Arnone III, 1997).

In the shrubland and meadow, a decrease of activity was also recorded in July, when the air temperature and rainfall were highest. Since the pedoclimatic conditions were similar for the three vegetation units during the wet season, discrepancies might be linked to microclimatic variations affecting the degradation casts at the three sites. While earthworm casts after drying exhibit high aggregate stability as compared to physical soil aggregates (Jouquet et a., 2008), recently formed casts are relatively unstable [START_REF] Mer | Age matters: Dynamics of earthworm casts and burrows produced by the anecic Amynthas khami and their effects on soil water infiltration[END_REF] and may be degraded by rainfall impact. In fact, casts were fully exposed to high intensity precipitations in shrubland and meadow, whereas in woodland they were partly protected by the tree canopy. From 2017 to 2018, we measured lower cast production in woodland (by 50%) and meadow (by 20%).

This may be partly explained by cast disintegration by rainfall impact, since the precipitations were higher by 454 mm in 2018 as compared to 2017. Moreover, contrasting impacts of precipitation events on earthworm abundance and activity under the three vegetation units cannot be excluded.

Surprisingly, the production of casts in shrubland increased by 100% from 2017 to 2018. This increase may result from the large spatial variability of earthworm activity in shrubland and the small size plots for the quantification of casts. In 2017, only half of the plots were colonised by earthworms, whereas in 2018 casts were found on the four plots. In this study we consequently showed that rainfall may have positive impacts and negative impacts on the production of casts depending on the vegetation. This combined effect might explain the relatively low accuracy of the random forest models to predict the production of casts. We suggest that the accuracy of the models may be improved in two ways: by measuring other variables (i) linked to the activity of earthworms, such as plants biomass and litter production, which in turn affect earthworm food supply and (ii) linked to the occurrence of splash erosion (e.g. the intensity of the precipitation and the size of the drops in each vegetation units). In other studies, mesh protection was used to measure only the production of casts by protecting the casts from rain splash [START_REF] Hauser | A continuous sampling technique to estimate surface cast production of the tropical earthworm Hyperiodrilus africanus[END_REF]. However, this method also influenced the measured variable (Le Bayon and Binet, 1999) by maintaining high soil moisture and most likely changing the amount of food supply for earthworms.

Earthworms generate more topsoil than they contribute to erosion

At the plot scale, the accumulation of casts for 2.5 years formed a vermic horizon of up to 6.5 cm height. This value is much higher than results found in England [START_REF] Butt | Action of earthworms on flint burial-A return to Darwin's estate[END_REF][START_REF] Darwin | The formation of vegetable mould through the action of worms: With observations on their habits, The Formation of Vegetable Mould Through the Action of Worms: With Observations on their Habits[END_REF] showing that the rate of flint burial by the surface casting activity of anecic earthworms varied from 0.21 to 0.96 cm year -1 . We assessed that the degradation of casts accounted between 2 and 34 kg m -2 for the 2.5 years of experiment. Comparing these figures with the amount of cast transported as sediment in water runoff (i.e., maximum of 100 g m -2 for four months during the rainy season in 2018), it is clear that casts once degraded were not transported by runoff out of the plots. It is therefore surprising that despite the huge production of casts and obvious cast degradation, the amount of sediments in the runoff water was similar in plots with and without casts. Thus, we hypothesize that high density of water stable casts acted as physical barriers toward the runoff, which promoted the redistribution of eroded casts on a very short distance through successive deposition/suspension processes. Our results were also in agreement with studies carried out at the same study site in 1m 2 plots showing a negative relationship between the density of earthworm casts and soil detachment [START_REF] Jouquet | Influence of earthworms and termites on runoff and erosion in a tropical steep slope fallow in Vietnam: A rainfall simulation experiment[END_REF][START_REF] Podwojewski | Land-use impacts on surface runoff and soil detachment within agricultural sloping lands in Northern Vietnam[END_REF].

At the catchment scale, the mass of casts reached at several locations large amounts of up to 6.5 kg m -2 . This led to the formation of a vermic horizon as it was observed in the 0.25 m 2 plots under all three vegetation units. However, the high percentage of broken casts covering the soil indicated that casts were rapidly degraded. Several factors might be involved such as the raindrop impacts, the trampling of buffalo and the activity of soil organisms. We assessed by a linear relationship, that for 2017, the amount of casts produced in the catchment was 5 times higher than what we measured during the survey in February in 2017. This indicates that up to 79% of all cast produced may be degraded within one year. If all the degraded materials were washed away contributing to material leaving the catchment, the annual erosion rate would have been much higher than 3 t ha -1 . Therefore, we assume that the material composing earthworm casts remained within the catchment after their degradation.

Our results provide evidence that earthworms are strongly involved in redistribution of soil material in this tropical catchment and do not promote soil erosion. In tropical catchments with high bioturbation, earthworms' influence on soil fertility and climate change mitigation may thus be considerable through the formation of organic matter rich horizons. The depth of vermic horizons occurring throughout the catchment may be depending on the landscape position. While the earthworm casting activity may be influenced to some extent by weather conditions, vegetation cover could control the accumulation of the vermic horizon because it may influence casting activity by impacting earthworm abundance and their activity.

Moreover, vegetation also influenced the disintegration of casts after their deposition.

Therefore, the occurrence of vegetation types and soil accrual may be related with higher soil formation rates occurring under meadow as compared to shrubland and woodland. Vegetation type is thus important to consider at the landscape scale in addition to (micro-) pedoclimate when assessing the earthworms' impact on soil accrual.

production by earthworms under controlled laboratory conditions. Biol. Fertil. Soils 39, 287-291. https://doi.org/https://doi.org/10.1007/s00374-003-0715-1 Zaller, J.G., Arnone III, J.A., 1997. Activity of surface-casting earthworms in a calcareous grassland under elevated atmospheric CO2. Oecologia 111,[249][250][251][252][253][254] https://doi.org/10.1007/PL00008817 The white square represented the average mass of the different types of casts. Abbreviation: wet, stand cast in production; dry, stand cast abandoned but not degraded; broken, cast not stand; total, sum of all the types of casts.

Figure captions

  This study was carried out from 2016 to 2018 in a tropical catchment located in Dong Cao village in the North of Vietnam (20° 57'N, 105° 29'E). The site belongs to the M-Tropics long-term observatory. The Dong Cao catchment with an area of 46 ha is located on sloping land with an average slope of 40%, locally exceeding 100%. Annual rainfall varied from 1770 mm in 2017 to 2224 mm in 2018. The annual air temperature was 23°C in 2017 and 2018. The dry season lasting from November to April and was characterized by the lowest precipitation (222 mm) and air temperature (19 ºC). The wet season lasting from May to October and was characterized by the greatest precipitation (1688 mm in 2017 and 1922 mm in 2018) and air temperature (27 ºC). The dominant soil type in the catchment is Acrisol

Figure 1 :

 1 Figure 1: Accumulation of globular earthworm casts on soil surface in the woodland.

Figure 2 :

 2 Figure 2: Weekly cast production measured in plots of 0.25 m 2 in woodland, shrubland and

Figure 3 :

 3 Figure 3: (a) Predicted vs. measured weekly production of earthworm casts of validation

Figure 4 :

 4 Figure 4: (a) photos depicting the formation of a new soil horizon caused by the accumulation of earthworm casts after 2.5 years in plots of 0.25 m 2 in woodland (W), shrubland (S) and

Figure 5 :

 5 Figure 5: Box plots presenting (a) runoff coefficient and (b) cumulative soil detachment measured in plots of 0.25 m 2 with earthworm casts and without cast from the 18 of August to

Figure 6 :

 6 Figure 6: Mass of earthworm casts sampled in 195 points covering the catchment of 50 ha.

Figure 7 :

 7 Figure 7: linear relationship between the number of earthworm casts collected after 1 week of production in plots of 0.25 m 2 in February 2017 in woodland (circle), shrubland (triangle) and

Acknowledgments

This project was financially supported by the CNRS/INSU (VINAWORM) research program under the framework of the EC2CO program. Dong Cao catchment is part of Multiscale Tropical Catchments (M-Tropics) project (https://mtropics.obs-mip.fr/) supported by the French national research institute for sustainable development (IRD). We are very grateful to Kim Yen Do for her invaluable help in organizing the experiment. We also acknowledge the technical assistance of the farmers in Dong Cao village for their participation in the sampling of soil samples.


Estelle Camizuli

Magali Rossi
email: magali.rossi@univ-smb.fr

Dominique Gasquet

years of mining activity in the northern French Alps

Keywords: Streambed sediments, Trace metal anomalies, Ancient mines, northern French Alps, Geostatistical methods. 1. Introduction

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

2016; Camizuli et al., 2018), when not based on atmospheric fluxes and TM deposition in ice or lake sediments [START_REF] Rosman | Lead from Carthaginian and Roman Spanish Mines Isotopically Identified in Greenland Ice Dated from 600 B.C. to 300 A.D. †[END_REF][START_REF] West | Heavy Metals in Holocene Peats from South West England: Detecting Mining Impacts and Atmospheric Pollution[END_REF]Arnaud et al., 2005;[START_REF] Giguet-Covex | Changes in erosion patterns during the Holocene in a currently treeless subalpine catchment inferred from lake sediment geochemistry (Lake Anterne, 2063 m a.s.l., NW French Alps): The role of climate and human activities[END_REF]. Rather scarce studies investigate long-term (i.e., from a hundred years to the Middle Ages or Antiquity) TMs dispersal from streambed sediment data [START_REF] Garcia-Ordiales | Stream bottom sediments as a means to assess metal contamination in the historic mining district of Almadén (Spain)[END_REF][START_REF] Monna | Tracking archaeological and historical mines using mineral prospectivity mapping[END_REF]. Persistence of trace metals in streambed sediments may reach potentially toxic levels for biota over time (Camizuli et al., 2014b;[START_REF] Garcia-Ordiales | Stream bottom sediments as a means to assess metal contamination in the historic mining district of Almadén (Spain)[END_REF][START_REF] Lecce | Floodplain storage of sediment contaminated by mercury and copper from historic gold mining at Gold Hill, North Carolina, USA[END_REF][START_REF] Beane | Abandoned metal mines and their impact on receiving waters: A case study from Southwest England[END_REF][START_REF] Gutiérrez | Abandoned Pb-Zn mining wastes and their mobility as proxy to toxicity: A review[END_REF][START_REF] Elbaz-Poulichet | The environmental legacy of historic Pb-Zn-Ag-Au mining in river basins of the southern edge of the Massif Central (France)[END_REF]Camizuli et al., 2018).

This paper investigates trace metals dispersal in the northern French Alps, where hundreds of small Fe, Cu, Pb-Zn and other metal deposits [START_REF] Meloux | Alpes-Nord. Bilan synthétique de dix années de prospection[END_REF][START_REF] Durand | The long-term interaction of mine tailings with soils and the wider environment: Examples from Mont Chemin, Switzerland[END_REF] have been mined from the Bronze Age to the 20th century [START_REF] Schibler | Miners and mining in the Late Bronze Age: a multidisciplinary study from Austria[END_REF][START_REF] Moulin | Quatre années de prospection sur les extractions de cuivre de l'âge du Bronze ancien dans le massif des Rousses en Oisans (Isère et Savoie, France)[END_REF][START_REF] Py | Interdisciplinary characterisation and environmental imprints of mining and forestry in the upper Durance valley (France) during the Holocene[END_REF]. Anthropogenic lead records in high-altitude lake sediments (Arnaud et al., 2005;[START_REF] Giguet-Covex | Changes in erosion patterns during the Holocene in a currently treeless subalpine catchment inferred from lake sediment geochemistry (Lake Anterne, 2063 m a.s.l., NW French Alps): The role of climate and human activities[END_REF] and archeological and historical investigations [START_REF] Meloux | Alpes-Nord. Bilan synthétique de dix années de prospection[END_REF]Bailly-Maître, 2019) indicate peak activity during the Middle Ages and the 19th century. Apart from Macôt-La Plagne and St-Georges d'Hurtières mines, which closed in the 1970s and the 1930s respectively, mining activity essentially stopped in mid-19th century or during the Middle Ages. Therefore, the northern French Alps is a good area to investigate the long-term dispersal of trace metals into the environment resulting from mining-activity, over several centuries. Several statistical methods derived from modern mining prospecting and ecological management plan [START_REF] Garcia-Ordiales | Stream bottom sediments as a means to assess metal contamination in the historic mining district of Almadén (Spain)[END_REF]Monna et al., 2014;Afzal et al., 2016;Agharezaei and Hezarkhani, 2016;[START_REF] Hoang | Identification of Spatial Distribution of Geochemical Anomalies Based on GIS and C-A Fractal Model -A Case Study of Jiurui Copper Mining Area[END_REF][START_REF] Mortazavi | Heavy metals assessment of surface sediments in Mighan wetland using the Sediment Quality Index[END_REF] were applied to the streambed sediment database provided by the Bureau de Recherches Géologiques et Minières (BRGM; French Geological Survey). Regional prospectivity maps were produced using Exploratory Data Analysis (EDA) and a fractal-based model to highlight As, Cu, Ni, Pb, and Zn anomalies. Sediment Quality Index (SQI, [START_REF] Marvin | Application of a Sediment Quality Index to the Lower Laurentian Great Lakes[END_REF] was also computed from seven TM concentrations (As, Cd, Cr, Cu, Ni, Pb, and Zn) to assess the overall quality of streambed sediments based on ecotoxicological thresholds. In order to investigate the local impact of ancient mining activity, this paper focuses on three major Pb-Ag mining districts that were mined at different historical periods from the Middle Ages to Modern Times: (i) Brandes, which was an important medieval Pb-Ag exploitation, (ii) Peisey-Nancroix, which was mined during the 17th-19th c. AD, and (iii) Macôt-La Plagne, which was a huge modern Pb-Ag mining area before becoming a major ski resort in the 1970s [START_REF] Durand | The long-term interaction of mine tailings with soils and the wider environment: Examples from Mont Chemin, Switzerland[END_REF][START_REF] Durand | Un siècle dans les mines de Savoie, Gap, Challes-les-Eaux[END_REF]. In the French Alps, silver was mined since at least the Middle Ages as it was used for mint [START_REF] Bailly-Maître | Un coron du Moyen-Age, le site de Brandes-en-Oisans (Huez, Isère)[END_REF], 2019). Pb-Ag mines are thus good candidates to study the local impact of past mining activity from streambed sediments analyses considering Pb is a stable element that remains for a long time in the environment, and is associated with other TMs such as arsenic, antimony or copper.

Geological and historical contexts 2.1. Geological setting

This study focuses on the northern French Alps, in an area of c. 4700 km² (Fig. 1a), with the altitude ranging from 300 m a.s.l. up to more than 3000 m a.s.l., resulting in a complex hydrographical network.

The current structure of the Alpine belt reflects its complex history. The external Alps represent the old European margin whereas the internal Alps correspond to the old Apulian margin. The internal and external Alps, separated each other by the major Penninic Front, consist of magmatic and metamorphic rocks formed during the Variscan orogeny that were exhumed during the Alpine collision, and extruded into their Mesozoic sedimentary cover during the Alpine orogeny [START_REF] Coward | Alpine tectonics[END_REF]Agard and Lemoine, 2003;[START_REF] Marthaler | The Alps and our Planet : The African Matterhorn, a geological story[END_REF]. The metamorphic Variscan basement is mainly composed of various gneisses, micaschists, migmatites and mafic rocks (amphibolites and metagabbros) intruded by Variscan and late-Variscan granitoids. These rocks essentially crop out in the external crystalline massifs (Ecrins-Pelvoux, Belledonne and Lauzière) of the external Alps. In the internal Alps, Variscan micaschists, Carboniferous schists and Cretaceous schists are abundant. The Mesozoic sedimentary cover is mainly composed of Triassic dolomites (internal Alps) and Jurassic and Cretaceous limestones and marls (external Alps). These lithologies correspond to the carbonates in Figure 1. Eocene flysches are the most abundant detrital rocks. Superficial deposits (Quaternary tills and alluviums) are found in the main valleys.

The complex Variscan and Alpine orogeneses were characterized by major fluid flows that were locally responsible for metal concentrations that have been mined since Antiquity [START_REF] Py | Interdisciplinary characterisation and environmental imprints of mining and forestry in the upper Durance valley (France) during the Holocene[END_REF]. However, the distribution and size of ore deposits is very heterogeneous depending on the geological setting. The Eastern Alps hosted some world class Fe-Cu-Sn deposits that have been mined since the Bronze Age [START_REF] Schibler | Miners and mining in the Late Bronze Age: a multidisciplinary study from Austria[END_REF]; some being still operated as open pits. On the opposite, the French Alps are characterized by numerous and disseminated small size ore deposits; about 350 ore deposits have been reported by the BRGM, but this number has been reassessed by [START_REF] Durand | The long-term interaction of mine tailings with soils and the wider environment: Examples from Mont Chemin, Switzerland[END_REF][START_REF] Durand | Un siècle dans les mines de Savoie, Gap, Challes-les-Eaux[END_REF].

Base-metal deposits (Fe, Cu, Pb-Zn-Ag) are the most common in the French Alps. With the exception of some sedimentary iron ores, most deposits occur as vein-type and stratiform mineralization that were operated underground. The distribution and size of ore deposits are clearly depending on the geological setting. Most ore deposits are located in the external Alps, in the External Crystalline Massifs and in the Permo-Triassic sedimentary formations. Some ore deposits are also located in the internal zone and are mainly hosted in the carboniferous schists (Fig. 1a; [START_REF] Meloux | Alpes-Nord. Bilan synthétique de dix années de prospection[END_REF][START_REF] Rossi | Panorama géologique des exploitations minières dans les Alpes françaises[END_REF][START_REF] Gasquet | Environmental Impacts and Management of Ancient Abandoned Mines in the French Alps[END_REF].

Mining activity in the northern French Alps

The mining industry was quite important in the French Alps. Since the Bronze Age, mining activity contributed to human settlement and economic development [START_REF] Py | Interdisciplinary characterisation and environmental imprints of mining and forestry in the upper Durance valley (France) during the Holocene[END_REF]. If many alpine mines were operated since Antiquity, mining activity reached its maximum during the Middle Ages and the 19th century, as evidenced by Pb records in lake sediments and peat bogs (Arnaud et al., 2005;[START_REF] Giguet-Covex | Changes in erosion patterns during the Holocene in a currently treeless subalpine catchment inferred from lake sediment geochemistry (Lake Anterne, 2063 m a.s.l., NW French Alps): The role of climate and human activities[END_REF][START_REF] Py | Interdisciplinary characterisation and environmental imprints of mining and forestry in the upper Durance valley (France) during the Holocene[END_REF]. The activity started to decline at the end of the 19th century [START_REF] Gasquet | Environmental Impacts and Management of Ancient Abandoned Mines in the French Alps[END_REF] due to the discovery of giant ore deposits worldwide. Metallurgical industries developed since the 19th century, due to the mining and processing industry and to the development of hydroelectricity. These industries are still in activity in the alpine valleys.

In the study area, 137 mining sites have been reported by the BRGM (Fig. 1b). Data and complementary information were freely available online at http://sigminesfrance.brgm.fr (last accessed 19 July 2019; now moved to http://infoterre.brgm.fr, Gisements, gîtes et indices and Mines-Substances principales databases). Most mines were small-scale underground exploitations, sometimes located above 1500 m, a.s.l. A few districts were more important, such as the St-Georges-d'Hurtières district which provided 1.5 Mt of siderite and more than 10 000 t of chalcopyrite between the Middle Ages and the 20th century. This paper focuses on trace metals dispersal around the three main Pb-Zn-Ag districts of the northern French Alps that were mined at different time periods from the Middle Ages to Modern Times. The Macôt-La Plagne mine, located at 1800 m a.s.l., has been one of the most important Pb-Ag mine in France: 126 600 t Pb and 350 t Ag were extracted from 1810 to 1973 [START_REF] Meloux | Alpes-Nord. Bilan synthétique de dix années de prospection[END_REF]. Mining activity was stopped due to the discovery of huge deposits worldwide making extraction no more profitable, and to the development of a major ski resort at La Plagne. Historical and archaeological records indicate that Pb-Ag may have been episodically extracted in Macôt-La Plagne since Roman Period [START_REF] Gimard | Les travaux romains des mines de plomb argentifère de Macôt-la-Plagne[END_REF]. The Pb-Ag mine of Peisey-Nancroix, located ~10 km east of Macôt-La Plagne at 1500 m a.s.l., was mined from 1644 to 1865, with a peak activity in the 18th century [START_REF] Nicolas | La Savoie au XVIIIe siècle, tome 2, Inflexions au siècle des Lumières[END_REF]: 22 000 t of Pb and 53 t of Ag were extracted during this period [START_REF] Meloux | Alpes-Nord. Bilan synthétique de dix années de prospection[END_REF]. Episodic mining activity occurred during the Middle Ages [START_REF] Durand | Un siècle dans les mines de Savoie, Gap, Challes-les-Eaux[END_REF]. Opposite to Macôt-La Plagne, the site has been unmodified since mine closure in 1865 and waste deposits (tailings, slags and washing residues) have been preserved on site. The Brandes (Huez, Oisans) district, located at high altitude (between 1800 and 2700 m a.s.l) was the most important medieval Pb-Ag exploitation in the French Alps. Archeological surveys show rational mining activity (metallurgical wastes, workshops, hydraulic equipments) from the 11th to the 17th centuries [START_REF] Bailly-Maître | Un coron du Moyen-Age, le site de Brandes-en-Oisans (Huez, Isère)[END_REF], 2019). There is no estimation of the extracted ore volume since the 11th century, but 900 000 tons of dumps (mainly barite) are bordering the Brandes plateau (pers. com. M.Ch. Bailly-Maître).

Materials and Methods

Geochemical data

In the study area, the BRGM collected 6892 samples of streambed sediments (Fig. 1c). The analyses were undertaken during the 1980s through an inventory of mineral substances at the national scale [START_REF] Lambert | Les données géochimiques et alluvionnaires de l'Inventaire minier du territoire national. Constitution d'une base de données exhaustive[END_REF]. Sediments were sieved at 125 µm, dried, and analyzed after a near-total dissolution for As, Cd, Cr, Cu, Ni, Pb, and Zn by Direct-Current Plasma spectroscopy. Data and complementary information were freely available online at http://sigminesfrance.brgm.fr (now moved to http://infoterre.brgm.fr, Géochimie-Inventaire Minier). Less than 15 % of the dataset display values below the limit of detection (LOD) for all elements except Cd, which display 93% data below LOD (Table 1). These data contain useful information and were thus substituted by a value of 75 % of the limit of detection as recommended in environmental sciences [START_REF] Baxter | Sample size and related issues in the analysis of lead isotope data[END_REF][START_REF] Helsel | Statistics for censored environmental data using Minitab and R, Wiley series in statistics in practice[END_REF].

A focus on the three Pb-Ag mining districts was realized at the same spatial scale (1:40 000, ~36 km²) from the BRGM dataset: 64 samples are available in Brandes, 50 samples in Peisey-Nancroix, and 74 samples in Macôt-La Plagne. At this scale, more samples would be required in order to provide better statistical results.

Data processing and statistical treatment

The Quantum GIS free software was used for mapping (QGIS Development Team, 2017), and the rgr, akima, segmented, and StatDA packages of the R free software (version 3.3.2--Sincere Pumpkin Patch) were used for statistical treatment (R Core Team, 2017).

Single-element approaches

Exploratory data analysis

Exploratory Data Analysis (EDA) is a robust method for the study of patterns and the identification of outliers [START_REF] Reimann | Statistical data analysis explained: applied environmental statistics with R[END_REF][START_REF] Carranza | Analysis and mapping of geochemical anomalies using logratiotransformed stream sediment data with censored values[END_REF], using descriptive statistics and graphical tools [START_REF] Tukey | Exploratory data analysis, Addison-Wesley series in behavioral science: quantitative methods[END_REF]. EDA-based mapping allows subdivision of the dataset into seven categories: (i) far low outliers, (ii) near low outliers, (iii) low background, (iv) background, (v) high background, (vi) near high outliers, and (vii) and far high outliers (see detailed methodology in [START_REF] Camizuli | Exploratory Data Analysis (EDA)[END_REF].

As geochemical analyses may vary according to the lithological units, a transformation procedure was applied to the dataset before computing threshold by EDA (see Pb in Fig. 2a), using the following equation [START_REF] Bonham-Carter | Background corrections to stream geochemical data using digitized drainage and geological maps: application to selwyn basin, yukon and northwest territories[END_REF][START_REF] Carranza | A catchment basin approach to the analysis of reconnaissance geochemical-geological data from Albay Province, Philippines[END_REF][START_REF] Yusta | Anomaly threshold estimation and data normalization using EDA statistics: application to lithogeochemical exploration in Lower Cretaceous Zn-Pb carbonate-hosted deposits, Northern Spain[END_REF]: Zi,j = (Xi,j-medi,j)/MADi,j

(1) where Zi,j represents the standardized data value for the element i of the lithology j. The MAD -Median Absolute Deviation -is equivalent to the standard deviation [START_REF] Tukey | Exploratory data analysis, Addison-Wesley series in behavioral science: quantitative methods[END_REF], and is defined as MADi,j = med [|Xi,j-medi,j|], where med represents the median. In the study region, the complex geological units have been simplified into eight lithological categories: (i) granitoids ( n=635), (ii) superficial deposits (n=160), (iii) detrital rocks (n=239), (iv) carbonates (n=2700), (v) gneisses ( n=286), (vi) schists (n=2267), (vii) mafic rocks (n=351), and (viii) migmatites ( n=254).

With such a transformation, however, half the values are negative and cannot be logtransformed, as recommended by [START_REF] Reimann | Statistical data analysis explained: applied environmental statistics with R[END_REF] for strongly skewed data. Therefore, the following transformation was performed to make all values greater or equal to 1 (see Pb in Fig. 2b): zi,j = 1 + Zi,j -min(Zi,j)

(2) Once the data became comparable, six Tukey thresholds were computed, for each considered element, defining the seven anomalies categories mentioned above for EDA mapping. The thresholds were calculated using the map.eda7 function of the rgr package [START_REF] Garrett | The 'rgr' package for the R Open Source statistical computing and graphics environment -a tool to support geochemical data interpretation[END_REF][START_REF] Garrett | rgr: Applied Geochemistry EDA[END_REF], supplementing the R software (R Core Team, 2017).

Fractal-based model

The fractal approach consists in an improvement with respect to the empirically-driven EDAmethod, as it integrates the spatial correlation of geochemical data and the scaleindependent features of geochemical landscapes (i.e., objects presenting a fractal dimension; [START_REF] Carranza | Geochemical anomaly and mineral prospectivity mapping in GIS[END_REF]. Fractal models are based on the geochemical distribution patterns resulting from numerous geological processes that took place at different rates and different temporal and geographical scales [START_REF] Zuo | Fractal/multifractal modeling of geochemical data: A review[END_REF].

The concentration-area plot, developed by [START_REF] Cheng | The separation of geochemical anomalies from background by fractal methods[END_REF], assesses whether the data are spatially multi-fractal, which is useful to determine the occurrence of multiple populations that are spatially dependent in the dataset. It can be used to determine the practical limits, upper or lower bounds, of the influence of the biogeochemical processes behind the spatial distribution of the data.

The caplot function of the rgr package [START_REF] Garrett | The 'rgr' package for the R Open Source statistical computing and graphics environment -a tool to support geochemical data interpretation[END_REF][START_REF] Garrett | rgr: Applied Geochemistry EDA[END_REF], supplementing the R software (R Core Team, 2017) was used to compute the concentration-area plot (see detailed methodology in [START_REF] Monna | Tracking archaeological and historical mines using mineral prospectivity mapping[END_REF]. The breakpoints (thresholds) of the straight lines, fitted by ordinary least square regression, were estimated through the segmented function of the segmented package. This function allows to express regression model with segmented relationships. The slope, intercept and adjusted-R-squared have been computed for each subset.

Multi-element approach: sediment quality index

The Sediment Quality Index (SQI) developed by the Canadian Council of Ministers of the Environment (CCME, 2002) and modified later by [START_REF] Marvin | Application of a Sediment Quality Index to the Lower Laurentian Great Lakes[END_REF], was chosen to describe the overall quality of the streambed sediments. Originally, the guideline values, used by the Canadian government, were the Probable Effect Levels (PELs, [START_REF] Marvin | Application of a Sediment Quality Index to the Lower Laurentian Great Lakes[END_REF]. However, [START_REF] Macdonald | Development and Evaluation of Consensus-Based Sediment Quality Guidelines for Freshwater Ecosystems[END_REF] have developed consensus-based Sediment Quality Guidelines (SQGs) for several contaminants, including a Threshold Effect Concentration (TEC) and a Probable Effect Concentration (PEC). For this study, as there is a lack of clear directives in France, and even in the European Union [START_REF] Besten | Biological effects-based sediment quality in ecological risk assessment for European waters[END_REF], the PEC values of As, Cd, Cr, Cu, Ni, Pb, and Zn (Table 1) were selected as reference values to calculate the SQIs. These values correspond to the levels above which harmful effects on sediment-dwelling organisms are expected to occur frequently [START_REF] Macdonald | Development and Evaluation of Consensus-Based Sediment Quality Guidelines for Freshwater Ecosystems[END_REF]. An index score was computed for each sample, taking into account the number of variables exceeding the guideline values, as well as the magnitude of their respective excess. Calculation of SQIs is detailed in Camizuli et al. (2014b). Final SQIs are expressed as a percentage, thus permitting inter-sample comparisons. To facilitate data representation, the five SQI categories of sediment quality defined by [START_REF] Marvin | Application of a Sediment Quality Index to the Lower Laurentian Great Lakes[END_REF] have been used: poor quality (0-45), marginal quality (45-60), fair quality (60-80), good quality (80-95), and excellent quality (95-100).

Results

Geochemical data

Maximum chemical element concentrations in streambed sediments vary widely and can reach extremely high values: 1076 mg.kg -1 for As, 45 mg.kg -1 for Cd, 1700 mg.kg -1 for Cr, 1530 mg.kg -1 for Cu, 1571 mg.kg -1 for Ni, 8900 mg.kg -1 for Pb, and 8190 mg.kg -1 for Zn. Except for Cd, which displays 93% data below LOD, the other elements exhibit asymmetric distributions, which tend to be more compatible with normal law once log10-transformed, as shown by changes in skewness (Table 1).

Single-element approaches 4.2.1. Exploratory Data Analysis

The percentage of data per EDA category is comparable for each element, with 49-50 % of the sample (n ~ 3400 samples) in the background category (Table 2). The distribution in the high outlier categories (near and far) vary according to the elements. Most geochemical anomalies are related to arsenic, then to lead and zinc, and finally to copper, which presents the smallest anomaly number.

Anomalies can be very localized; their spatial distribution is different for each element (Fig. 3a, Fig. S1a). 5.2% of streambed sediments present a high As anomaly. Most As anomalies are geologically controlled, as evidenced (i) from a cluster of more than 200 samples in the high outliers categories localized in the Série Satinée micaschists around Albertville, and (ii) from the alignment of anomalies along the major tectonic structures on the SE side of the External Crystalline Massifs (contact of the Variscan basement and the Mesozoic sedimentary cover). Some other anomalies, which are less spread, are associated with former mining activities (e.g. Chalanches mine) or to the occurrence of unmined hydrothermal veins.

Lead and zinc anomalies display a similar spatial pattern, which is more heterogeneous than for arsenic, and well distributed all over the study area, suggesting no lithological control. 5.1% of streambed sediments present a high Pb anomaly. Opposite to As, Pb anomalies are mainly related to former mining activity (e.g., Macôt-La Plagne and Peisey-Nancroix district; Brandes district; Rocheray district near St-Jean de Maurienne; Chalanches district, etc.). However, a good correlation is observed between Pb anomalies and major tectonic structures or unmined hydrothermal veins. Nickel anomalies are less numerous (3%) and correspond to several limited clusters that are directly related to ultramafic or mafic rocks in the Belledonne, Grandes Rousses and Vanoise massifs.

Copper anomalies (2.3%) are mostly dispersed within the External Crystalline Massifs (Grandes-Rousses, Belledonne, Beaufortain). Even though copper mines and deposits are located in this domain (e.g., St-Georges d'Hurtières), there is no clear relationship between Cu anomalies and mining.

Fractal-based model

With the fractal approach, the number of samples presenting a high or very high anomaly ranges between 0.5 % (n=33) for Ni and 2 % (n=143) for As (Table 3). This method discriminates few clusters, but the spatial distribution and the intensity of the anomalies are different for each element (Fig. 3b, Fig. S1b).

The distribution of As anomalies is similar to that of EDA ones: the lithological control of the Série Satinée micaschists remains evident (Albertville area; ~ 100 samples out of 143), but clusters are highlighted along the internal contact of the ECM, which are related to the occurrence of Au veins in the Glandon and Madeleine passes.

Lead and zinc anomalies display similar pattern; they are mainly localized in the vicinity of former Pb-Zn-Ag mines. Nickel anomalies are very rare (0.5% of samples) and localized downstream mafic rocks, providing evidence of a strong lithological control. Copper anomalies are scattered in the ECM displaying a pattern similar to EDA one, controlled essentially by geology but not by former mining activity.

Multi-element approach

More than 60 % of the 6892 samples belong to the good or excellent SQI categories (Table 4). This percentage reaches ~95 % if the fair SQI category is included. Several clusters with low SQIs (5 % of the samples) can be distinguished in Fig. 4. The largest cluster is located in around Albertville, in the Série Satinée micaschists of Beaufortain and Lauzière. Another striking cluster follows the Arc river on approximately 7 km near Bonneval-sur-Arc, in the extreme east of the study area, were mafic and ultramafic rocks outcrop. The other anomalies are mainly located in the Belledonne, Grandes Rousses, and Vanoise massifs. The SQI pattern doesn't provide clear evidence of a control by former mines.

Local anomalies

Figure 5 shows Pb anomalies obtained from EDA and fractals methods and SQI maps for Brandes, Peisey-Nancroix and Macôt-La Plagne Pb-Ag districts that were mined at different historical periods. Different numbers of anomalies are obtained from each method. The EDA method is the most sensitive to Pb anomalies, despite normalization to the geological background, with 22 anomalies in Brandes (Pb contents ranging from 147 mg.kg -1 to 8900 mg.kg -1 ), 5 anomalies in Peisey-Nancroix (from 139 mg.kg -1 to 583 mg.kg -1 ) and 14 anomalies in Macôt-La Plagne (from 103 mg.kg -1 to 3200 mg.kg -1 ). With the fractal approach, 9 anomalies are highlighted in Brandes, only 1 in Peisey-Nancroix and 8 in Macôt-La Plagne.

In Peisey-Nancroix, streambed sediments still have a locally high Pb content resulting in localized Pb anomalies with the EDA method, despite mine closure more than 120 years ago. However, these anomalies are not evidenced using the fractal method, and the SQIs indicate at least a fair quality of the sediments. Therefore, the spatial impact of mining activity seems very limited nowadays. In Macôt-La Plagne, higher extracted volumes and more recent mine closure are responsible for higher anomalies evidenced from the three methods: 5 samples present a poor or a marginal quality. Even though Brandes is the oldest mining district, it presents the largest Pb anomalies and the most samples with poor or marginal quality (n=9).

Discussion

Comparison of the statistical methods

The anomaly maps are slightly different according to the method used, as each method has its own definition of the geochemical background [START_REF] Reimann | Statistical data analysis explained: applied environmental statistics with R[END_REF]. The two single-element methods investigated provide similar results concerning the clusters with the highest anomalies, but the fractal-based model discriminates more the dataset than the EDA approach. The multi-element method (SQI) combines all single-element clusters into a single map. However, if an element is too dominant, the SQI value is driven by this element, as for arsenic in the Série Satinée micaschists, near Albertville (Fig. 3). Areas displaying positive anomalies with single-element approaches are not always characterized by poor SQI, thus suggesting limited dispersal of trace metals for these areas. In Peisey-Nancroix, this can be related to the rather small size of the exploitation compared to the two other sites.

At local scale, similar results are observed: (i) the fractal method is more discriminant than the EDA approach, and (ii) SQIs are generally driven by a single element (Fig. 5). In Brandes, poor and marginal sediments are clearly related to the Pb anomalies evidenced from the EDA and fractal methods. In Macôt-La Plagne, the poor and marginal samples are either related to Pb or to Ni anomalies. In Peisey-Nancroix, all sediments have a fair to good quality, despite high contents in some trace metals.

To study the spatial distribution of TMs in sediments, the choice of the statistical method must be based on the purpose of the study [START_REF] Yousefi | Introduction to the special issue: GIS-based mineral potential targeting[END_REF]. Some elements better represent the rock types, such as Fe, whereas others are better appropriated to study the occurrence of mineral deposits (e.g., Pb, Zn; [START_REF] Monna | Tracking archaeological and historical mines using mineral prospectivity mapping[END_REF][START_REF] Kirkwood | Stream sediment geochemistry as a tool for enhancing geological understanding: An overview of new data from south west England[END_REF]. The EDA method is well suited to track archaeological sites and plan field prospection [START_REF] Camizuli | Exploratory Data Analysis (EDA)[END_REF], as it focuses on a single element and as more anomalies are evidenced. This makes it easier to highlight contamination gradients downstream a hot spot, which could be related to an ancient mining site. The fractal method can be used for screening the main clusters, as only the highest anomalies are highlighted. The concentration-area model was used for this paper, but other fractal models can be successfully used to delineate geochemical anomalies (Afzal et al., 2016;[START_REF] Parsa | Multifractal interpolation and spectrum-area fractal modeling of stream sediment geochemical data: Implications for mapping exploration targets[END_REF]. The SQI introduces the probable effect on the aquatic biota and therefore can be used complementary to other environmental analyses to help determining the impact of mining activity on the ecosystem (Camizuli et al., 2014b) or more generally for river plan management [START_REF] Gerber | Application of a Sediment Quality Index for the assessment and monitoring of metals and organochlorines in a premier conservation area[END_REF]Ahn et al., 2019).

Influencing parameters for geochemical anomalies at regional scale

Geochemical anomalies are either driven by natural or by anthropogenic parameters. Some specific lithologies have high background contents in one or several chemical elements, thus highlighting natural anomalies. For example, the largest As cluster is clearly related to the As-rich Série Satinée micaschists in this region of Albertville, and the Ni clusters are related to outcropping mafic rocks. As major tectonic contacts localize intense fluid flow, they are favorable zones for ore deposition (e.g., the contact between the Variscan crystalline basement and the Mesozoic sedimentary cover). If some of these mineralizations have been mined, many veins that were too small to be mined are contributing to local anomalies and trace metals dispersal through erosion. Mining activities enhance spreading of TMs in the environment by modifying the landscape and the hydrological network [START_REF] Stöllner | Mining and economy: a discussion of spatial organisations and structures of early raw material exploitation[END_REF][START_REF] Bertrand | Active channel width as a proxy of sediment supply from mining sites in New Caledonia[END_REF]. Downstream the largest mines and mining districts, Pb anomalies are systematics (Fig. 3). Comparison with a more precise mining database (unpublished data, R. Durand) indicates that most Pb anomalies are related to mining activity. In areas where the anomalies are disconnected with any known ore deposits or specific lithology, field prospection is required in order to determine their origin.

At regional scale, this study showed that TMs spatial distribution is either due to metal-rich lithologies, to the occurrence of mineralized veins or to former mining activity. The size and the shape of the clusters can be linked to the origin of the anomalies: mining activity and the occurrence of mineralized veins result in localized clusters, whereas metal-rich lithologies result in large clusters.

The sampling mesh size must also be taken into consideration together with the hydrographic network, the topography [START_REF] Shahrestani | Dilution correction equation revisited: The impact of stream slope, relief ratio and area size of basin on geochemical anomalies[END_REF] and the climate [START_REF] Montes-Avila | Chemical mobility of inorganic elements in stream sediments of a semiarid zone impacted by ancient mine residues[END_REF] to study the spatial distribution of geochemical data. Moreover, the sediment grain size fraction is also an important parameter as metal content increases with decreasing grain size [START_REF] Mackenzie | Investigation of contaminant metal dispersal from a disused mine site at Tyndrum, Scotland, using concentration gradients and stable Pb isotope ratios[END_REF]. In this study, the sampling density is ~1 sample per km², whereas in the Morvan massif, the same database shows a sampling density of 2-3 samples per km² [START_REF] Lambert | Les données géochimiques et alluvionnaires de l'Inventaire minier du territoire national. Constitution d'une base de données exhaustive[END_REF][START_REF] Monna | Tracking archaeological and historical mines using mineral prospectivity mapping[END_REF], providing a better resolution to investigate geochemical anomalies at the district scale than in the Alps. This difference in the sampling mesh is due to the larger scale and the topography of the current study area. Nonetheless, in other studies, the sampling densities are closer to 1 sample per 10 or 20 km² [START_REF] Lapworth | Geochemical mapping using stream sediments in west-central Nigeria: Implications for environmental studies and mineral exploration in West Africa[END_REF][START_REF] Zuluaga | Stream sediment geochemical mapping of the Mount Pinatubo-Dizon Mine area, the Philippines: Implications for mineral exploration and environmental risk[END_REF], which is not dense enough to allow focusing on a district scale. Such low sampling densities help to the recognition of mineral systems during the early stages of exploration [START_REF] Reimann | Low density geochemical mapping and mineral exploration: application of the mineral system concept[END_REF].

A sampling mesh of 1 sample per km 2 is thus a minimum for detailed studies at the district scale.

1000 years of trace metals dispersal from mining activity

Focusing on a single mining district provides information concerning the ore type, the extent of the exploitation and the extracted volume. This study focused on three Pb-Ag districts that were mined at different historical periods: the Middle Ages for Brandes, the 17th-19th c. in Peisey-Nancroix, and Modern Times in Macôt-La Plagne. At this local scale, with the EDA approach it is easier to follow the TMs path in the streambed sediments using the EDA approach. Indeed, as shown in this study, lead concentration decreases with increasing distance downstream mining areas [START_REF] Mackenzie | Investigation of contaminant metal dispersal from a disused mine site at Tyndrum, Scotland, using concentration gradients and stable Pb isotope ratios[END_REF]Camizuli et al., 2014b). As the three mining districts are situated at 1500-2000 m a.s.l. in a mountainous region, torrential streaming contributes to trace metal dispersal. The largest anomaly clusters are observed in Brandes and in Macôt-La Plagne, which were the most important mining sites during the Middle Ages and Modern Times respectively, resulting in the extraction of several hundred thousand tons of ore for each site [START_REF] Meloux | Alpes-Nord. Bilan synthétique de dix années de prospection[END_REF]pers. com. M.Ch. Bailly-Maître). Peisey-Nancroix, which was a smaller exploitation, shows much more localized number of anomalies in figure 5. Multi-element approach for this site does not allow highlighting the impacts on the aquatic ecosystem with the current dataset. Comparison of these three mining sites suggest that, at local scale, the anomalies are related to the volume of extracted ore rather than to the timing and duration of mining. This supports the idea that old mining sites still release trace metals over centuries and up to a thousand years (Camizuli et al., 2018).

As shown from SQIs, 95 % sediments analyzed in the northern French Alps have a fair to good quality, suggesting only 5% sediments could present a threat to biota. This percentage increases with more recent mining activities [START_REF] Lecce | Floodplain storage of sediment contaminated by mercury and copper from historic gold mining at Gold Hill, North Carolina, USA[END_REF]. In areas showing marginal to poor sediment quality, further investigations on the bioavailable fraction are needed to evaluate this possible effect on biota [START_REF] Casiot | Hydrological and geochemical control of metals and arsenic in a Mediterranean river contaminated by acid mine drainage (the Amous River, France); preliminary assessment of impacts on fish (Leuciscus cephalus)[END_REF]Camizuli et al., 2014b;[START_REF] Garcia-Ordiales | Stream bottom sediments as a means to assess metal contamination in the historic mining district of Almadén (Spain)[END_REF][START_REF] Gutiérrez | Abandoned Pb-Zn mining wastes and their mobility as proxy to toxicity: A review[END_REF] and on human health [START_REF] Varrica | Lead isotopic fingerprint in human scalp hair: The case study of Iglesias mining district (Sardinia, Italy)[END_REF].

Concluding remarks

Three statistical methods (EDA, fractals and SQI) derived from modern mining exploration and ecological management plan were applied on geochemical data from streambed sediments of a former mining territory (the northern French Alps). These methods were used in order to highlight geochemical anomalies and to investigate the long-term dispersal of trace metals and their possible impacts on the environment up to a thousand years after mine closure. This study investigated two different scales: the regional scale (the northern French Alps) and the district scale (focusing on the three most important Pb-Ag mining sites -Brandes, Peisey-Nancroix and Macôt-La Plagne). Geochemical anomalies appear to be driven either by natural parameters, such as a metal-rich lithologies or the occurrence of hydrothermal veins related to major tectonic contacts, or by anthropogenic concentration related to former mining activity (extraction and ore processing). Some geochemical anomalies may be related to unknown ore deposits that may have a future economic potential or may have an archeological interest if complementary field surveys provide evidence of ancient mining activity.

In former mining districts, despite mine closure happened tens to a thousand years ago, local geochemical anomalies are clearly related to the remobilization of trace metals from the extraction and ore processing sites, providing evidence of very long-term remanence and dispersal of trace metals. However, the SQIs provide evidence of a very limited residual environmental impact.

Studying trace metal dispersal in ancient mining territories is thus a good tool to understand and predict long-term environmental impacts of future mining activity. Multiplying such studies in various geological and climatic contexts would contribute to improving knowledge on the parameters that control long-term trace metal dispersal. This is necessary to develop strategies to prevent them, and thus to minimize the long-term environmental impacts of mining activity. Only the high anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1. [START_REF] Kabata-Pendias | Trace elements in soils and plants[END_REF][START_REF] Macdonald | Development and Evaluation of Consensus-Based Sediment Quality Guidelines for Freshwater Ecosystems[END_REF].

Tables and table captions

Figure captions Figure 1 .

 1 Figure captions Figure 1. Overview of the study area, showing the location of the three targeted sites (Br for Brandes, MP for Macôt-La Plagne, PN for Peisey-Nancroix) and other mines cited in the text (Ch for Chalanches, Ro for Rocheray, SGH for St-Georges d'Hurtières). (a) Simplified lithological map of the northern French Alps modified from the 1:1 000 000 BRGM geological map and OneGeology data, (b) mining sites and ore deposits from the BRGM database (topography from BD ALTI® 25 m; hydrographic network from BD Carthage®), and (c) BRGM sampling location of streambed sediments. The contoured area in (a) and (b) delimits the BRGM streambed sediment database.Figure 2. Example of data standardization according to the lithological units for Pb. (a) Data without standardization, (b) data after standardization, UIF: upper inner fence, UOF: upper outer fence. Sup. deposits stands for superficial deposits.Figure 3. Anomaly maps for As, Cu, and Pb resulting from EDA (a) and fractal (b) approaches. Only the highest anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.Figure 4. Sediment Quality Index maps highlighting the poor and marginal anomalies, and showing location of the mines. See legend of the lithological map in Fig. 1.Figure 5. Comparison of Pb anomaly distributions from the EDA and fractal methods, and SQI maps for the three studied mining districts: (a) the medieval site of Brandes, (b) the historical Peisey-Nancroix mine, and (c) the modern Macôt-La Plagne mine. Supplementary data Figure S1. Anomaly maps for Ni and Zn resulting from EDA (a) and fractal (b) approaches. Only the high anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.

Figure 2 .

 2 Figure captions Figure 1. Overview of the study area, showing the location of the three targeted sites (Br for Brandes, MP for Macôt-La Plagne, PN for Peisey-Nancroix) and other mines cited in the text (Ch for Chalanches, Ro for Rocheray, SGH for St-Georges d'Hurtières). (a) Simplified lithological map of the northern French Alps modified from the 1:1 000 000 BRGM geological map and OneGeology data, (b) mining sites and ore deposits from the BRGM database (topography from BD ALTI® 25 m; hydrographic network from BD Carthage®), and (c) BRGM sampling location of streambed sediments. The contoured area in (a) and (b) delimits the BRGM streambed sediment database.Figure 2. Example of data standardization according to the lithological units for Pb. (a) Data without standardization, (b) data after standardization, UIF: upper inner fence, UOF: upper outer fence. Sup. deposits stands for superficial deposits.Figure 3. Anomaly maps for As, Cu, and Pb resulting from EDA (a) and fractal (b) approaches. Only the highest anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.Figure 4. Sediment Quality Index maps highlighting the poor and marginal anomalies, and showing location of the mines. See legend of the lithological map in Fig. 1.Figure 5. Comparison of Pb anomaly distributions from the EDA and fractal methods, and SQI maps for the three studied mining districts: (a) the medieval site of Brandes, (b) the historical Peisey-Nancroix mine, and (c) the modern Macôt-La Plagne mine. Supplementary data Figure S1. Anomaly maps for Ni and Zn resulting from EDA (a) and fractal (b) approaches. Only the high anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.

Figure 3 .

 3 Figure captions Figure 1. Overview of the study area, showing the location of the three targeted sites (Br for Brandes, MP for Macôt-La Plagne, PN for Peisey-Nancroix) and other mines cited in the text (Ch for Chalanches, Ro for Rocheray, SGH for St-Georges d'Hurtières). (a) Simplified lithological map of the northern French Alps modified from the 1:1 000 000 BRGM geological map and OneGeology data, (b) mining sites and ore deposits from the BRGM database (topography from BD ALTI® 25 m; hydrographic network from BD Carthage®), and (c) BRGM sampling location of streambed sediments. The contoured area in (a) and (b) delimits the BRGM streambed sediment database.Figure 2. Example of data standardization according to the lithological units for Pb. (a) Data without standardization, (b) data after standardization, UIF: upper inner fence, UOF: upper outer fence. Sup. deposits stands for superficial deposits.Figure 3. Anomaly maps for As, Cu, and Pb resulting from EDA (a) and fractal (b) approaches. Only the highest anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.Figure 4. Sediment Quality Index maps highlighting the poor and marginal anomalies, and showing location of the mines. See legend of the lithological map in Fig. 1.Figure 5. Comparison of Pb anomaly distributions from the EDA and fractal methods, and SQI maps for the three studied mining districts: (a) the medieval site of Brandes, (b) the historical Peisey-Nancroix mine, and (c) the modern Macôt-La Plagne mine. Supplementary data Figure S1. Anomaly maps for Ni and Zn resulting from EDA (a) and fractal (b) approaches. Only the high anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.

Figure 4 .

 4 Figure captions Figure 1. Overview of the study area, showing the location of the three targeted sites (Br for Brandes, MP for Macôt-La Plagne, PN for Peisey-Nancroix) and other mines cited in the text (Ch for Chalanches, Ro for Rocheray, SGH for St-Georges d'Hurtières). (a) Simplified lithological map of the northern French Alps modified from the 1:1 000 000 BRGM geological map and OneGeology data, (b) mining sites and ore deposits from the BRGM database (topography from BD ALTI® 25 m; hydrographic network from BD Carthage®), and (c) BRGM sampling location of streambed sediments. The contoured area in (a) and (b) delimits the BRGM streambed sediment database.Figure 2. Example of data standardization according to the lithological units for Pb. (a) Data without standardization, (b) data after standardization, UIF: upper inner fence, UOF: upper outer fence. Sup. deposits stands for superficial deposits.Figure 3. Anomaly maps for As, Cu, and Pb resulting from EDA (a) and fractal (b) approaches. Only the highest anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.Figure 4. Sediment Quality Index maps highlighting the poor and marginal anomalies, and showing location of the mines. See legend of the lithological map in Fig. 1.Figure 5. Comparison of Pb anomaly distributions from the EDA and fractal methods, and SQI maps for the three studied mining districts: (a) the medieval site of Brandes, (b) the historical Peisey-Nancroix mine, and (c) the modern Macôt-La Plagne mine. Supplementary data Figure S1. Anomaly maps for Ni and Zn resulting from EDA (a) and fractal (b) approaches. Only the high anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.

Figure 5 .

 5 Figure captions Figure 1. Overview of the study area, showing the location of the three targeted sites (Br for Brandes, MP for Macôt-La Plagne, PN for Peisey-Nancroix) and other mines cited in the text (Ch for Chalanches, Ro for Rocheray, SGH for St-Georges d'Hurtières). (a) Simplified lithological map of the northern French Alps modified from the 1:1 000 000 BRGM geological map and OneGeology data, (b) mining sites and ore deposits from the BRGM database (topography from BD ALTI® 25 m; hydrographic network from BD Carthage®), and (c) BRGM sampling location of streambed sediments. The contoured area in (a) and (b) delimits the BRGM streambed sediment database.Figure 2. Example of data standardization according to the lithological units for Pb. (a) Data without standardization, (b) data after standardization, UIF: upper inner fence, UOF: upper outer fence. Sup. deposits stands for superficial deposits.Figure 3. Anomaly maps for As, Cu, and Pb resulting from EDA (a) and fractal (b) approaches. Only the highest anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.Figure 4. Sediment Quality Index maps highlighting the poor and marginal anomalies, and showing location of the mines. See legend of the lithological map in Fig. 1.Figure 5. Comparison of Pb anomaly distributions from the EDA and fractal methods, and SQI maps for the three studied mining districts: (a) the medieval site of Brandes, (b) the historical Peisey-Nancroix mine, and (c) the modern Macôt-La Plagne mine. Supplementary data Figure S1. Anomaly maps for Ni and Zn resulting from EDA (a) and fractal (b) approaches. Only the high anomaly categories were reported to facilitate reading. See legend of the lithological map in Fig. 1.

Table 1 .

 1

		Fe2O3	MnO	As	Cd	Cr	Cu	Ni	Pb	Zn
		%	%	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1
	Limit of detection	0.1	0.01	20	2	10	10	10	10	10
	Min.	<0.1	<0.01	<20	<2	<10	<10	<10	<10	<10
	Q1	4.4	0.09	29	<2	45	22	28	32	86
	Median (Q2)	5.7	0.13	42	<2	60	32	40	43	112
	Mean	5.8	0.16	59	<2	71	39	49	66	138
	Q3	7	0.18	62	<2	80	45	55	60	149
	Max.	25	2.7	1076	45	1700	1530	1571	8900	8190
	Skewness	0.71 (-1.87)	5.07 (-0.11)	5.66 (0.72)	nc	8.98 (0.07)	14.33 (0.05)	11.92 (0.07)	23.93 (1.58)	23.99 (-0.64)
	nb. cen	3	21	942	6411	56	397	215	45	93
	(\%)	(0.04)	(0.3)	(13.7)	(93)	(0.8)	(5.8)	(3)	(0.7)	(1.3)
	European
	topsoils			11.6	0.28	94.8	17.3	37	32	68.1
	(mean) 1
	TEC 2			9.79	0.99	43.4	31.6	22.7	35.8	121
	PEC 3			33	4.98	111	149	48.6	128	459
	1 Kabata-Pendias (2011)
	2 MacDonald et al. (2000)
	3 MacDonald et al. (2000)

Table 1 .

 1 Descriptive statistics for Fe2O3, MnO, As, Cd, Cr, Cu, Ni, Pb and Zn concentrations on the BRGM Alps dataset. Skewness after log10-transformation is indicated between parentheses, nb. cen are the number of data below the limit of detection and between parenthesis their corresponding percentages. nc stands for not calculated, TEC for Threshold Effect Concentration and PEC for Probable Effect Concentration. Data from http://sigminesfrance.brgm.fr. Concentrations found in the literature are also indicated

Table 3 .

 3 Summary of the parameters for the concentration-area plot of each study element (fractal approach).

	Table 3.
								nb. of samples
	As	range	slope	intercept r-squared	interpretation	Akima BRGM dataset	%
		$<$ 25	-1.163	116.79	0.998	Very low background	602	1400	20.3 %
		25-32 32-41	-1.733 -2.327	131.62 150.53	0.998 1	Low background Low background	590 1022	776 1151	28.0 %
		41-60 60-82	-1.61 -0.824	119.58 73.43	0.993 0.985	Background Background	1494 839	1738 852	37.6 %
		82-103	-0.136	19.23	0.995	High background	132	309
		103-122 122-213	-0.063 -0.033	11.79 7.81	0.961 0.925	High background High background	57 138	151 284	12.1 %
		213-273	-0.009	3.28	0.879	High background	31	88
	10	273-463	-0.004	1.76	0.859	High anomaly	34	108	1.6 %
	11	$>$ 463	-0.002	1.07	0.968	Very high anomaly	5	35	0.5 %
								nb. of samples
	Cu	range	slope	intercept r-squared	interpretation	Akima BRGM dataset	%
		$<$ 13	-0.774	105.71	0.996	Very low background	212	661	9.6 %
		13-19	-1.227	111.83	0.998	Low background	371	684
		19-24	-2.328	132.44	0.999	Low background	590	791	39.4 %
		24-31	-3.109	151.41	0.999	Low background	1062	1240
		31-43	-2.799	138.99	0.978	Background	1616	1635
		43-66	-0.765	52.84	0.941	Background	820	1230	46.6 %
		66-90	-0.165	15.94	0.955	Background	194	348
		90-126	-0.044	5.35	0.902	High background	65	197	2.9 %
		$>$ 126	-0.001	0.33	0.649	High anomaly	14	106	1.5 %
								nb. of samples
	Ni	range	slope	intercept r-squared	interpretation	Akima BRGM dataset	%
		$<$ 15	-0.615	104.91	0.986	Very low background	227	540	7.8 %
		15-24 24-33	-1.337 -2.027	115.68 132.05	0.995 0.999	Low background Low background	604 890	801 1156	28.4 %
		33-45 45-72	-2.235 -1.083	138.18 84.08	0.999 0.958	Background Background	1326 1363	1672 1901	51.8 %
		72-125	-0.146	20.32	0.946	High background	372	613
		125-256	-0.016	4.97	0.947	High background	100	165	11.4 %
		256-325	-0.003	2.12	0.955	High background	13	11
	10	325-372 372-694	-0.004 -0.002	2.17 1.41	0.957 0.93	High anomaly High anomaly	8 28	0 26	0.4 %
	11	$>$ 694 -3.99E-04	0.5	0.854	Very high anomaly	13	7	0.1 %
								nb. of samples
	Pb	range	slope	intercept r-squared	interpretation	Akima BRGM dataset	%
		$<$ 20	-0.21	102.46	0.826	Very low background	107	356	5.2 %
		20-28 28-37	-1.161 -2.58	122.02 161.26	0.974 0.999	Low background Low background	458 1147	852 1372	32.3 %
		37-55 55-77	-2.153 -0.944	142.9 77.79	0.989 0.969	Background Background	1876 925	2272 1085	48.7 %
		77-162 162-369	-0.092 -0.006	14.39 2.11	0.888 0.914	High background High background	363 56	719 160	12.8 %
		369-890 -2.98E-04	0.35	0.955	High anomaly	7	42	0.6 %
		$>$ 890 -8.94E-05	0.17	0.822	Very high anomaly	5	34	0.5 %
								nb. of samples
	Zn	range	slope	intercept r-squared	interpretation	Akima BRGM dataset	%
		$<$ 21 21-43	-0.043 -0.059	100.2 100.76	0.947 0.89	Very low background Very low background	33 67	132 137	3.9 %
		43-64 64-87	-0.257 -0.755	110.21 142.67	0.928 0.988	Low background Low background	238 855	422 1108	22.2 %
		87-111 111-218	-1.204 -0.463	181.67 91.6	0.998 0.882	Background Background	1383 2130	1584 2867	64.6 %

Table 4 .

 4

	Category	Range	Percentage Nb of samples
	Poor	]-1,45]	1.38	95
	Marginal	]45,60]	3.93	271
	Fair	]60,80]	31.36	2161
	Good	]80,95]	39.61	2730
	Excellent ]95,100]	23.72	1635

Table 4 .

 4 Percentage and number of samples in each SQI category

Table 1 .

 1

		Fe2O3	MnO	As	Cd	Cr	Cu	Ni	Pb	Zn
		%	%	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1	mg.kg -1
	Limit of detection	0.1	0.01	20	2	10	10	10	10	10
	Min.	<0.1	<0.01	<20	<2	<10	<10	<10	<10	<10
	Q1	4.4	0.09	29	<2	45	22	28	32	86
	Median (Q2)	5.7	0.13	42	<2	60	32	40	43	112
	Mean	5.8	0.16	59	<2	71	39	49	66	138
	Q3	7	0.18	62	<2	80	45	55	60	149
	Max.	25	2.7	1076	45	1700	1530	1571	8900	8190
	Skewness	0.71 (-1.87) 5.07 (-0.11) 5.66 (0.72)	nc	8.98 (0.07) 14.33 (0.05) 11.92 (0.07) 23.93 (1.58) 23.99 (-0.64)
	nb. cen (\%)	3 (0.04)	21 (0.3)	942 (13.7)	6411 (93)	56 (0.8)	397 (5.8)	215 (3)	45 (0.7)	93 (1.3)
	European topsoils (mean) 1			11.6	0.28	94.8	17.3	37	32	68.1
	TEC 2			9.79	0.99	43.4	31.6	22.7	35.8	121
	PEC 3			33	4.98	111	149	48.6	128	459
	1 Kabata-Pendias (2011)
	2 MacDonald et al. (2000)
	3 MacDonald et al. (2000)

Table 2 .

 2

	Categories	zAs				zCu		zNi			zPb		zZn
		Cut level	n	%	Cut level	n	%	Cut level	n	%	Cut level	n	%	Cut level	n	%
	Far low outliers	< 0.51	0	0.0	< 0.48	0	0	< 0.64	0	0	< 0.75	0	0	< 0.93	0	0
	Near low outliers	0.51 -1.06	29	0.4	0.48 -1	2	0.03 0.64 -1.25	59	0.9 0.75 -1.41	50	0.7 0.93 -1.66 170 2.5
	Low background	1.06 -2.21 1711 24.8	1 -2.11	1764 25.6 1.25 -2.43 1724 25.0 1.41 -2.65 1674 24.3 1.66 -2.98 22.8
	Background (Q1-Q3)	2.21 -3.6 3439 49.9 2.11 -3.46 3448 50.0 2.43 -3.79 3394 49.2 2.65 -4.03 3459 50.2 2.98 -4.39 49.9
	High background	3.6 -7.49 1354 19.6 3.46 -7.26 1516 22.0 3.79 -7.39 1507 21.9 4.03 -7.56 1361 19.7 4.39 -7.84 20.4
	Near high outliers	7.49 -15.59 298 4.3 7.26 -15.25 139 2.0 7.39 -14.41 164 2.4 7.56 -14.17 231 3.4 7.84 -14.03 236 3.4
	Far high outliers	> 15.59	61	0.9	> 15.25	23	0.3	> 14.41	44	0.6	> 14.17	117 1.7	> 14.03	67	1.0

Table 3 .

 3

								nb. of samples
	As	range	slope	intercept r-squared	interpretation	Akima BRGM dataset %
		$<$ 25	-1.163	116.79	0.998	Very low background	602	1400	20.3
		25-32 32-41	-1.733 -2.327	131.62 150.53	0.998 1	Low background Low background	590 1022	776 1151	28.0
		41-60 60-82	-1.61 -0.824	119.58 73.43	0.993 0.985	Background Background	1494 839	1738 852	37.6
		82-103	-0.136	19.23	0.995	High background	132	309
		103-122 122-213	-0.063 -0.033	11.79 7.81	0.961 0.925	High background High background	57 138	151 284	12.1
		213-273	-0.009	3.28	0.879	High background	31	88
	10	273-463	-0.004	1.76	0.859	High anomaly	34	108	1.6
	11	$>$ 463	-0.002	1.07	0.968	Very high anomaly	5	35	0.5
								nb. of samples
	Cu	range	slope	intercept r-squared	interpretation	Akima BRGM dataset %
		$<$ 13	-0.774	105.71	0.996	Very low background	212	661	9.6
		13-19	-1.227	111.83	0.998	Low background	371	684
		19-24	-2.328	132.44	0.999	Low background	590	791	39.4
		24-31	-3.109	151.41	0.999	Low background	1062	1240
		31-43	-2.799	138.99	0.978	Background	1616	1635
		43-66	-0.765	52.84	0.941	Background	820	1230	46.6
		66-90	-0.165	15.94	0.955	Background	194	348
		90-126	-0.044	5.35	0.902	High background	65	197	2.9
		$>$ 126	-0.001	0.33	0.649	High anomaly	14	106	1.5
								nb. of samples
	Ni	range	slope	intercept r-squared	interpretation	Akima BRGM dataset %
		$<$ 15	-0.615	104.91	0.986	Very low background	227	540	7.8
		15-24 24-33	-1.337 -2.027	115.68 132.05	0.995 0.999	Low background Low background	604 890	801 1156	28.4
		33-45 45-72	-2.235 -1.083	138.18 84.08	0.999 0.958	Background Background	1326 1363	1672 1901	51.8
		72-125	-0.146	20.32	0.946	High background	372	613
		125-256	-0.016	4.97	0.947	High background	100	165	11.4
		256-325	-0.003	2.12	0.955	High background	13	11
	10	325-372 372-694	-0.004 -0.002	2.17 1.41	0.957 0.93	High anomaly High anomaly	8 28	0 26	0.4
	11	$>$ 694 -3.99E-04	0.5	0.854	Very high anomaly	13	7	0.1
								nb. of samples
	Pb	range	slope	intercept r-squared	interpretation	Akima BRGM dataset %
		$<$ 20	-0.21	102.46	0.826	Very low background	107	356	5.2
		20-28 28-37	-1.161 -2.58	122.02 161.26	0.974 0.999	Low background Low background	458 1147	852 1372	32.3
		37-55 55-77	-2.153 -0.944	142.9 77.79	0.989 0.969	Background Background	1876 925	2272 1085	48.7
		77-162 162-369	-0.092 -0.006	14.39 2.11	0.888 0.914	High background High background	363 56	719 160	12.8
		369-890 -2.98E-04	0.35	0.955	High anomaly	7	42	0.6
		$>$ 890 -8.94E-05	0.17	0.822	Very high anomaly	5	34	0.5
								nb. of samples
	Zn	range	slope	intercept r-squared	interpretation	Akima BRGM dataset %
		$<$ 21 21-43	-0.043 -0.059	100.2 100.76	0.947 0.89	Very low background Very low background	33 67	132 137	3.9
		43-64 64-87	-0.257 -0.755	110.21 142.67	0.928 0.988	Low background Low background	238 855	422 1108	22.2
		87-111 111-218	-1.204 -0.463	181.67 91.6	0.998 0.882	Background Background	1383 2130	1584 2867	64.6
		218-294 294-465	-0.032 -0.011	11.56 5.25	0.966 0.893	High background High background	125 89	293 217	7.4
		465-597	-0.003	1.61	0.968	High anomaly	18	53	0.8
	10	$>$ 597 -1.81E-04	0.22	0.975	Very high anomaly	6	79	1.1

Acknowledgements

The authors are grateful to the University of Savoie Mont-Blanc for funding the one-year postdoctoral fellowship for Estelle Camizuli and for funding the TraMines research project, directed by Magali Rossi and Marie Forget (EDYTEM laboratory). This study was also funded by the French national research program EC2CO-Biohefect of the CNRS (TRACES project, grant number 183446). The authors warmly thank Robert Durand for his knowledge concerning local mines and his personal data and particularly Bruno Tourlière, in the BRGM, for providing the sediment database used for this study.

 2. Cut levels obtained from the computation of the log-transformed standardized data (EDA approach). The number of samples in each category is specified. The seven categories are based on the Tukey boxplot [START_REF] Camizuli | Exploratory Data Analysis (EDA)[END_REF][START_REF] Tukey | Exploratory data analysis, Addison-Wesley series in behavioral science: quantitative methods[END_REF].

Chambéry Chambéry

Grenoble