Joachim Tremosa
email: j.tremosa@brgm.fr

Mathieu Debure

Sathya Narayanasamy

Paul-Olivier Redon

Diederik Jacques

Francis Claret

Jean-Charles Robinet

Shale weathering: A lysimeter and modelling study for flow, transport, gas diffusion and reactivity assessment in the critical zone

Keywords: weathering, critical zone, vadose zone, gas diffusion, chemical reactivity, shale

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.

Introduction

The vadose zone have attracted more attention over the past decade as a privileged place for reactions between water, atmosphere, rock and soil. This part of the earth system called the critical zone (CZ) extends from the vegetation to the groundwater [START_REF] Li | Expanding the role of reactive transport models in critical zone processes[END_REF] and is a system that combines chemical, biological, physical, and geological processes supporting life all together [START_REF] Brantley | Crossing disciplines and scales to understand the critical zone[END_REF]. The CZ is subject to hydraulic, geomorphological and climatic events (wind, rainfall, flood, drought) that finally fracture, grind and dissolve rocks [START_REF] Anderson | Physical and chemical controls on the critical zone[END_REF][START_REF] Worthington | Enhancement of bedrock permeability by weathering[END_REF]. The weathering of shales is of particular interest since clay minerals are ubiquitous on the earth's surface and play an important role in chemical element cycles (CO2, O2, nutrients...) in the atmosphere, water and soils. Shales also serve as caprock for waste or subsurface gas storage and for which it is necessary to anticipate shale alteration. The weathering of the rocks occurs over a variety of time and space scales (from mineral grains to watersheds) that makes it difficult to quantify the processes and their rates [START_REF] Littke | Quantification of loss of calcite, pyrite, and organic matter due to weathering of Toarcian black shales and effects on kerogen and bitumen characteristics[END_REF][START_REF] Pacheco | Weathering of plagioclase across variable flow and solute transport regimes[END_REF][START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the Earth sciences[END_REF]. To circumvent such scaling issues, one approach consists of acquiring data in the laboratory and in the field and analysing the results with reactive transport models (RTM) in order to decipher the mechanisms occurring at the earth surface and subsurface [START_REF] Goddéris | Reactive transport models of weathering[END_REF][START_REF] Hadermann | The Grimsel (Switzerland) migration experiment: integrating field experiments, laboratory investigations and modelling[END_REF][START_REF] Mayer | Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage[END_REF][START_REF] Steefel | Reactive transport modeling: An essential tool and a new research approach for the Earth sciences[END_REF].

CZ studies cover a large range of fields including regolith formation that contributes to important processes such as nutrient cycling, carbon sequestration, erosion, and acid rain mitigation [START_REF] Brantley | Crossing disciplines and scales to understand the critical zone[END_REF][START_REF] Jin | Mineral weathering and elemental transport during hillslope evolution at the Susquehanna/Shale Hills Critical Zone Observatory[END_REF]. Reactions with atmospheric O2 and CO2 induce nested alteration fronts in the subsurface from the atmosphere [START_REF] Brantley | Probing deep weathering in the Shale Hills Critical Zone Observatory, Pennsylvania (USA): the hypothesis of nested chemical reaction fronts in the subsurface[END_REF][START_REF] Lerouge | Origin of dissolved gas (CO2, O2, N2, alkanes) in pore waters of a clay formation in the critical zone (Tégulines Clay, France)[END_REF]. In general, the depths of these weathering reactions are unknown, with a close link between the reactivity of the CZ systems and the water table. The water table can act as a buffer zone for contaminants [START_REF] Molins | Coupling between geochemical reactions and multicomponent gas and solute transport in unsaturated media: A reactive transport modeling study[END_REF] as well as reducing the O2 migration rate and preventing mineral oxidization as observed in mining tailings [START_REF] Ouangrawa | Reactive transport modelling of mine tailings columns with capillarity-induced high water saturation for preventing sulfide oxidation[END_REF]. Acid rock drainage is a major concern when shales or mining waste deposits contain reactive sulphide minerals that are unstable in contact with atmospheric O2 [START_REF] Blowes | The pore-water geochemistry and the mineralogy of the vadose zone of sulfide tailings, Waite Amulet, Quebec, Canada[END_REF][START_REF] Littke | Quantification of loss of calcite, pyrite, and organic matter due to weathering of Toarcian black shales and effects on kerogen and bitumen characteristics[END_REF]. In addition to the acidification of the pore water, the release of Fe, SO4 and toxic metals (Cu, Pb, Zn and other trace metals) has also been reported [START_REF] Mayer | Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage[END_REF]. Besides water table fluctuations, carbonate dissolution was one of the mechanisms reported to buffer the system and prevent pore water acidification [START_REF] Brantley | Probing deep weathering in the Shale Hills Critical Zone Observatory, Pennsylvania (USA): the hypothesis of nested chemical reaction fronts in the subsurface[END_REF][START_REF] Ouangrawa | Reactive transport modelling of mine tailings columns with capillarity-induced high water saturation for preventing sulfide oxidation[END_REF]. The prediction of the impact of sulphide-rich tailings on the environment require an understanding of the coupling between hydraulic (flow, transport) and geochemical (mineral dissolution/precipitation, oxido-reduction reactions) processes in the system.

Reactive transport modelling is a powerful tool to investigate biogeochemical systems in the fields of hydrogeology, geochemistry and biogeochemistry [START_REF] Li | Expanding the role of reactive transport models in critical zone processes[END_REF][START_REF] Macquarrie | Reactive transport modeling in fractured rock: A state-of-thescience review[END_REF]. The simulation of metal elements release in the vadose zone, taking into account water-rock-gas interactions, flow and transport processes, is the focus of an increasing number of studies [START_REF] Acero | Multiphase flow and reactive transport model in vadose tailings[END_REF][START_REF] Mayer | Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage[END_REF][START_REF] Mayer | Multicomponent reactive transport modeling in variably saturated porous media using a generalized formulation for kinetically controlled reactions[END_REF][START_REF] Molson | Numerical simulations of pyrite oxidation and acid mine drainage in unsaturated waste rock piles[END_REF]. However, studies have focused on either laboratory experiments (< dm) or the field scale (> 10 m) without considering the intermediate scale that is reached through lysimeter monitoring.

The objective of this study is to provide insight into the chemical evolution of shale tailings containing small amounts of sulphide minerals by considering flow, transport, gas diffusion and water table variations (rainfall, evaporation) at a lysimeter scale. This study may help to reconcile the discrepancies that usually occur between laboratory data and field measurements [START_REF] Guo | Reconciling chemical weathering rates across scales: Application of uranium-series isotope systematics in volcanic weathering clasts from Basse-Terre Island (French Guadeloupe)[END_REF][START_REF] Maher | The mineral dissolution rate conundrum: Insights from reactive transport modeling of U isotopes and pore fluid chemistry in marine sediments[END_REF][START_REF] Pacheco | Role of fractures in weathering of solid rocks: narrowing the gap between laboratory and field weathering rates[END_REF]) by more accurately controlling boundary conditions. The quantification of the oxidation of the tailings material and its associated effects (pH and redox state modifications, in particular) can be used to assess the possible release of sulphates and trace metals into the environment.

Materials and methods

Tailings of Callovian-Oxfordian (COx) shale

The studied site stems from the research conducted in the eastern part of France to build up a deep underground radioactive waste disposal facility in a deep geological formation. The Callovian-Oxfordian shale (e.g. a clay-rich sedimentary formation) layer has been studied as a potential host rock for more than 20 years. The thermo-hydro-mechanical and chemical behaviour of Callovian-Oxfordian (COx) shale was notably investigated in an underground research laboratory located 500 m deep. COx shale is an indurated shale that exhibits a fined grained texture characterised by a porosity between 15 and 20 % [START_REF] De Craen | High sulphate concentrations in squeezed Boom Clay pore water: evidence of oxidation of clay cores[END_REF] and pore diameters mainly ranging between 10 and 100 nm [START_REF] Goddéris | Reactive transport models of weathering[END_REF]. These nanometre-sized pores induce low hydraulic conductivity (10 -13 to 10 -14 m s -1 ) and low diffusion coefficients [START_REF] Descostes | Diffusion of anionic species in Callovo-Oxfordian argillites and Oxfordian limestones (Meuse/Haute-Marne, France)[END_REF][START_REF] Harrington | Gas flow in Callovo-Oxfordian claystone (COx): results from laboratory and field-scale measurements[END_REF]. With a Biot coefficient of 0.6 and a Young modulus of 4 GPa, COx shale is a cohesive and stiff material [START_REF] Charlier | An unsaturated hydromechanical modelling of two in-situ experiments in Callovo-Oxfordian argillite[END_REF]. The mineralogy of COx is constituted by clay minerals (mainly mixed-layered illite-smectite minerals and illite), carbonates (mainly calcite) and quartz. Minor amounts of sulphur minerals (1 to 2 % of pyrite, galena…) and sulphate minerals (celestite) are also found, as well as iron oxide and (oxy-)hydroxides (goethite). Pyrite is known to contain trace amounts of As, Pb, Co, Cu, Ni and Zn [START_REF] Lerouge | Mineralogical and isotopic record of biotic and abiotic diagenesis of the Callovian-Oxfordian clayey formation of Bure (France)[END_REF]. In the conditions of the geological formation, the pristine pore water is in equilibrium with the shale-forming minerals [START_REF] Gaucher | A robust model for pore-water chemistry of clayrock[END_REF] and exhibits a pH of about 7.2. Reductive conditions are found in the shale and pore water, forced by an equilibrium with pyrite.

Several million cubic meters of the Callovo-Oxfordian rocks are expected to be excavated during the building operation of the facility. This material will be stored as tailings at the surface where it will form a technosol [START_REF] Scholtus | Expected evolution of a Technosol derived from excavated Callovo-Oxfordian clay material[END_REF]. Therefore, it is subject to meteorological variations that modify the COx chemical and mechanical properties due to decompaction, disaggregation and leaching.

Mineralogical alterations of various tailings samples were identified by ERM (Etudes Recherches Matériaux, Poitiers, France) company using a JEOL JSM 5600LV scanning electron microscope (SEM) combined with energy dispersive X-Ray spectroscopy (EDS) examinations. It enabled giving a chemical mapping together with the crystallographic information on these tailing samples.

The weathering of the rock stems from its mechanical deconstruction and from pyrite oxidation [START_REF] De Windt | Modeling of oxygen gas diffusion and consumption during the oxic transient in a disposal cell of radioactive waste[END_REF][START_REF] Vinsot | Oxidation front and oxygen transfer in the fractured zone surrounding the Meuse/Haute-Marne URL drifts in the Callovian-Oxfordian argillaceous rock[END_REF] that eventually increase the permeability of the rock from the surface to the depths and water infiltration will increase leaching of chemical elements due to dissolution/precipitation of minerals. A part of the tailings stored at the surface is also expected to be reused as backfill materials in the repository. Understanding and predicting the physico-chemical evolution of the Callovo-Oxfordian tailings is thus required to prevent any potential environmental impacts, to design the tailing management scenarios (geometry, water cycle…) and to evaluate the physico-chemical state of the tailings used as a backfill materials.

Experimental setup

2.2.1.

Lysimeter station

The weathering of the Callovian-Oxfordian shale was studied since 2013 using four stainless steel lysimeters (Figure 1) with a cross-sectional area of 1 m² and a height of 2 m (Andra lysimeter station, Osne-le-Val, France, 48°30'33" N, 5°12'43"E; https://www.ugt-online.de/en/products/lysimetertechnology/). Two lysimeters, called L5 and L6, were filled with COx shale tailings already stored at the surface for 10 years. The other two lysimeters, called L7 and L8, were filled with recently excavated (< 6 months) COx shale from loose tailings. Monoliths were directly cored from aged tailings to fill lysimeters L5 and L6 while the lysimeters L7 and L8 were hand-filled and compacted layer by layer (30 cm) to reach a density of 1,700 kg m -3 , close to the density measured in the aged tailings (1,800 kg m -3 ). Since July 2013, each lysimeter has been continuously weighed by precision scales and is equipped with soil hydrology devices such as tensiometers, water content sensors, temperature sensors and suction probes at depths of 20, 50, 100 and 150 cm. Lysimeters are also equipped with tipping counters for the measurement of basal water leakage.

The combined soil moisture and temperature sensors used (UMP-1, Umwelt Geräte Technik GmbH) are based on frequency domain reflectometry (FDR) method and measure local changes in dielectric permittivity around the probe that can be related to changes in water content. The soil moisture sensors measure the local conductivity around the sensor, hence, the measured water content corresponds both to the water in the large space between the grains of argillaceous rock and in the porosity of the grains. Matrix potential was measured through suction in ceramic filters by tensiometers (Tensio 160, Umwelt Geräte Technik GmbH). At soil saturation, matrix potential is zero and is positive with soil desaturation. The pore size of the ceramic filters is higher than the argillaceous rock pore size and then the measured suction potential corresponds to the one of the intergrain spaces of the tailing material. Water sampling was also performed at different depths using ceramic suction probes (Umwelt Geräte Technik GmbH) and an air-pump operated by the matrix potential value. The sampled water was stored in 500 mL glass bottles where water accumulated between two manual sampling times, every two months. Tipping counters (Umwelt Geräte Technik GmbH) were used for basal drainage monitoring and for leaking water sampling. Each lysimeter was placed on three high-sensitivity weighting cells (10 g resolution). A datalogger has recorded and transmitted the weight of the lysimeters, the matrix potential, the water content, the drainage at the bottom and the temperature at an hourly frequency since January 2014.

The lysimeter station has been equipped since July 2014 with a weather station measuring rainfall, temperature, relative humidity, wind speed, atmospheric pressure and incident solar radiation data.

2.2.2.

Water sample analyses

The water in the lysimeters was sampled every two months from January 2014 to January 2017. The water samples were then filtrated at 0.45 µm and analysed at the Laboratoire Interdisciplinaire des Environnements Continentaux (LIEC, Nancy, France) for total carbon and inorganic carbon contents using a VCSH TOC analyser (Shimadzu Scientific Instruments) and for anions (Cl -, NO2 -, NO3 -2 , PO4 -2 , SO4 -2 , Br -) by ionic chromatography (Dionex). Major and trace elements (Al, As, Ba, Ca, Cd, Co, Cr, Cu, Fe, K, Mg, Mn, Mo, Na, Ni, P, Pb, S, Sb, Se, Zn) were measured at the Laboratoire Sols et Environnement (LSE, Nancy, France) using ICP-OES (iCAP Duo Thermo Scientific).

Additional water samples were also taken from L6 and L7 in March 2016 and November 2017 to determine the 34 S/ 32 S isotopic ratio. Immediately after sampling, the water samples were filtrated at 0.45 µm, and cadmium acetate at 50 g.L -1 was added to trap sulphurs and prevent their oxidation in sulphates. The samples were filtrated at 0.22 µm in the laboratory and sulphates were precipitated with BaCl2. The BaSO4 precipitate was then mixed with vanadium pentoxide (V2O5) and heated at 1000 °C in tin capsules under O2 flux. The gaseous SO2 produced during combustion was purified by gas chromatography (GC) and analysed at the Bureau de Recherches Géologiques et Minières (BRGM, Orléans, France) with an elemental analyser coupled to a continuous flow isotope ratio mass spectrometer (Delta-Plus CF-IRMS, Thermo Fisher Scientific).

Water balance

Water balance calculations allows for estimating the flux of water between the atmosphere, the COx shale in the lysimeters and the drainage at the base of the lysimeter at the lysimeter scale. To reproduce flow in the lysimeter and to determine the hydraulic parameters, it is necessary respecting the water balance to verify the coherence of the recorded data.

In a lysimeter, the water balance can be calculated by considering the precipitation (P), evaporation (EP), drainage at the base of the lysimeter (DR) and the variation in water content (ΔS), according to equation 1 (e.g., [START_REF] Cui | Water Balance and Evapotranspiration Monitoring in Geotechnical and Geoenvironmental Engineering[END_REF][START_REF] Tarantino | Laboratory and field testing of unsaturated soils[END_REF]):

∆ -= -(eq.1) Compared to a water balance usually calculated in soil, the surface run-off is not considered in the lysimeter water balance because the top ledge of the lysimeter impedes such run-off. Rain interception and water uptake by vegetation can also be disregarded because of the absence of vegetation in our lysimeter study. Furthermore, only drainage is considered at the base of the lysimeter, since water can only flow out of the lysimeter at its basis.

The lysimeter and local meteorological data were used to calculate the water balance in each lysimeter on a daily basis from the 01/01/2014 to 09/04/2017. These data include the lysimeter weights, i.e. the water content variation (ΔS), and the volume of drained water measured directly at the lysimeter base (DR).

Precipitation (P) and evaporation (EP) are also involved in the water balance. The local meteorological data measured onsite were completed with data from the Météo France weather stations of Chevillon (station 52123003, located 7 km from the lysimeters) for precipitation and temperature records and Houdelaincourt (station 55248001, located about 20 km away) for evaporation evolution, determined using the Monteith method. Besides some minor differences during rain events between the Osne-le-Val and Chevillon weather stations, the cumulated amount of precipitation was similar between the two stations. Data from Chevillon weather station can therefore reasonably be used to cover the 6-month period before the installation of the weather station at the lysimeter site.

By combining the different available data, various water balances were calculated, to assess the coherence of the data. The first one aimed at calculating the evaporation by water balance on the lysimeters, using the direct measurement of drainage and the variations in lysimeter weight as the variation in water content at the lysimeter scale. Another water balance in the lysimeters was made to calculate the basal drainage by considering the water content variation obtained from the lysimeter mass evolution, the precipitation series (Météo France station for the beginning of 2014 and then from Osne-le-Val weather station) and the evaporation calculated with the Penman-Monteith formula using Osne-le-Val weather station data. A correction factor of 0.2 was applied to the Penman-Monteith evaporation to account for the difference between the potential evaporation and the actual evaporation.

Reactive transport model in unsaturated conditions

A reactive transport model that aims reproducing the data collected in the lysimeters has been set up to evaluate the weathering phenomena occurring in the COx tailings. The approach undertaken here aimed at numerically simulating the water content and the evolution of solute concentrations in the lysimeters during the 2014 to 2017 period, using the meteorological data and the basal drainage in the lysimeters during the same period as boundary conditions. The calculation code HP1 v2.4 [START_REF] Jacques | User manual of the multicomponent variably-saturated flow and transport model HP1. Description, Verification and Examples[END_REF][START_REF] Jacques | The HPx software for multicomponent reactive transport during variably-saturated flow: Recent developments and applications[END_REF]Jacques et al., 2008a;[START_REF] Šimůnek | Multicomponent geochemical transport modeling using HYDRUS-1D and HP1[END_REF], which results from the coupling of Hydrus-1D [START_REF] Simunek | The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media[END_REF] and PHREEQC [START_REF] Parkhurst | Description of Input and Examples for PHREEQC Version 3-a Computer Program for Speciation, Batch-reaction, One-dimensional Transport, and Inverse Geochemical Calculations[END_REF] calculation codes, was used. HP1 can be used to simulate water flow, gas diffusion and solute transport in an unsaturated medium and can consider the set of chemical reactions expected in shale tailings by a thermo-kinetic approach.

While the geometry of the lysimeter was represented using a 1D model, advection and diffusion transport mechanisms occurring in the lysimeter were accounted for by a dual-porosity model.

Because of the mechanical deconstruction of the COx shale during its excavation, the hydraulic properties of the tailings are greatly changed compared to the initial shale and a preferential flow develops between the shale aggregates with little transfers with the water in the porosity of the aggregates. The preferential flow between the soil aggregates was modelled as the mobile porosity zone while the immobile zone corresponded to the porosity inside the COx shale aggregates. The mobile zone corresponded to 10% of the volume of the lysimeter with a porosity of 1 while the immobile zone was considered for the remaining 90% of the volume of the lysimeter, with a porosity of 0.15. The 2 m-high lysimeter was represented as 101 nodes, each 2 cm long, and both a mobile and immobile domain was associated with each of these 101 nodes. A schematic representation of the modelled lysimeter is shown in Figure 2 where the main flow, transport and reactive processes are indicated. Simulation aim to reproduce the data that have been gathered between January 2014 and April 2017.

2.4.1.

Hydraulic and transport model

To simulate water flow and solute transport in the lysimeters, a dual porosity model was used [START_REF] Šimůnek | Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone[END_REF][START_REF] Šimůnek | Modeling Nonequilibrium Flow and Transport Processes Using HYDRUS[END_REF]. This model considered vertical preferential water and gas flow in the zone between the rock aggregates (mobile zone) and a water and solute transfer between the aggregates (immobile zone) and the mobile zone.

Water flow in the mobile zone and water transfer between the mobile and immobile zones in a dual porosity model are described in Hydrus according to the following equations, based on the Richards equation (Gerke and van Genuchten, 1993a;Gerke and van Genuchten, 1993b;[START_REF] Šimůnek | Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone[END_REF]:

= , 1 - (eq.2) = (eq.3) = ! ! (eq.4)
where the indices m and im refer to the mobile and immobile zones, respectively, h is the hydraulic head (L), t and z are the time (T) and vertical space (L) coordinates, respectively, K(h,z) is the relative hydraulic conductivity (L T -1 ), θ is the volumetric water content and Γw is the water transfer between mobile and immobile domains (T -1 ). The water transfer depends on the head difference between the mobile and the immobile domains and is written as = $ ! -! , where αw is the water mass coefficient transfer between the mobile and the immobile zones (T -1 ).

The transfer of chemical components between the two zones of the porous medium is described by the two following mass balance equations (Gerke and van Genuchten, 1993a;[START_REF] Šimůnek | Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone[END_REF]:

% & ' ( % ' ) ! * + = ! ! % , + ! , % ' - -% -! -+ (eq.5) % 1 -) ! * + = -! + (eq.6)
where c is the solute concentration (M L -3 ), Dm is the dispersivity coefficient in the mobile zone (L 2 T - 1 ), s in the sorbed or precipitated concentration of the solute (M M -1 ), fm is the fraction of sorption site or solid in contact with water in the mobile domain, Sm is the source/sink term for geochemical reactions (M L -3 T -1 ) and Γs is the solute transfer between mobile and immobile domains (M L -3 T -1 ).

This solute transfer is calculated by the expression + = $ + .1 -/ 0 1 2 ! -2 ! 2 * , introducing a dependence to the water mass transfer. αs is the solute mass transfer coefficient between the mobile and the immobile zones (T -1 ) and the concentration c* is cim if Γw > 0 and c* is cm if Γw < 0. wf is the relative volumetric proportion of the immobile zone in the considered medium.

In equation 5, solute transport through the gas phase is considered by an additional term on the right hand side where Sg is the saturation of the gas phase in the mobile domain, θs m is the water content at saturation of the mobile domain, Dg is the pore gas diffusion coefficient (L 2 T -1 ) and cg is the solute concentration in the gas phase. The pore gas diffusion coefficient is calculated from the local saturation and porosity, according to the formula proposed by [START_REF] Millington | Gas diffusion in porous media[END_REF]:

, = 6 ! 7 8 ⁄ , : 8 ⁄ , ;
(eq.7)

where Dg 0 is the gas diffusion coefficient in a free medium. The coupling of chemical species transport as solute in the soil water and as gas species in the soil atmosphere allows the reactivity of the gases with minerals and the feeding of the gas species from the atmosphere to be considered depending on the soil saturation. This coupling is of particular interest when considering the oxidation of pyrite by atmospheric O2 since the transport of O2 in water is slow and the dissolved O2 is rapidly consumed by pyrite. Without O2 gas transport the pyrite oxidation is restricted to the first cm of the soil column.

Water retention in the mobile and immobile domains was defined by the van Genuchten-Mualem relationships [START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils[END_REF], linking water content and hydraulic head in a porous medium and calculating the effective hydraulic conductivity. They are written as follows:

< = = > ( = > = 7 7?|A | B < 0 (eq.8.a) < = 1 ≥ 0 (eq.8.b) F = 1 - 7 G 6 > 1 (eq.8.c) < = + < I [1 -K1 -< L M ! ] O (eq.8.d)
where Se is the effective saturation, θs is the water content at saturation, θr is the residual water content, α, n and l are empirical parameters of the retention functions and Ks is the hydraulic conductivity at saturation.

Reactive model

The reactivity considered in the model is based on the alteration mechanisms identified in the COx tailings ( § 3.2). Major (Ca, Na, Mg, K, Fe, Al, S, C, Si, Cl) and trace elements (Pb, Co, Ni, Zn) measured in the lysimeter solutions and involved in the reactive processes were considered as solute in the model. O2 and CO2 were considered as gas species that can be dissolved in water. A selection of primary phases (pyrite, calcite and goethite, such as iron (oxy-)hydroxide) can dissolve or precipitate through a kinetic control, according to the kinetic parameters selected by [START_REF] Marty | A database of dissolution and precipitation rates for clay-rocks minerals[END_REF].

Thermochimie v.9b0 [START_REF] Giffaut | Andra thermodynamic database for performance assessment: ThermoChimie[END_REF] was used as the thermodynamic database for these calculations. The composition of the considered pyrite phase includes Pb, Ni, Co and Zn trace metals at proportions given by [START_REF] Lerouge | Mineralogical and isotopic record of biotic and abiotic diagenesis of the Callovian-Oxfordian clayey formation of Bure (France)[END_REF]. Gypsum completes the selected mineralogical assemblage as secondary phases whose occurrence is controlled by thermodynamic equilibrium. For the sake of simplicity, clay minerals were not introduced because of their slow dissolution rate compared to the simulation duration. Cation exchange involving Na + , Ca 2+ , Mg 2+ , K + and Fe 2+ which occurs at the surface of clay minerals was nevertheless considered given its influence on the solution composition during transient chemical and hydraulic conditions [START_REF] De Craen | High sulphate concentrations in squeezed Boom Clay pore water: evidence of oxidation of clay cores[END_REF]Jacques et al., 2008b). The cation exchange capacity varying between 13.3 and 18.3 meq kgrock -1 was measured by cobaltihexammine extraction on tailings sampled at different depths at the installation of lysimeters L5 and L6. The cation exchange selectivity coefficients determined for the pristine COx claystone [START_REF] Gaucher | A robust model for pore-water chemistry of clayrock[END_REF] were used. Sorption of trace metals can take place at the surface of goethite minerals that locally precipitate in the vicinity of altered pyrite crystals and at the surface of clay minerals. Retention coefficients were fitted on the results for each trace element Co (log K = 10), Ni (log K = 9.7), Pb (log K = 11.7) and Zn (log K = 9.9), considering a site density of 0.28 mol of sites per mol of goethite [START_REF] Dixit | Comparison of Arsenic(V) and Arsenic(III) Sorption onto Iron Oxide Minerals: Implications for Arsenic Mobility[END_REF]. K is the half-reaction selectivity constant of the mass-action equation describing the distribution of species between the exchange sites and the solution, as defined in PHREEQC code.

2.4.3.

Hydraulic and chemical initial and boundary conditions

The hydraulic boundary conditions used in the reactive transport model corresponded to the precipitation, evaporation and drainage presented or recalculated in section 3.1.1. Precipitation and evaporation series were then applied as upper boundary conditions on the tailing column while lysimeter basal drainage was used as lower boundary condition. Hydraulic boundary conditions were defined at a daily frequency. The initial water content profile in the model was set to that measured in January 2014 on the lysimeters.

The rainwater composition was the chemical boundary condition at the lysimeter surface and was based on SO4 (0.017 mmol L -1 ) and Cl (0.021 mmol L -1 ) concentrations measured on site and on rainwater compositions reported in [START_REF] Appelo | Geochemistry, groundwater and pollution[END_REF] for the remaining composition. Since the tailings were already oxidized at the beginning of the simulation, the initial chemical system assumed for the model showed the effects of oxidation. Initial solutions in the lysimeter column corresponded to the composition of solutions at -20 cm, -50 cm, -100 cm, -150 cm, -200 cm in lysimeter L6. Since the Fe concentration was generally below the quantification limit (5 10 -8 mol L -1 ), it was determined by equilibrating with goethite. The pH of the solutions was altered until the saturation index of calcite was close to zero. These six solutions were discretized further by linear interpolation into 26 solutions in order to reduce the difference in concentration of the solution between two nodes. The cation distribution on the exchanger was calculated to be at equilibrium with the initial solutions. The initial solutions contained trace metal concentrations based on data.

Wherever data for the trace metal were not available at that particular depth at the Jan-14 L6 sampling campaign, the trace metal concentration for that depth was chosen based on the range of concentrations in solutions measured at a different sampling time.

Regarding the gas phase, the CO2 and O2 atmospheric composition was set as the initial gas composition in the unsaturated porosity of the modelled lysimeters. The initial model parameters are summarised in Supplementary file A.

Results and discussion

3.1. Water fluxes calculated from lysimeter data Due to the release of the mechanical stress, the flow properties in the tailings of Callovian-Oxfordian (COx) shale are higher compared to the properties of the sound rock for in-situ conditions. Indeed, the rock has been milled and uncompacted during the excavation operations, and has been partially

recompacted when stored at the ground surface as tailings. As a result, the porosity and the pore size have increased compared to its initial state. With an apparent density of 1700 kg m -3 after recompaction of the tailings, the porosity can be estimated at 55 to 60 %, while it was of about 10 to 20 % initially. The tailings constitute an assemblage of µm to cm blocks or aggregates of rock. Hence, the pore space of the tailings becomes highly heterogeneous with the occurrence of large spaces (macropores, fissures…) between the rock blocks and the porosity of these blocks is made of pores of very little size. Consequently, zones of preferential flow are expected in such a porous medium, as

shown by the range of hydraulic conductivity from 10 -3 to 10 -8 m s -1 measured on tailings samples (Andra, internal report). In addition, depending on the meteorological conditions the COx shale tailings are subject to variations in water saturation.

3.1.1.

Water balance in the lysimeters

Water content and drainage in the lysimeters showed a seasonal evolution with a dry period during summer. The measured daily weight variations could be irregular where a sharp weight variation was often followed by an opposite weight variation.

The lysimeters' weight seemed to increase during rain events and to decrease shortly after. On the other hand, the drainage at the base of the lysimeters seemed to be less dependent on precipitation, since heavy rains were not followed by a drainage outflow, while heavy drainage could be observed without direct links to rain events. This observation suggests a noticeable delay in the water flow within the lysimeters.

The evaporation calculated from the water balance in the lysimeters was not coherent because the expected seasonal variation with higher evaporation in summer than in winter was not found (Figure 3.a and b). In addition, negative evaporation was sometimes calculated, following days with high positive evaporation, and this inconsistency could not be solved by averaging the data over several days. The drainage recalculated using the evaporation obtained from meteorological data, with the daily values averaged over three days, stayed in the range of the drainage measured from the lysimeters (Figure 3.c), but the instantaneous volume of water drained from the lysimeters' base was lower and more spread over time.

The precipitation data, the evaporation calculated using the Penman-Monteith formula and the recalculated drainage that will be used as boundary conditions for the lysimeter flow simulations are summarized in Figure 4.

Insights on water flow in the lysimeters

Water content showed a seasonal evolution with, for the sensor at -20 cm, a strong link with the rain events. A delay was observed at depth in the seasonal water content evolution, where the wet period at -150 cm was delayed by several months compared to the wet period close to the surface.

Rock tailings at depth remained close to saturation (with a water content of 60 %). The water content measurements showed a non-negligible difference between the different lysimeters for a given depth. However, even if the water content value was different between the lysimeters, a temporal similarity was observed between the water content increases and decreases.

The suction potential is null in a saturated porous medium and increases with the desaturation of this medium. In agreement, the highest recorded suction potentials corresponded to the more pronounced desaturation events, close to the surface and during the summer period. At depth, suction potentials remained close to 0, indicating that the medium remained saturated or close to saturation. Suction potential gradients gave insights on the water flow in the lysimeters. These suction potential gradients indicated a water flow from the top to the bottom of the lysimeters during the wet season and a reverse gradient during the summer period, suggesting an upward capillary flow. The potential gradients were generally higher between the tensiometers close to the surface than between the deep tensiometers. At the lysimeter scale, the mean gradient was about 0.05 cm cm -1 but much higher (and reversed) gradients of some cm cm -1 were occasionally recorded during dry periods.

Weathering mechanisms

Elemental chemistry in solution

The concentrations measured in solution were similar in the four lysimeters. We chose to present only the results for L6 (Figure 5).

The pH varied between 7 and 8 throughout the depth of the lysimeters. The pH decrease observed with depth in January 2017 compared to the other measurements is not characteristic of an acidification from this date as such a pH decrease was also observed in November 2015. The pH values in the lysimeters were higher than the pH in the pristine COx shale pore water, before oxidation, in agreement with the carbonate ions that were lower. The data were consistent with each other as the higher the carbonate ions the lower the pH [START_REF] Millero | The Marine Inorganic Carbon Cycle[END_REF]. Nonetheless, carbonate ions concentration might be expected to have been higher because of the calcite dissolution and biological activity. However, the degassing during sampling due to the exposure of the lysimeter water to the atmosphere lowers the carbonate content in solution. The dissolved carbonate concentration was lower close to the surface, increased until 50 cm depth and then remained stable.

Cl in solution is two orders of magnitude below the values reported in the pristine COx pore water [START_REF] Vinsot | In situ characterization of the Callovo-Oxfordian pore water composition[END_REF]. The value was slightly higher than the rainwater concentration measured on site

(2.1 10 -5 mol L -1 ). This concentration supported the assumption of percolating rainwater in the lysimeter with interaction with pore water. Note that the concentration decreased from January 2014 to January 2015 and then remained stable highlighting a possible leaching of Cl that decreased with time.

The sulphate concentration increased with depth and remained in the same range throughout the study period. Sulphate levels in the lysimeter were two to three times higher than in pristine COx pore water for which redox reduced condition prevails and two orders of magnitude higher than in rainwater (1.75 10 -5 mol L -1 ). This observation is consistent with a dissolution of pyrite, oxidized by the oxygen from the atmosphere. This observation is also consistent with the gypsum saturation indices that were at equilibrium from 50 to 200 cm deep and the observation of gypsum precipitates in the lysimeters.

Na, K and Mg concentrations increased with depth in all lysimeters during the monitored period. K and Mg concentrations remained similar over time while Na concentration clearly decreased between January 2014 and January 2017. At the surface of the lysimeters, Na reached in 2017 the concentrations observed in natural soils (forest and meadow) formed on Callovian-Oxfordian shale [START_REF] Scholtus | Expected evolution of a Technosol derived from excavated Callovo-Oxfordian clay material[END_REF]. In addition, K and Mg concentrations were higher in the lysimeter solutions than in COx pore water while the opposite was observed for Na due to its leaching.

Ca increased from the surface to 50 cm deep and then decreased until the bottom of the lysimeter.

This element was not only buffered by the ion exchanger but mainly by highly reactive minerals such as calcite that controlled its concentration in the Callovian-Oxfordian [START_REF] Gaucher | A robust model for pore-water chemistry of clayrock[END_REF]. Such a concentration higher than the COx pore water could display dissolution of this phase. However, the saturation indices calculated on the basis of the solution analyses showed that calcite, like dolomite, was oversaturated and apparently ruled out this hypothesis (Figure 5). Nonetheless, the solutions were collected in the lysimeters every 2 months and were not preserved from the atmosphere during that period. The pCO2 measured in the COx is close to 10 -2 atm [START_REF] Gailhanou | Effects of a thermal perturbation on mineralogy and pore water composition in a clay-rock: An experimental and modeling study[END_REF][START_REF] Vinsot | In situ characterization of the Callovo-Oxfordian pore water composition[END_REF], in agreement with the pCO2 calculated from the solution analyses (Figure 5). CO2 degassing no doubt occurred and modified the carbonate concentration measured by the TOC analyser [START_REF] Lerouge | Origin of dissolved gas (CO2, O2, N2, alkanes) in pore waters of a clay formation in the critical zone (Tégulines Clay, France)[END_REF]. This may explain the calcite and dolomite oversaturation.

Figure 5. Variation of the major elements in lysimeter 6 during the 3 years of monitoring (January 2014(January , 2015(January , 2016(January and 2017) ) and saturation indices of calcite, gypsum and dolomite. Water at depth 0 cm corresponds to water that accumulates at the top of the lysimeter. The red line corresponds to the value in the pristine COx pore water, before any oxidation.

3.2.2.

Isotopic ratio of sulphates (δ 34 S)

All the results were reported in δ units relative to international standards, defined by: δ = (RSample/RStandard -1) x 1000‰, where R is the measured isotopic ratio in the sample (δ 34 S/δ 32 S) and in the standard, which was CDT (Canyon Diablo Troilite) for sulphur. The isotopic signature of the SO4 from the lysimeter samples ranged from -41 to -18 ‰ CDT (Table 1 andsupplementary file B). These values correspond to the pyrite domain in the pristine COx formation that ranges between -40 and -20 ‰ CDT [START_REF] Lerouge | Mineralogical and isotopic record of biotic and abiotic diagenesis of the Callovian-Oxfordian clayey formation of Bure (France)[END_REF]. Isotopic analysis of well-preserved COx pore water in anaerobic conditions shows high δ 34 S values, between 18.7 and 22.3 ‰ CDT [START_REF] Vinsot | In situ characterization of the Callovian-Oxfordian seepage water composition while exposing the rock to oxygen gas[END_REF] which was in line with the δ 34 S values of celestite [START_REF] Lerouge | Mineralogical and isotopic record of biotic and abiotic diagenesis of the Callovian-Oxfordian clayey formation of Bure (France)[END_REF]. After introducing oxygen into the preserved COx pore water samples, the SO4 concentration in the oxidized pore water increased twofold and the δ 34 S value dropped down to -16.4‰ CDT. This decrease was interpreted to be due to release of SO4 depleted in δ 34 S, i.e. diagenetic pyrite present at the deeper parts of COx (Tremosa et al., 2015;[START_REF] Vinsot | In situ characterization of the Callovian-Oxfordian seepage water composition while exposing the rock to oxygen gas[END_REF]. Based on these results, the isotopic analysis of the sampled solutions from the lysimeter indicated that the SO4 concentration of the solutions contained a substantial amount of SO4 released directly from diagenetic pyrite. These results confirm the oxidation of pyrite in the tailings, although they are not sufficient to quantify the extent of pyrite dissolution or the amount of SO4 released by the pyrite. Mineralogical characterizations of COx tailings evidenced the presence of pyrite, calcite and the precipitation of secondary phases: gypsum and iron (oxy-)hydroxides (Figure 6). The iron (oxy-

)hydroxides formation was linked to pyrite alteration and was described as a three-step process involving the oxidation of disulphide (S2 2-) and Fe 2+ that led to the formation of Fe 3+ [START_REF] Appelo | Geochemistry, groundwater and pollution[END_REF]. Finally, iron (oxy-)hydroxides precipitated unless the pH remained below 3 (Debure et al., 2017). The overall reaction is described by equation 9:

4RS O 15 U O 14V O U → 4RS UV 8 8 U Y O= 16V ?
(eq.9)

The pyrite dissolution released protons in solution and thus led to an acidification of the system, which was buffered by the dissolution of calcite towards the equilibrium pH of calcite (≈8). Indeed, the proton production caused calcite dissolution in agreement with equation 10.

[

\[U 8 2V ? → [\ O? V O U [U O (eq.10)
This calcite dissolution released CO2 and Ca that may eventually precipitate with the sulphates coming from the pyrite dissolution to form gypsum according to equation 11:

[\ O? U Y O= 2V O U → [\ U Y . 2V O U (eq.11)
Note that gypsum precipitation is highly dependent on the water saturation. Indeed lysimeters faced weather changes (rain, snow, wind …) over the years and seasons with water saturation and desaturation periods. Water loss led to an increase in the concentration. Evaporation of half the pristine COx pore water [START_REF] Gaucher | A robust model for pore-water chemistry of clayrock[END_REF] led to gypsum precipitation. Such water loss induced an increase in the sulphate concentration by a factor of 1.7 compared to the pristine pore water. This factor reached 3.5 for an evaporation of 70% of the COx pore water. This mechanism highlights the possible formation of gypsum at the top of the lysimeter that was more subject to saturation/desaturation cycles. Nonetheless, leaching and infiltration caused an increase in the sulphate and calcium concentrations from the surface towards the bottom of the lysimeters (Figure 5). Therefore, gypsum precipitated at the bottom of the lysimeter as shown by the saturation indices calculation (Figure 5) even if there was no water loss.

Trace element release in the lysimeters

In the COx shale, pyrite contains several trace elements: As, Co, Cu, Ni, Pb and Zn [START_REF] Grangeon | The influence of natural trace element distribution on the mobility of radionuclides. The exemple of nickel in a clay-rock[END_REF][START_REF] Lerouge | Mineralogical and isotopic record of biotic and abiotic diagenesis of the Callovian-Oxfordian clayey formation of Bure (France)[END_REF]. Most of these elements were found at low concentrations in the solutions collected in the lysimeters (Figure 7). Among them, As was below the quantification limit (3.2 10 -8 mol L -1 ) most of the time. As was one of the main traces incorporated in pyrite with Co, Ni and Pb, so this behaviour could be explained by the high affinity of As for iron (oxy-)hydroxides that precipitate as a result of the pyrite dissolution [START_REF] Daus | Arsenic speciation in iron hydroxide precipitates[END_REF][START_REF] Debure | Retention of arsenic, chromium and boron on an outcropping clay-rich rock formation (the Tégulines Clay, eastern France)[END_REF][START_REF] Giménez | Arsenic sorption onto natural hematite, magnetite, and goethite[END_REF]. Cu and Zn are infra-traces in pyrite [START_REF] Lerouge | Mineralogical and isotopic record of biotic and abiotic diagenesis of the Callovian-Oxfordian clayey formation of Bure (France)[END_REF] explaining the low Cu concentration while Zn was provided by sphalerite (ZnS) dissolution. Furthermore, Zn and Ni were present in feldspar and detrital minerals thus accounting for their higher quantity than Co and Cu. The values detected for Pb were explained by the presence of galena (PbS) in the pristine COx too [START_REF] Lerouge | Mineralogical and isotopic record of biotic and abiotic diagenesis of the Callovian-Oxfordian clayey formation of Bure (France)[END_REF]. With the exception of Cu concentrations that remained in the range of the COx pore water, the concentrations of Co, Ni, Pb and Zn were two orders of magnitude higher in the tailings water than in the COx pore water.

Ba was detected as well. Mica, quartz and feldspar were the Ba carrier and their dissolution is the most relevant explanation for its presence in the solution. Ba released at the bottom of the lysimeters reached a concentration of 2 10 -7 mol L -1 in March 2014 but stabilized and did not exceed 7.3 10 -8 mol L -1 after March 2015 at the bottom of the lysimeters. Regarding Se, pyrite and fossil organic matter were the Se-bearing phases in the initial COx. In sampled solutions, Se was under the quantification limit (5.3 10 -8 mol L -1 ) in the older tailings (~ 10 years in 2014) while it was still detected in the new tailings (~3 years in 2014) but decreased from 3.9 10 -7 mol L -1 in 2014 to 1.3 10 -7 mol L -1 in 2017 at the bottom of the lysimeters. This disappearance could be due to the leaching of this element or to its sorption on the iron (oxy-)hydroxides [START_REF] Jordan | Competition between selenium (IV) and silicic acid on the hematite surface[END_REF]. Mn was detected as well at the bottom of the lysimeters and varied from 5.8 10 -6 mol L -1 in the aged tailings to 1.8 10 -8 mol L -1 in the new tailings. Mn was detected especially in the solutions of the aged tailings. Mn was carried by kaolinite, micas and carbonates and MnO represented initially 0.04 wt% of clayrock. Such a Mn concentration seemed to show a stage of dissolution that had not yet occurred in the new tailings. Other elements like Cd, Mo and Sb were always below the quantification limit during the three years of monitoring (3.6 10 -9 mol L -1 , 2.1 10 -7 mol L -1 and 8 10 -8 mol L -1 respectively) as well as Cr most of the time (QL = 2.3 10 -8 mol L -1 ), which was detected only occasionally but did not exceed 1.6 10 -7 mol L -1 . Finally, Al and Fe reached their highest concentrations (2.2 10 -6 mol L -1 and 1.8 10 -7 mol L - 1

) at the bottom of the lysimeters in the first samples; they were not detected anymore after 5 months and remained below their limit of quantification (QLAl = 1.3 10 -7 mol L -1 and QLFe = 1.8 10 -8 mol L -1 ).

3.2.5.

Gases: O2 and CO2

The most important mechanism driving the COx tailing alteration was pyrite oxidation, which is highly dependent on O2 available in the lysimeters and on the exposed surface of pyrite to O2. The lysimeters were in contact with the atmosphere, thus initially O2 was at the air value and decreased with depth due to its consumption by pyrite and because of the lower accessibility of the deeper facies (permeability, saturation). However, oxidation of the water samples collected every two months was not prevented between sampling times and during sampling and, then, no redox state differences with the depth were revealed by measurements on the water. The redox conditions of the system can neither be assessed by the nitrite/nitrate couple nor the sulphide/sulphate couple.

The CO2 amount was assessed by the measurement of inorganic carbon in solution that was compared occasionally with the alkalinity measurements. Both parameters were in agreement each time. However as mentioned earlier the interaction with the atmosphere during water sampling led to CO2 outgassing that lowered the analysed CO2 compared to its quantity in the lysimeters.

3.2.6.

Synthesis on the mechanisms involved in the COx alteration and release of elements Microscope observations on tailing samples showed pyrite alteration in the tailings that finally led to the formation of iron (oxy-)hydroxides. As oxidizing conditions prevailed in the lysimeters, pyrite dissolution led to an increase in sulphates in solution which was confirmed by δ 34 S isotopic analyses.

In the meantime, calcite dissolved and buffered the system close to pH 8 and drove the increase in Ca in solution. The increasing amount of Ca and sulphates may lead to gypsum precipitation as already seen for outcropping or shallow clay-rich rock formations [START_REF] Debure | Retention of arsenic, chromium and boron on an outcropping clay-rich rock formation (the Tégulines Clay, eastern France)[END_REF][START_REF] Lerouge | A Deep Alteration and Oxidation Profile in a Shallow Clay Aquitard: Example of the Tégulines Clay[END_REF]. In addition, desaturation phenomena due to warm and dry weather conditions could enhance this secondary phase's precipitation by drying the tailings. In addition to its influence on saturation indices, water content will affect the gas diffusion, with a higher gas mobility under unsaturated conditions. CO2 migration influences the pH and the carbonate phases while O2 diffusion will be the driver of the pyrite oxidation which releases trace elements such as As, Co, Cu, Ni, Pb and Zn in solution. Such metal elements have an affinity for clays and iron (oxy-)hydroxides and their behaviour is expected to be controlled by sorption on these phases.

Modelling of lysimeter evolution and Discussion

3.3.1.

Water content and transport a) Water content

Water content coupled with pressure head were used to determine hydraulic parameters through an inversion. A single set of parameters were determined for the four lysimeters considering their similar evolution with time. The inversion criteria, by visual assessment, were based on the ability of the model to describe the water content variations and delay at depth, depending on rainfall and evaporation at the top of the lysimeters and drainage at the bottom. The lower measuring depth (-200 cm) was usually saturated and so the reproduction of the data focused on the shallower depths (Figure 8). Overall, the modelled water content reproduced the most significant measurements especially the water content increase in July and August 2014 as well as the decrease in October 2016. Finally, the higher water contents were consistent with the rainfall at the shallower depth while at depth a retardation was observed. Given the scattering of the water contents measured between the different lysimeters, the absolute water content value seemed less important to simulate than the water content variations over time at the different depths. The water content evolution was modelled for a hydraulic conductivity at saturation of 100 cm day -1 within the mobile domain, which is in the range of the hydraulic conductivities measured in tailings and in the lysimeters (Andra, pers. comm.). Water and solute transfer coefficients between the mobile and immobile zone as well as dispersivity are the parameters describing flow and transport in the double porosity model. These parameters were determined by fitting the Cl concentration over time in the lysimeters. Cl is usually used as an inert tracer (no sorption, exchange or precipitation as a mineral). The initial concentration profile in the mobile zone implemented in the model was the profile measured in January 2014 (i.e. equal to or lower than 4 10 -4 mol L -1 ). The concentrations in the immobile zone were set to the pristine COx pore water concentration, 4.1 10 -2 mol L -1 [START_REF] Vinsot | In situ characterization of the Callovo-Oxfordian pore water composition[END_REF]. The Cl in rainwater measured on site was 2.1 10 -5 mol L -1 .

Cl concentrations increased during dry seasons while they decreased during wet seasons when chlorine was diluted by the rainwater, resulting in concentrations 20 times lower than the initial chlorine concentration in the system (Figure 9). Despite the two orders of magnitude of difference between the concentration in the immobile and in the mobile zone, the Cl in the immobile zone remained unchanged throughout the simulation. Indeed, the higher Cl loss was observed at 20 cm and represented only 2.6 % of the initial Cl concentration highlighting the small transfer coefficient between the mobile and immobile zone. Chlorine in the immobile zone was initially 4.1 10 -2 mol L -1 , some hundred times higher than chlorine in the mobile zone.

c) Gas transport

Seasonal desaturation of the lysimeters enables gas infiltration throughout the whole lysimeter. In this study, gas diffusion focused on O2 and CO2 that were identified as the main sources of reactivity according to the solution and solid analyses (Section 3.2). According to the results of the model (Figure 10), the higher O2 temporal variations were reported near the surface (20 and 50 cm deep) while CO2 was not present at a depth of 20 cm and displayed similar values between 50 and 150 cm deep to finally disappear at the bottom of the modelled lysimeter with the water saturation. O2 was consumed by pyrite oxidation that limits its diffusion to larger depths, together with the fact that gas diffusion is lower when water saturation increases. CO2 was released by calcite in response to pyrite oxidation. Maximal CO2 release in the soil atmosphere occurred at 20 cm and was lower at depth because of the lower pyrite oxidation and the water saturation of the lysimeter. The range of pH was well reproduced by the RTM model at the different depths of the lysimeter profile. Nonetheless, the modelled pH was slightly lower than the measured pH at 50 cm deep and increased at the surface (20 cm deep), because of CO2 outgassing. At larger depth, the modelled pH was higher because of the lower reactivity of the rock. At 150 cm deep, the variations of pH over time were related to the presence or not of a gas phase and when O2 can enter at this depth, pyrite was oxidized and pH was lower than under saturated conditions. Under conditions closed to the atmosphere, the oxidation of the tailings and the dissolution of calcite led to in-situ pH in the COx shale tailings of about 9 [START_REF] Marty | Weathering of an argillaceous rock in the presence of atmospheric conditions: a flow-through experiment and modeling study[END_REF]. Lower pH of 7.5 was obtained by equilibration with the atmosphere when air entered in the tailing porosity or during sampling.

K and Na concentrations in solution were only controlled by the cation exchange at clay surfaces since we considered no minerals other than pyrite, calcite, and secondary gypsum and goethite in the model. Both elements increased at depth because of leaching from the surface of the lysimeter, where weathering is more pronounced, and an accumulation of solute elements in the deeper part of the lysimeter. It is clear for Na at 150 cm depth that the displayed concentration decreased over the time (Figure 12). Such leaching was not observed for Ca and SO4, which displayed seasonal variations at 50 cm, well reproduced by the model. The measurements as well as the model reported fewer fluctuations at 150 cm deep (Figure 13). The variations depended on the water content that influences element solubility and thus mineral dissolution and precipitation. Pyrite in the model decreased by 0.10 wt% at 50 cm deep and by 0.02 wt% at 150 cm deep. This dissolution was congruent with calcite dissolution that was greater at 50 cm than at 150 cm (Figure 14). The dissolutions were consistent with the O2 availability in the lysimeter that enabled greater pyrite oxidation close to the surface than at depth (Figure 10). In addition, pyrite oxidation released SO4 in solution while calcite released Ca. The two elements were seasonally able to form gypsum when the water desaturation was enough to allow for gypsum oversaturation. This phenomenon occurred at shallow depths, down to 100 cm deep, while the water content was too high below to permit gypsum precipitation. Ferrous iron released by pyrite dissolution was oxidized and continuously precipitated as goethite. was simulated together with their sorption on clays and iron (oxy-)hydroxide phases. The behaviour of trace metals was well described by the model at 50 cm deep while Pb was overestimated at 150 cm deep (Figure 15 and Figure 16). This overestimation is due to the saturation of the sorption site at this depth. It is worth noting that the model satisfactorily described the trace metals' behaviour despite their low amount compared to major elements, suggesting an appropriate description of the different reactive processes within the tailings.

From these modelled results, it was possible to determine a solid/liquid distribution coefficient (RD) that varies with time and depth and is corrected from water content variations according to equation 12:

_ = `(ab c `(abde aB (eq.12)

where Csolution and Csolid are, respectively, the element concentrations at a given time in the solution (mol L -1 ) and on the solid (mol kg -1 ).

Because of the water content and concentration variations with time, the distribution coefficients also varied with time in a large range (Table 2). Nonetheless, the distribution coefficients of each element were within the range of the values determined from laboratory experiments (Andra, pers.

comm.) for Co (225 L kg -1 ), Ni (160 L kg -1 ) and Zn (2000 L kg -1 ) but exceeded the value given for Pb (1000 L kg -1 ). However, these values were determined for saturated conditions. Desaturation as observed in the lysimeter had a high impact on the distribution coefficients that were lower at the bottom of the lysimeter, under saturated conditions. This can be related to the solute concentration in water which increases after evaporation, increasing consequently the sorption (Jacques et al., 2008a, b). The water content and its variation are finally the key parameters controlling the whole reactivity of the system. Indeed, saturation state rules O2 transport through the lysimeters and therefore pyrite oxidation. This study confirms the importance of considering variable unsaturated water flow and the transport of reactive gas to simulate the weathering of rock piles [START_REF] Mayer | Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage[END_REF][START_REF] Molins | Coupling between geochemical reactions and multicomponent gas and solute transport in unsaturated media: A reactive transport modeling study[END_REF]. Consequently, Fe, SO4 and protons were released in solution. Protons acidified the lysimeter but the carbonates (mostly calcite) buffered the system by dissolving [START_REF] Debure | Study of Iron-Bearing Dolomite Dissolution at Various Temperatures: Evidence for the Formation of Secondary Nanocrystalline Iron-Rich Phases on the Dolomite Surface[END_REF][START_REF] Mayer | Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage[END_REF][START_REF] Ouangrawa | Reactive transport modelling of mine tailings columns with capillarity-induced high water saturation for preventing sulfide oxidation[END_REF] and released Ca in the solution. This chemical evolution recorded in the lysimeters was corresponding to the one identified during shale weathering [START_REF] Littke | Quantification of loss of calcite, pyrite, and organic matter due to weathering of Toarcian black shales and effects on kerogen and bitumen characteristics[END_REF] and in laboratory oxidation experiments on COx shale tailings [START_REF] Marty | Weathering of an argillaceous rock in the presence of atmospheric conditions: a flow-through experiment and modeling study[END_REF]. A decrease of the water content increases the element concentration in solution and enables the precipitation of secondary phases such as gypsum and goethite. The modelled pyrite dissolution was higher near the surface (10 wt% of dissolution) than at the bottom (1 wt% of dissolution) of the lysimeter because of the higher O2 availability. Therefore, goethite precipitation was higher near the surface than at the bottom, which highlighted the spatial difference in the lysimeter. A weathering decreasing with depth was expected, in link with the availability in O2, its consumption and its resupply [START_REF] Brantley | Probing deep weathering in the Shale Hills Critical Zone Observatory, Pennsylvania (USA): the hypothesis of nested chemical reaction fronts in the subsurface[END_REF][START_REF] Littke | Quantification of loss of calcite, pyrite, and organic matter due to weathering of Toarcian black shales and effects on kerogen and bitumen characteristics[END_REF][START_REF] Vriens | Poregas distributions in waste-rock piles affected by climate seasonality and physicochemical heterogeneity[END_REF]. Besides influencing the reactivity, water content also induces water composition changes by dilution, as previously observed by Jacques et al. (2008b).

Transport of trace metals release from the COx tailings minerals (mainly sulphurs) is limited thanks to their sorption on clay minerals and on secondary minerals. Iron (oxy-)hydroxide phases and particularly goethite are well known to complex heavy metal on their surface [START_REF] Debure | Retention of arsenic, chromium and boron on an outcropping clay-rich rock formation (the Tégulines Clay, eastern France)[END_REF][START_REF] Dzombak | Surface complexation modeling. Hydrous ferric oxide[END_REF][START_REF] Lynch | Environmental Risk of Metal Mining Contaminated River Bank Sediment at Redox-Transitional Zones[END_REF][START_REF] Schultz | Adsorption and desorption of metals on ferrihydrite: reversibility of the reaction and sorption properties of the regenerated solid[END_REF]. The prediction of the trace metals' behaviour in the tailings could be improved by dedicated sorption data.

The geometry and the thickness of the tailings could also have an influence on the sulphate content in the leakage water since pyrite oxidation is lower at depth and, then, a thick tailing could dilute the shallow release of SO4. It is worth noting that the model was able to reproduce SO4 evolution during the investigated period might thus be used to predict its behaviour over a longer period of time.

e) Uncertainties of the model and variability between lysimeters

The model was adjusted to describe the measured water content, fluxes and transport well, but no fitting was done for the geochemical part. However, from one lysimeter to another the concentrations of elements varied by a factor of 2. This difference was partly due to the hydraulic conductivity that is different between the lysimeters. Nevertheless, this hypothesis is not sufficient to explain the variations in concentration since a change in hydraulic conductivity by several orders of magnitude would be necessary to overcome this difference. The presence of cracks and preferential pathways cannot be discarded particularly because of the size heterogeneity of the lysimeter filling material.

Another source of variability may be the kinetic parameters chosen for pyrite and especially the reactive surface area. The same surface area was considered for the entire profile while the alteration was greater near the surface than at depth. The pyrite nanoparticles are major contributors to reactivity and dissolved first, lowering the reactive specific surface [START_REF] Marty | Weathering of an argillaceous rock in the presence of atmospheric conditions: a flow-through experiment and modeling study[END_REF]. This reactivity may partly explain why the modelling curves were in the upper part of the data.

Thus, the older lysimeters should be less reactive. However, the lysimeter with aged COx (L5 and L6) did not release less SO4 than the lysimeter filled with more recent tailings (L7 and L8).

In addition to the sorption model with a global surface site that was used here, the [START_REF] Dzombak | Surface complexation modeling. Hydrous ferric oxide[END_REF]) surface complexation model was tested to describe the retention of trace elements.

Although the modelling results were consistent for Co, Ni and Zn, the complexation constants of Pb were too high and led overestimations of the sorbed Pb on the shale by a factor of 1000.

Complexation constants were usually determined in the laboratory under well-controlled and shortterm experimental conditions that can explain why they can be different from complexation constants determined from field data.

Conclusion

Flow, transport, gas diffusion and reactivity in the critical zone were investigated through a lysimeter and modelling study of the Callovian-Oxfordian shale deposited as tailings. The pyrite reactivity inside the shale tailing was highly dependent on the water saturation and the distribution of water over the lysimeter height. The water saturation slowed down the downward migration of O2 gas inside the lysimeter and therefore pyrite oxidation. As the lysimeters were desaturated from the surface, preferential weathering was observed in the upper layers compared to the deeper parts. Calcite dissolution stemmed from pyrite oxidation and buffered the system, preventing acid drainage at the outlet. The CO2 produced from this reaction then released to the atmosphere. Pyrite and calcite dissolution led to gypsum and iron (oxy-)hydroxide precipitation as the main secondary phases. The two phases were an indicator of the alteration and of the water desaturation (gypsum). During its alteration, pyrite released sulphates and heavy metals that finally sorbed onto goethite and clay minerals. Such sorption is efficient for maintaining low concentrations of trace metals in drainage water and then to strongly limit their release in the environment. However, high sulphate concentrations (up to 5 10 -2 mol L -1 ) remained in the drainage water.

The difference of scale between lysimeters and tailings is not likely to affect the hydrodynamic and geochemical parameters determined for the lysimeters to the same extent as the discrepancies observed between laboratory and field measurements. The lysimeter parameters are immediately applicable to tailings in order to predict metal and sulphate releases in the surface and groundwater.

However, the geometry and the landscape design of the tailings may highly affect the elemental releases especially due to runoff and water travel time in the lysimeter. A low runoff will limit the transport of suspended matter and promote infiltration inside the tailings. By considering a more representative geometry of tailings in the predictive calculations, it would be possible to accurately quantify the oxidation front in the tailings. In addition, it could be used to predict the migration and dilution in the deeper non-altered part of the tailings of elements such as sulphates released by pyrite oxidation. The sulphate and metal concentrations in the effluent water at the base of the tailings could then be estimated for environmental evaluation.

This study resulted in the quantification of the mineralogical evolution of COx shale and the drainage of major and trace elements in water in atmospheric conditions, accounting for saturation and desaturation that can also be extrapolated using reactive transport modelling to quantify the weathering and evolution of clay-rich formations at outcrops or in sulphide-bearing tailings.

Figure 1 .

 1 Figure 1. Diagram and illustration of the principle of the lysimeters for the collection of runoff water, water inside the lysimeter and water leaching through the lysimeter.

Figure 2 .

 2 Figure 2. Schematic representation of the lysimeter reactive transport model in unsaturated conditions.

Figure 3

 3 Figure 3: a) Evaporation calculated from the water balance for each lysimeter, given by Meteo France (Houdelaincourt station) and calculated by the Penman-Monteith formula using Osne-le-Val weather station data; b) evaporation showed as cumulated evaporation; c) Cumulated basal drainage measured for the four lysimeters and calculated from the water balance.

Figure 4 :

 4 Figure 4: Daily precipitation recorded by the local weather station from January 2014 to April 2017, and cumulated precipitation, drainage and evaporation in the lysimeters considered as boundary conditions for the model.

Figure 6 .

 6 Figure 6. Highly oxidized framboïdal pyrite, containing gypsum; iron (oxy-)hydroxides and preserved calcite were identified as well;. (Zoom x100). Sample EST50652 (VER1014; 2,00m).

Figure 7 .

 7 Figure 7. Trace elements released in the lysimeter over the 3 years of monitoring. As was below the quantification limit (QL) = 3.2 10 -8 mol L -1 . Co, Pb and Zn were below their QL in the COx pore water (QLCo = 8.5 10 -10 mol L -1 , QLPb = 2.4 10 -10 mol L -1 , QLZn = 7.7 10 -9 mol L -1 ).

Figure 8 .

 8 Figure 8. Comparison of the water content measured at different depths in the lysimeters to the water content modelled with HP1. b) Chlorine transport

  Sensitivity tests made on the solute mass transfer coefficient revealed that only a small range of solute mass transfer coefficient values could reproduce Cl variations. In addition, the contribution of the excess Cl content from the porosity of claystone aggregates immobile zone to the Cl content in the inter-aggregate porosity of the tailings mobile zone.

Figure 9 .

 9 Figure 9. Chlorine evolution in the mobile zone over time in the four lysimeters (L5, L6, L7 and L8), at different depths (20 cm, 50 cm, 100 cm, 150 cm and 200 cm), compared to the chlorine modelled with HP1 and chlorine loss in the immobile zone.

Figure 10 .

 10 Figure 10. Calculated evolution of gases (O2 and CO2, in mol per L of tailing material) in the lysimeters at 50 and 150 cm deep.

Figure 11 .

 11 Figure 11. Evolution of pH and bicarbonates in the four lysimeters (L5, L6, L7 and L8) at 50 and 150 cm deep. Symbols represents measurements in the lysimeter and solid lines the results of HP1 modelling.

Figure 12 .

 12 Figure 12. K, Na and Ca evolution with time in the four lysimeters (L5, L6, L7 and L8) at 50 and 150 cm deep. Symbols represents measurements on the lysimeter solution and solid lines the results of HP1 modelling.

Figure 13 .

 13 Figure 13. Ca and SO4 evolution with time in the four lysimeters (L5, L6, L7 and L8) at 50 and 150 cm deep. Symbols represents measurements on the lysimeter solution and solid lines the results of HP1 modelling. b) Mineralogical evolution

Figure 14 .

 14 Figure 14. Evolution of pyrite, calcite and gypsum modelled by HP1 at 50 and 150 cm deep. The data are represented as weight percent variations.

Figure 15 .

 15 Figure 15. Evolution of trace elements (Co, Ni, Pb and Zn) with time in the four lysimeters (L5, L6, L7 and L8) at 50 cm deep. Symbols represent measurements on the lysimeter solution and solid lines the results of HP1 modelling.

Figure 16 .

 16 Figure 16. Evolution of trace elements (Co, Ni, Pb and Zn) with time in the four lysimeters (L5, L6, L7 and L8) at 150 cm deep. Symbols represent measurements on the lysimeter solution and solid lines the results of HP1 modelling.

Table 1

 1

	Sample	Depth (cm)

. Sulphate isotopic measurements on lysimeter 6 (aged tailings) sampled in

March 2016 and November 2017 (CDT: Canyon Diablo Troilite)

. The water collected at the surface formed no precipitates when BaCl2 was added to the solution (i.e. no sulphates in solution).

Sulphate δ 34 S‰ CDT (+/-0.3‰) March 2016 Sulphate δ 34 S‰ CDT (+/-0.3‰) November 2017

		L6-0	0	-23.8	-37.3
		L6-1	-20	-35.0	-37.8
		L6-2	-50	-41.2	-38.6
		L6-3	-100	-37.8	-37.7
		L6-4	-150	-32.6	-33.7
		L6-5	-200	-34.1	-34.9
	3.2.3.	Mineralogical changes

Table 2 .

 2 Solid/liquid distribution coefficients calculated for Pb, Co, Ni and Zn from the modelling results obtained with HP1.

	Depth (cm)	RD Pb (L

kg -1 ) RD Co (L kg -1 ) RD Ni (L kg -1 ) RD Zn (L kg -1 )

		min	2833	57	28	45
	50	max	96010	1916	960	1522
		average	40800	814	408	647
		min	2793	56	28	44
	150	max	25018	499	250	396
		average	12152	242	122	193

Acknowledgements

This study was supported by the French Radioactive Waste Management Agency (Andra) and the French Geological Survey (BRGM) in the framework of their scientific partnership. LIEC (Laboratoire Interdisciplinaire des Environnements Continentaux -UMR 7360 Université de Lorraine/CNRS), LSE (Laboratoire Sols et Environnement -UMR 1120 Université de Lorraine/INRA), Laboratoire GeoRessources (UMR 7359 Université de Lorraine/CNRS), LRGP (Laboratoire Réactions et Génie des Procédés -UMR 7274 Université de Lorraine/CNRS) and GISFI (Groupement d'Intérêt Scientifique sur les Friches Industrielles) are acknowledged for providing data from the lysimeters.

ERM (Etudes Recherches Matériaux) is aknowledged for their mineralogical analysis and observations on tailing samples. The authors are grateful to the two reviewers and the editor for their fruitful comments, remarks and proposals that helped to improve the manuscript.