# Leaching of pollutant metals (Pb, Zn) from abandoned mine tailings: A multicomponent reactive transport model of a pilot-scale experiment

Samuel Mertz, Nicolas Devau, Hugues Thouin, Fabienne Battaglia-Brunet, Marie-Paule Norini, Marc Crampon, Lydie Le Forestier

<figure><img src="image_1.png" /></figure>

## To cite this version:

Samuel Mertz, Nicolas Devau, Hugues Thouin, Fabienne Battaglia-Brunet, Marie-Paule Norini, et al. Leaching of pollutant metals (Pb, Zn) from abandoned mine tailings: A multicomponent reactive transport model of a pilot-scale experiment. Science of the Total Environment, 2025, 960, pp.178248. 10.1016/j.scitotenv.2024.178248. insu-04887690

HAL Id: insu-04887690

https://insu.hal.science/insu-04887690v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à 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.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Leaching of pollutant metals (Pb, Zn) from abandoned mine tailings: A multicomponent reactive transport model of a pilot-scale experiment

Samuel Mertz<sup>a,b</sup>, Nicolas Devau<sup>a</sup>, Hugues Thouin<sup>a</sup>, Fabienne Battaglia-Brunet<sup>a,b</sup>, Marie-Paule Norini<sup>a</sup>, Marc Crampon<sup>a</sup>, Lydie Le Forestier<sup>b,*</sup>

## HIGHLIGHTS

* A multicomponent reactive transport model was built using data from a previous pilot-scale experiment performing leaching on former tailings.

* The main process controlling the Pb/Zn reactivity was the dissolution/precipitation of Pb/Zn-bearing mineral phases.

* The geochemical model previously developed at centimetric scale has proved effective in describing the behavior of Pb and Zn at metric scale.

* The release of metal contaminants from tailings continues beyond a century so natural attenuation is not a realistic solution for remediation.

## GRAPHICAL ABSTRACT

<figure><img src="image_5.png" /></figure>

## ARTICLE INFO

Editor: Frederic Coulon

Keywords:

Lead

Zinc

Sorption

Dissolution – Precipitation

Vadose zone

## ABSTRACT

Mine tailing deposits pose a global problem, as they may contain metal contaminants in various geochemical forms and are likely to be leached from the surface into the underlying groundwater, which can result in health and/or environmental risks. Unfortunately, little is currently known regarding the water flow and mass balance related to leaching in the vadose zone as these factors are still difficult to measure at the field scale. A pilot-scale experiment was run in a 1 m³ instrumented column for 6 months to address this issue. This 70 cm high column was filled with highly Pb-contaminated tailings and then watered regularly. The top half remained unsaturated while the bottom half was kept saturated. Continuous water flow and water saturation measurements were recorded and the physico-chemical properties of the soil solutions were monitored weekly at different levels in the column. A 1D multicomponent reactive transport model was built to simulate the fate and transport of Pb as well as other elements. The variably saturated water flow was simulated using the Richards equation, while the van Genuchten analytic form was used to describe the unsaturated soil hydraulic properties. The main processes considered to simulate the reactive transport were advection-dispersion, thermodynamic equilibrium, and kinetically-controlled dissolution-precipitation reactions. The most reactive Pb-bearing phases accounted for in the simulation were anglesite (PbSO₄) and plumbojarosite Pb₀.₅Fe₃(SO₄)₂(OH)₆. The simulations accurately reproduced the water flow and mass balance as well as the drop of 2 pH units experimentally measured in the pore solution. This trend resulted from plumbojarosite dissolution followed by ferrihydrite precipitation. The

<footnote><sup>a</sup> BRGM, BP 36009, 45060 Cedex 2 Orléans, France</footnote>

<footnote><sup>b</sup> Univ. Orléans, CNRS, BRGM, ISTO, UMR 7327, F-45071, Orléans, France</footnote>

<footnote>* Corresponding author.
E-mail address: lydie.leforestier@univ-orleans.fr (L. Le Forestier).</footnote>
increased Pb concentration in the soil solution induced by the dissolution of Pb-bearing phases was partially offset by Pb-sorption onto newly formed iron oxides and the precipitation of secondary mineral phases, e.g. anglesite. The modeling results could be used to assess potential risks of groundwater contamination by mine tailings.

# 1. Introduction

Metal mining is an ancient practice dating back several thousand years. This activity increased in importance over time, especially during the 19th century industrial revolution (Harris, 1988). Mining currently represents an essential industrial sector, with annual world production of certain metals reaching millions or even billions of tons including 4.5 Mt. of Pb, 13 Mt. of Zn and 2.5 Gt of Fe (Bernhard and Reilly, 2020). However, nearly 5 Gt of tailings are generated each year via these mining activities (Lu and Wang, 2012). Tailings consist of a heterogeneous mixture of fine particles from the deposit, the host rock or the ore refining process. The amounts of metals in these tailings are not economically viable, but may have environmental impacts. These contaminants can ultimately be dispersed throughout the environment, thereby threatening humans and the environment (Plumlee and Morman, 2011; Csavina et al., 2012; Atanacković et al., 2013). There are three major dispersion processes: wind erosion, water erosion and leaching. The influence of these processes on the extent of mobility of metals is dependent on many parameters such as the average slope of the site, the texture of the tailings, the vegetation cover and the climatic conditions (Chaney and Swift, 1984; Cerdà, 1998; Morbidelli et al., 2018). Regarding the leaching process, the geochemical properties of the tailings have a critical role in the metal contaminant transfer process.

The tailings mineralogy is a parameter influencing the mobility of metals and is dependent on the degree of mineral weathering and therefore the tailings age. Metal ore deposits are often associated with sulfide accessory minerals such as pyrite (FeS₂). Mining deposit leads to the excavation of these minerals and their storage on the surface. Pyrite exposure to atmospheric oxygen leads to its dissolution by oxidation and H⁺ release in solution. This acidification process, or so-called acid mine drainage (AMD), generally takes place in recent tailings containing high quantities of sulfides. These tailings are characterized by acidic conditions (pH 2–5) which induce high mobility of metals. The effects of AMD associated with recent tailings on metal mobility have been widely studied and discussed in the literature (Evangelou and Zhang, 1995; Akcil and Koldas, 2006; Simate and Ndlovu, 2014), while fewer studies have focused on the processes involved in older tailings. The latter are typically oxidized. Sulfide oxidation first occurs close to the surface, before progressively going deeper. Blowes and Jambor (1990) defined a sulfide alteration scale ranging from 0 to 10 to assess the progress of tailing oxidation. From 0 to 1, only a few pyrite grains are weakly altered along the edges, whereas only traces of pyrite remain at a scale of 10. Greater alteration of older tailings will generate a substantial amount of dissolved SO₄²⁻. The mobility of metals in this type of tailings is controlled by sorption and/or precipitation/dissolution of secondary phases. Iron derived from sulfide dissolution can lead to the precipitation of iron as iron oxide, leading to the adsorption of the metals (McCarty et al., 1998; Moncur et al., 2005; Gemici, 2008). In the case of Pb, sulfate bearing phases such as anglesite or plumbojarosite have often been detected in altered tailings (Blowes and Jambor, 1990; Moncur et al., 2005; Hayes et al., 2012; Courtin-Nomade et al., 2016). Sequential extractions carried out in tailings with different degrees of alteration showed that the proportion of Pb in the oxidizable fraction was lower in older tailings, suggesting that dissolution of the primary sulfide phases had occurred by oxidation (Anju and Banerjee, 2010). The degree of alteration of the tailings, and therefore by extension the age of this waste material, has a marked impact on the mobility of metals. The processes that are relevant for recent tailings are not necessarily so for older

tailings. It is therefore important to conduct in-depth studies on the behavior of metals in these systems, particularly during the leaching of these ancient tailings.

Several approaches have been developed for this purpose. A first approach involves the use of batch leaching tests for quick small-scale studies on several conditions (e.g. pH, solid/liquid ratio, duration, concentration) (Lee et al., 2012; Akhavan and Golchin, 2021). However, these protocols do not allow for solution renewal. The use of centimetric- to decametric-scale microcosms allows researchers to closer replicate real conditions via solution renewal (Jurjovec et al., 2002; Hiller et al., 2016; Thouin et al., 2019), but the spatiotemporal scales of these experiments remain well below those prevailing in the field. Another approach involves in situ experimentation using lysimeters, for instance, to obtain a direct record of the system dynamics (Shaw et al., 1998; Sand et al., 2007). However, several limitations are inherent to this type of experimentation. The latter is dependent on meteorological conditions for the water supply, which can affect the representativeness of the observations. Moreover, a single lysimeter is generally not sufficient to be able to gain a global understanding of the site. Otherwise, setting-up and implementing a lysimeter network is extremely expensive, so this is relatively unfeasible for large-scale studies. Each of these different approaches has its limitations that have to be taken into close account when assessing the results. The metric-scale mesocosm experimental approach offers a relevant trade-off for studying physical, chemical and biological processes involved in metal leaching at spatiotemporal scales close to those found in the field, while enabling precise and regular monitoring of water and mass fluxes (Thouin et al., 2017, 2018). Unlike centimetric- or decametric-scale, this experimental approach commonly involves the use of a metric to multi-metric column to be able to set up controlled recharge and drainage conditions and control the proportion between the unsaturated and saturated zones. Another advantage of this protocol is that it enables precise measurement of water and mass fluxes at the column inlet and outlet, so as to obtain a precise water and mass balance. However, this type of experimental protocol also has shortcomings, mainly related to the difficulty of identifying and quantifying the processes controlling metal leaching within the column. Moreover, as the experimental duration is limited (~1 year), long-term prediction of the fate of the system is difficult to achieve, hence other tools are needed to overcome these limitations.

Among these tools, an integrated modeling approach involving the use of a flow model under variable saturation conditions combined with a reactive transport model (RTM) would be relevant (Steefel et al., 2005). RTM models involve the coupling of two domains, i.e. integration of mechanistic descriptions of the different physical and biogeochemical processes controlling the global reactivity of the system, combined with a rigorous description of flows in complex media (variably saturated porous media) and of solute transport. This method has the advantage of being able to identify the different processes involved (physical, chemical and biological), and more importantly to account for their influence on a given property. Over the last 20 years, RTM models have been regularly used to simulate experimental data to gain insight into the dynamics of metals in tailings (Jurjovec et al., 2004; Acero et al., 2009; Ouangrawa et al., 2009). They also help to understand and quantify the long-term fate of the system, as to enhance management of these former mining sites (Muniruzzaman et al., 2018). These models have, however, not yet been widely applied to assess metal leaching dynamics in old tailings, while taking the continuum of unsaturated and saturated zones into account.

The present study was designed to identify and quantify water and
solute fluxes through former mine tailings in the unsaturated and satu-
rated zone continuum. An RTM was thus developed to simulate and
interpret the results of a leaching experiment performed at metric scale
(Thouin et al., 2022). This RTM combines a geochemical model with a
flow model under variable saturation conditions, along with a solute
transport model. The system processes and properties were identified
based on a previous study (Mertz et al., 2021), involving leaching ex-
periments carried out in microcosms with the same tailings. The com-
bination of results of a metric-scale experiment conducted under
controlled conditions (Thouin et al., 2022) with mechanistic modeling
of the processes should generate an accurate picture of the leaching
dynamics of Pb and Zn in former mine tailings.

## 2. Material and methods

2.1. Tailings

Tailings were sampled at a former Pb—Ag mine (Roure-les-Rosiers
site, Pontgibaud mining district, Puy-de-Dôme, France). The main
characteristics of these tailings were outlined in **Thouin et al. (2019):**
granulometry mainly composed of sandy particles (91 %), low organic
matter level (0.09 %), and a contamination by Pb (26,432 mg kg⁻¹), As
(1134 mg kg⁻¹), Ba (1063 mg kg⁻¹) and Zn (265 mg kg⁻¹). The other
main components were Si (37.86 wt%), Fe (0.52 wt%) and S (0.62 wt%).
The tailings mostly consisted of quartz and feldspar (orthoclase), and the
identified Pb-bearing phases were anglesite PbSO₄ and beudantite
PbFe₃(AsO₄)(SO₄)(OH)₆, a solid solution between plumbojarosite
Pb₀.₅Fe₃(SO₄)₂(OH)₆ (sulfate pure pole) and segnitite PbFe³⁻
³AsO₄(AsO₃OH)(OH)₆ (arsenate pure pole) (Thouin et al., 2019; Mertz

et al., 2021).

2.2. Pilot experiment

A laboratory pilot-scale device was set up to assess the mobility of
metallic contaminants (Pb, Zn) through old mine tailings (Thouin et al.,
2022) (Fig. 1). It consisted of a stainless steel instrumented column (100
cm dia. and 115 cm height) kept at a constant temperature of 23 ± 1 °C.
A more detailed description of the structure and characteristics of the
laboratory pilot-scale device may be found in Thouin et al. (2017, 2018).
The column was filled with 814 kg of homogenized tailings from the <2
cm fraction. Three depths (D1: 26.25 cm, D2: 43.75 cm and D3: 61.5 cm,
measured from the tailing surface) were instrumented: (i) with pore
water samplers (inert polytetrafluoroethylene/quartz porous probes;
three probes per depth, regularly spaced) to subsequently measure
physico-chemical parameters such as pH, Eh and electric conductivity,
and to analyze major ions and metallic contaminants; and (ii) with time
domain reflectometry (TDR; TRIME-PICO32, IMKO) probes to measure
soil moisture and temperature (Fig. 1). The lower half of the column,
including the D2 and D3 levels, was always saturated, while the upper
half, including D1, remained in unsaturated conditions. A layer of gravel
and a geotextile membrane were placed at the base of the column to
retain fine particles. The leachate (Lch) was collected and quantified in a
steel tank equipped with a balance at the exit of the column to calculate
the drainage water flows.

A first 3-month experimental period was devoted to the physical
stabilization (settlement) of the tailings (t<sub>3 months</sub> to t<sub>0</sub>). No watering was
done during the first month, while during the following 2 months water
was added several times, corresponding to a total of 71 L added. The

<figure><img src="image_6.png" /><figcaption>**Fig. 1.** Design of the pilot experiment with the instrumented device.</figcaption></figure>
saturated zone was setup using a U-tube that connects the stainless steel column outlet to the steel tank used for weighing the leachate. Mont Roucous® water was chosen because its chemical composition is close to that of rainwater (Table 1). After 3 months of stabilization, the tailings settled by 12.5 % in height, and the thickness of the tailings in the column reached a constant value of 70 cm. During the second 6-month period (t₀ to t₆ months), rainfall was simulated using a sprinkler system connected to a water reservoir and fed by pumps, with watering (5.5 L of Mont Roucous® water per application) carried out three times a week.

## 2.3. Reactive transport model

The reactive transport model was described using the HP1 code. This calculation code is based on the coupling of HYDRUS-1D 4.17 (Simünek et al., 2008) for one-dimensional solute flow and transport with the PHREEQC version 3 geochemical code (Parkhurst and Appelo, 2013). The geochemical reactions considered were: (1) aqueous speciation of solutes, (2) thermokinetic dissolution-precipitation reactions of primary and secondary mineral phases, and (3) surface complexation reactions on iron oxides. Most of the thermodynamic data used in the model were taken from the Thermoddem database (Blanc et al., 2012). The methodology used corresponded to that of a model based on continuous formulation of equations describing the physical, chemical and biological processes taking place on a small scale (micro- to nano-metric). This approach involved averaging the system properties within a representative elementary volume (Steefel et al., 2005). Combining these two codes allowed us to associate a mechanistic geochemical model with a flow and solute transport model in a variably saturated medium. In a previous study (Mertz et al., 2021), a multicomponent reactive model was developed to simulate and forecast the leaching of metal pollutants (Pb, Zn) through microcosms filled with the same tailings. This model accounted for the following geochemical processes: kinetically-controlled dissolution and precipitation reactions, sorption reactions (i.e. surface complexation reactions) and water-gas interactions. The initial conditions and geochemical processes considered for this RTM were configured to be as close as possible to those used in previous simulations to reproduce the leaching experiments at a centimetric scale.

### 2.3.1. Initial conditions

The leaching experiment was simulated by considering a one-dimensional column of 70 cm height, divided into 100 cells, with each corresponding to a representative elementary volume. The mineralogical assemblage included anorthite, anglesite, plumbojarosite and ferrihydrite. Ferrihydrite was determined through interpretation of the mineralogical (XRD) and chemical analyses (ICP-MS, ICP-AES) associated with thermodynamic saturation index (SI) calculations. This assemblage did not reflect all the minerals present in the tailings but rather just the reactive minerals, influencing the mobility of metals. The quantity of each mineral phase was determined based on the total concentrations of Pb, S, Zn, Fe, and Ca measured in the solid phase (Table 2). Pb-bearing phases were represented by anglesite and plumbojarosite. As the mineralogical analyses and SI calculations did not reveal the presence of Zn-bearing phases in the system, the assumption was to incorporate Zn into anglesite and plumbojarosite by substitution with Pb, a process that has already been observed in similar contexts (Dutrizac, 1984; Bril et al., 2008). Further details on the formulation and

<table><caption>**Table 2**

Mineralogical assemblage used in the model at t₃ months</caption><thead><tr><th>Minerals</th><th>Reaction</th><th>Log K<sup>a</sup></th><th>Amount (wt%)</th><th>Surface area (cm<sup>2</sup> g<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Anorthite</td><td>Ca(Al<sub>2</sub>Si<sub>2</sub>)O<sub>8</sub> + 8H<sup>+</sup> ⇌ 2Al<sup>3+</sup> + Ca<sup>2+</sup> + 2H<sub>4</sub>SiO<sub>4</sub></td><td>24.22</td><td>0.1<sup>b</sup></td><td>110</td></tr><tr><td>Anglesite</td><td>Pb<sub>0.98</sub>Zn<sub>0.02</sub>SO<sub>4</sub> ⇌ 0.98Pb<sup>2+</sup> + 0.02Zn + SO<sub>4</sub><sup>2-</sup></td><td>-7.85</td><td>3.25<sup>c</sup></td><td>300</td></tr><tr><td>Plumbojarosite</td><td>Pb<sub>0.455</sub>Zn<sub>0.045</sub>Fe<sub>3</sub>(SO<sub>4</sub>)<sub>2</sub>(OH)<sub>6</sub> + 6H<sup>+</sup> ⇌ 3Fe<sup>3+</sup> + 0.455Pb<sup>2+</sup> + 2SO<sub>4</sub><sup>2-</sup> + 6H<sub>2</sub>O + 0.045Zn<sup>2+</sup></td><td>11.46</td><td>1.4<sup>c</sup></td><td>720</td></tr><tr><td>Ferrihydrite</td><td>Fe(OH)<sub>3</sub> + 3H<sup>+</sup> ⇌ Fe<sup>3+</sup> + 3H<sub>2</sub>O</td><td>3.40</td><td>0.1<sup>c</sup></td><td>4 × 10<sup>5</sup></td></tr></tbody></table>

<sup>a</sup> Values from the Thermoddem database (Blanc et al., 2012).

<sup>b</sup> Fitted values.

<sup>c</sup> Estimated values based on chemical analyses.

implications of this hypothesis were provided in Mertz et al. (2021). The simulation was divided into the same two periods as those used in the pilot experiment. The first period corresponded to the physical stabilization phase from t₃ months to t₀. The duration of this stabilization phase was determined on the basis of time when the residues stopped compacting under the effect of watering. The chemical composition of the initial pore water in the simulations (t₃ months) corresponded to that used in a previous study at the microcosm scale (Mertz et al., 2021) (Table 1). During this first period, no flow was considered. Moreover, and as no monitoring was done, variations of the simulated geochemistry could not be compared to any measurements. The model was therefore not constrained during the first 3 months. At the end of this stabilization phase (t₀), the chemical compositions of the pore water and simulated solid phase were used as input data for simulation of the 6-month leaching phase (t₀ to t₆ months).

### 2.3.2. Boundary conditions

The flow model was designed to account for both the rainfall-like watering and the drainage at the base of the column. A Neumann condition evolving with time was thus applied. The actual flux infiltrating at the top of the column was dependent on: (1) an inflow related to watering, and (2) an outflow due to evaporation and soil water content at time t. If the amount of inflow exceeded the soil infiltration capacity, a maximum 2 cm water layer formed at the top of the column. This water layer infiltrated during the following time step. A drainage condition at the base of the column was set (Dirichlet condition). Moreover, a constant hydraulic charge was considered to keep the last 35 cm of the column in a saturated state. During the first stabilization phase of the experiment (t₃ months to t₀), a water input of 71 L was simulated, corresponding to an input of 9 cm of water. This was evenly distributed over the last 2 months of this period. Concerning the leaching period, three 5.5 L (i.e. 0.7 cm) waterings were applied weekly. Moreover, an evaporation of 0.01 cm day⁻¹ was taken into account in the simulation, this value was adjusted to simulate a water flux balance consistent with that measured.

For the solute transport boundary conditions, a mass flux (Neumann condition) was defined for each solute present in the leaching water at the top of the column. This flux was calculated based on the chemical composition of the Mont Roucous® water (Table 1) and the water flux

<table><caption>**Table 1**

Chemical composition of the watering mineral water (Mont Roucous®) and of the simulated pore water in the tailings at t₃ months</caption><thead><tr><th></th><th>pH</th><th>SO<sub>4</sub><sup>2-</sup><br>(mol L<sup>-1</sup>)</th><th>Pb<br>(mol L<sup>-1</sup>)</th><th>Zn<br>(mol L<sup>-1</sup>)</th><th>Ca<br>(mol L<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Mineral water</td><td>6</td><td>3.2 × 10<sup>-5</sup></td><td>–</td><td>–</td><td>7.2 × 10<sup>-5</sup></td></tr><tr><td>Simulated pore water at t<sub>3months</sub></td><td>6</td><td>4.6 × 10<sup>-4</sup></td><td>4.3 × 10<sup>-5</sup></td><td>7.9 × 10<sup>-6</sup></td><td>4 × 10<sup>-4</sup></td></tr></tbody></table>

<sup>a</sup> Modified concentration to obtain a saturation index close to the thermodynamic equilibrium for the minerals initially present in the tailings.
calculated by the water flow model.

Simulations based on the virtual presence of a non-reactive tracer in the leachate were run to determine the residence time of solutes within the column. In fact, no inert tracer test was performed during the pilot-scale experiment, because it would take too long to obtain a complete breakthrough curve for the bottom of the column. Moreover, the use of chlorine as tracer was overruled due to the high Cl concentrations measured in the pore waters, suggesting the presence of chlorine-bearing phases in the tailings. Indeed, Pb chlorinated phases such as mimetite have already been previously identified in tailings similar to those studied (Courtin-Nomade et al., 2016).

### 2.3.3. Water flow and solute transport

Water flow in the column was considered as being a uniform vertical and variably saturated flow. This phenomenon may be described by the Richards equation:

$$ \frac{\partial \theta}{\partial t} = \frac{\partial}{\partial z} \left[ K \left( \frac{\partial h}{\partial z} + 1 \right) \right] \quad (1) $$

where $\theta$ is the volumetric water content (cm$^3$ cm$^{-3}$), K the unsaturated hydraulic conductivity (cm day$^{-1}$), t the duration (day), h the water pressure head (cm), and z the vertical spatial coordinate. A hydraulic model with a single porosity was used due to the sandy texture of the tailings, and combined with Van Genuchten eqs. (van Genuchten, 1980) to plot the water retention curve:

$$ \theta(h) = \begin{cases} \theta_r + \frac{\theta_s - \theta_r}{[1 + |\alpha h|^n]^m}, & h < 0 \\ \theta_s, & h \ge 0 \end{cases} \quad (2) $$

with $\theta_r$ being the residual water content (cm$^3$ cm$^{-3}$), $\theta_s$ the saturated water content (cm$^3$ cm$^{-3}$), $\alpha$ (cm$^{-1}$), n and m empirical form coefficients, where m = 1–1/n. The unsaturated hydraulic conductivity was defined by its relationship to the hydraulic charge, according to the model of Mualem (1976):

$$ K(h) = \begin{cases} K_s Se^t \left[1 - (1 - Se^{1/m})^m\right]^2, & h < 0 \\ K_s, & h \ge 0 \end{cases} \quad (3) $$

where $K_s$ is the saturated hydraulic conductivity (cm day$^{-1}$), l a soil tortuosity parameter and $Se$ the effective saturation where $Se = (\theta - \theta_r)/ (\theta_s - \theta_r)$.

The hydraulic model parameter values were defined on the basis of the particle size properties of the tailings determined by wet processing sieving. Carsel and Parrish (1988) established a relationship between soil texture and hydrodynamic parameters, thereby offering a first approximation of the system properties. The tailings were composed of 93 wt% grains ranging from 2 mm to 20 µm in size. The hydrodynamic parameter values considered in the simulations are reported in Table 3. The longitudinal dispersivity and the water content at saturation were then adjusted to respect the water balance and the water contents measured at the different column depths. The transport of solutes in a variably saturated medium was described by an advection-dispersion equation:

$$ \frac{\partial \theta C}{\partial t} = \frac{\partial}{\partial z} \left[ \theta \cdot D \frac{\partial C}{\partial z} \right] - q \frac{\partial C}{\partial z} \quad (4) $$

<table><caption>**Table 3**

Parameters used in the simulation for water flow and solute transport.</caption><thead><tr><th>Parameters</th><th>θ<sub>r</sub></th><th>θ<sub>s</sub><sup>*</sup></th><th>α<sup>*</sup> (cm<sup>-1</sup>)</th><th>n</th><th>K<sub>s</sub> (cm day<sup>-1</sup>)</th><th>D<sub>t</sub><sup>*</sup> (cm)</th></tr></thead><tbody><tr><td>Values</td><td>0.045</td><td>0.38</td><td>0.1</td><td>2.68</td><td>712.8</td><td>8</td></tr></tbody></table>

with C being the concentration of solutes in the liquid (mol cm$^{-3}$), D the dispersion coefficient (cm$^2$ day$^{-1}$) and q the average velocity of the fluid in the pores (cm day$^{-1}$). The dispersion coefficient D was equal to DLq with DL being the longitudinal dispersivity (cm). The transport parameter values were determined to obtain the best representation of the measured data used for model calibration (Table 3). Diffusion processes were not taken into account in the simulations, as they are considered to be negligible with respect to the resulting high hydraulic conductivity of the sandy tailings texture.

### 2.3.4. Geochemical processes

The geochemical model considered in the simulations involved two main processes, i.e. surface complexation reactions on iron oxides and thermokinetic dissolution/precipitation reactions of primary and secondary mineral phases. The complexation reactions were simulated using the model of Dzombak and Morel (1990), a formalism developed to simulate the surface complexation of cations on the iron oxide surfaces. The sorption of three elements was considered in the simulations: Pb and Zn as major contaminants in our system, and Ca as the major cation that could have a competitive role with respect to the metallic cations.

These complexation reactions were described using the following generic formalism:

$$ S^{\alpha s} (A_i^{\alpha i})_{v_i} + v_j A_j^{\alpha j} = S^{\alpha i} (A_i^{\alpha i})_{v_j} + v_i A_i^{\alpha i} \quad (5) $$

with S being the surface chemical species, z the charge of the chemical species, $A_{i,j}$ the chemical species of the adsorbate, and $v_{i,j}$ the stoichiometric coefficients. A thorough description of the complexation reactions is presented in the Appendix (Table A.1/2).

The dissolution and precipitation of the mineral phases in the system were described by a kinetic law derived from the transition state theory (Lasaga, 1981, 1984):

$$ \frac{dm}{dt} = -A \cdot \sum_j \left[ k_j e^{\frac{-E_j}{R}} \left( \frac{1}{T} - \frac{1}{298.15} \right) f_j(a_{ij}) g_j(\Delta G_r) \right] \quad (6) $$

where $dm/dt$ is the reaction rate (mol s$^{-1}$), A the surface area of the mineral (m$^2$), k the rate constant (mol m$^{-2}$ s$^{-1}$), E the activation energy (J mol$^{-1}$), R the gas constant (J mol$^{-1}$ K$^{-1}$) and T the temperature (K). The unitless term $f_j(a_{ij})$ represents the activity of an $i_{th}$ aqueous chemical species participating in the $j_{th}$ process. The unitless term $g_j(\Delta G_r)$ reflects the distance from thermodynamic equilibrium. The kinetic parameters are described in Table 4.

### 2.3.5. Model performance criteria

To assess the model prediction accuracy, the simulated values were compared to the measured ones using the normalized root mean square error (NRMSE), with the equation:

$$ \text{NRMSE} = \frac{\sqrt{\frac{1}{n} \sum_{i=1}^{n} (p_i - m_i)^2}}{y} \quad (7) $$

where n is the number of observations, $p_i$ the predicted values, $m_i$ the measured values, and y the mean value of the measured data. A low NRMSE value implies a good fit of the model to the measured data.

## 3. Results

### 3.1.1. Flow and mass balances

The simulated cumulative column outflows were compared to those

<footnote>* Adjusted values.</footnote>
<table><caption>**Table 4**
Kinetic parameters used to describe dissolution and precipitation of mineral phases in the simulations.</caption><thead><tr><th rowspan="2">Minerals</th><th colspan="6">Kinetic parameters</th><th colspan="3">Neutral mechanism</th><th colspan="3">Base mechanism</th></tr><tr><th colspan="3">Acid mechanism</th><th colspan="3">Neutral mechanism</th><th rowspan="2">k<sub>diss</sub><sup>a</sup></th><th rowspan="2">k<sub>pre</sub><sup>b</sup></th><th rowspan="2">E<sup>c</sup></th><th rowspan="2">k<sub>diss</sub><sup>a</sup></th><th rowspan="2">k<sub>pre</sub><sup>b</sup></th><th rowspan="2">E<sup>c</sup></th><th rowspan="2">n<sup>g</sup></th></tr></thead><tbody><tr><td>Anorthite<sup>a</sup></td><td>2.73 × 10<sup>-3</sup></td><td>2.73 × 10<sup>-4</sup></td><td>16.6</td><td>1.41</td><td>6.54 × 10<sup>-9</sup></td><td>6.54 × 10<sup>-10</sup></td><td>17.8</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Anglesite</td><td>8.6 × 10<sup>-7</sup><sup>*</sup></td><td>6.9 × 10<sup>-7</sup><sup>*</sup></td><td>31.3<sup>d</sup></td><td>0.3<sup>d</sup></td><td>6.9 × 10<sup>-8</sup><sup>*</sup></td><td>5.2 × 10<sup>-8</sup><sup>*</sup></td><td>31.3<sup>d</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Plumbojarosite<sup>w</sup></td><td>2.6 × 10<sup>-4</sup></td><td>2.6 × 10<sup>-5</sup></td><td>79</td><td>0.89</td><td>1.7 × 10<sup>-5</sup></td><td>1.7 × 10<sup>-6</sup></td><td>79</td><td>9.5 × 10<sup>-10</sup></td><td>9.5 × 10<sup>-12</sup></td><td>79</td><td>0.39</td><td>—</td></tr><tr><td>Ferrihydrite<sup>*</sup></td><td>—</td><td>—</td><td>—</td><td>—</td><td>3.4 × 10<sup>-9</sup></td><td>3.4 × 10<sup>-10</sup></td><td>86.5</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

<sup>a</sup> Dissolution rate constant k<sub>diss</sub> in mol cm<sup>-2</sup> day<sup>-1</sup>.

<sup>b</sup> Precipitation rate constant k<sub>pre</sub> in mol cm<sup>-2</sup> day<sup>-1</sup>.

<sup>c</sup> Arrhenius activation energy E in kJ mol<sup>-1</sup>.

<sup>d</sup> Values from Dove and Czank, 1995

<sup>e</sup> Values from Palandri and Kharaka, 2004

<sup>f</sup> Reaction order n with respect to H<sup>+</sup>.

<sup>g</sup> Reaction order n with respect to OH<sup>-</sup>.

<sup>*</sup> Adjusted values.

measured (Fig. 2a). During the experiment, 401 L of Mont Roucous® water were added to the top of the column and 388 L were recovered at the column outlet, which was equivalent to a water layer of approximately 49.4 cm. The simulation well reproduced the experimental data: 51.1 cm of water were added over this 6-month period, 1.7 cm were evaporated, so there was 49.4 cm outflow at the outlet. The cumulative evolution of the simulated outflows was linear while the measured values exhibited slight temporal variations, however the model accurately predicted these data (NRMSE = 6.8 %).

The measured and simulated water contents at the different depths (D1, D2, D3) are shown in Fig. 2. In the unsaturated zone (D1), the measured water content averaged 28.1 % ± 1.1 %, while the simulated water content averaged 27.83 % and ranged from 27.75 % to 28.15 % depending on the watering cycles. In the saturated zone, the measured water contents were slightly different depending on the depth. The

average water content was 38.9 ± 0.6 % at D2, and 37.1 ± 0.6 % at D3. By imposing a hydraulic charge in the model for the saturated zone, the water content was constant at 38 %. This value was close to the measurements recorded during the leaching experiment. Overall, the simulated data closely matched the measured data (NRMSE = 6.1 %, 2.8 % and 2.9 % for D1, D2 and D3, respectively).

The solute transit time through the column during the leaching experiment was calculated on the basis of the results of the fate of a virtual non-reactive tracer from 5 days to 6 months (Fig. 3). In this simulation, 20 days were thus required for the first drops of the non-reactive tracer which was continuously added in the watering solution to reach the bottom of the column, and 6 months were necessary to obtain the same tracer concentration in the pore water at the base of the column and in the watering solution. A slight overconcentration was also observed in the first 2 cm of the column, reaching a maximum C/C₀

<figure><img src="image_7.png" /><figcaption>**Fig. 2.** Temporal patterns of cumulative outflows (a), volumetric water content at D1 (b), D2 (c) and D3 depths (d). Red lines represent the measured data while blue lines represent the simulated values.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 3.** Simulation of variations in the C/C₀ ratio with depth of the virtual non-reactive tracer (Trc) from 5 days to 6 months, where C is the concentration at a given time and C₀ its concentration initially present in the leaching water.</figcaption></figure>

ratio of 1.04. This could be attributed to the evaporation process, which slightly reduced the volume of water present in the unsaturated zone, thereby increasing the chemical element concentrations.

## 3.2. Variations in pore water geochemical properties and mineralogical assemblage

The simulated pH changes at the different depths was reported with the measured values in Fig. 4.a/c/e. Very close agreement between the measured and simulated data was observed (NRMSE = 10.4 %, 11.1 % and 11 % for D1, D2 and D3, respectively). The measured pH values remained constant at ~4 over the 6 months, with an average of 4.08 ± 0.45 for D2 and 3.98 ± 0.36 for D3. The measured pH values were slightly higher for D1 (4.32 ± 0.45), and > 5 in the first 40 days. The simulated pH values were also constant at ~4.

The modeled variations in Pb concentration in solution gave results very close to the measured values (Fig. 4.b/d/f) (NRMSE = 12.9 %, 23.2 % and 24.8 % for D1, D2, and D3, respectively). Dynamics similar to those of the measured Pb concentrations were highlighted at each depth, i.e. a slight temporal decrease from ~8 × 10⁻⁵ mol L⁻¹ to 6 × 10⁻⁵ mol L⁻¹. For the modeled values, an 11 % concentration decrease was observed in the first hours of the simulation, i.e. from 6.4 × 10⁻⁵ mol L⁻¹ to 5.7 × 10⁻⁵ mol L⁻¹. This last value continued to slightly but steadily decrease for different durations depending on the depth. The lowest Pb concentration was reached after 2.5 days for D1 with 5.6 × 10⁻⁵ mol L⁻¹, after 15 days for D2 with 5.4 × 10⁻⁵ mol L⁻¹, and after 35 days for D3 with 5.1 × 10⁻⁵ mol L⁻¹. Then the simulated values increased until reaching a plateau whose values decreased with depth: 6.8 × 10⁻⁵ mol L⁻¹, 6 × 10⁻⁵ mol L⁻¹ and 5.4 × 10⁻⁵ mol L⁻¹ for D1, D2, and D3, respectively.

The simulated variations in Zn concentration at different depths was compared to their measured values (Fig. 5.a/b/c). The predicted Zn concentration values were consistent with the measured values (NRMSE = 15.7 %, 7.6 % and, 28.3 % for D1, D2 and D3, respectively). The measured Zn concentrations at different depths remained constant at ~1 × 10⁻⁵ mol L⁻¹. The initially simulated values at 1.1 × 10⁻⁵ mol L⁻¹ decreased or increased until reaching a plateau whose value slightly increased with depth, reaching 0.84 × 10⁻⁵ mol L⁻¹, 1 × 10⁻⁵ mol L⁻¹, and 1.2 × 10⁻⁵ mol L⁻¹ at D1, D2 and D3, respectively. The model was able to quite accurately reproduce the measured SO₄²⁻ concentrations (NRMSE = 39.7 %, 45.5 % for D1, and D2, respectively), except for D3 which had a high NRMSE of 81.2 % (Fig. 5.d/e/f). The measured SO₄²⁻ concentrations did not show clearcut dynamics, with a mean concentration of 3.8 ± 1.2 × 10⁻⁴ mol L⁻¹, indicating substantial dispersion of the values, mainly during the first 40 days.

As for Zn, the predicted SO₄²⁻ concentrations reached a plateau, whose value increases with depth, with 4 × 10⁻⁴ mol L⁻¹, 4.8 × 10⁻⁴ mol L⁻¹, and 5.7 × 10⁻⁴ mol L⁻¹ for D1, D2 and D3, respectively. The greater discrepancy noted between the simulated and measured values for sulfates could have been due to the influence of surface complexation reactions that were not taken into account. Finally, regarding Ca concentration pattern, the simulation correctly predicted the measured data (NRMSE = 25.9 %, 37.3 % and 42 % for D1, D2, and D3, respectively) (Fig. 5.g/h/i). The measured Ca values did not show clearcut dynamics due to the scant number of available data and their dispersion. Nevertheless, the average values at each depth were ~ 1.45 × 10⁻⁴ mol L⁻¹, 1.4 × 10⁻⁴ mol L⁻¹ and 1.22 × 10⁻⁴ mol L⁻¹. The simulated Ca variation patterns were similar to those of Zn and SO₄²⁻ with a plateau and slightly increasing with depth, with 1.2 × 10⁻⁴ mol L⁻¹, 1.43 × 10⁻⁴ mol L⁻¹, and 1.65 × 10⁻⁴ mol L⁻¹ for D1, D2 and D3, respectively. All measured
<figure><img src="image_9.png" /><figcaption>**Fig. 4.** Temporal variations in pH at D1 (a), D2 (c) and D3 (e), and the Pb concentration in solution at D1 (b), D2 (d) and D3 (f) depths. Red dots correspond to the experimental data while the lines refer to the simulated values. For the inset plots, the x-axis is expressed in days.</figcaption></figure>

Fe concentrations in pore water at D1, D2 and D3 levels were below the detection limit of 0.1 mg L⁻¹ (i.e. [Fe] < 1.8.10⁻⁶ mol L⁻¹). These very low Fe concentrations could be explained by the oxidizing conditions, which promoted Fe precipitation.

The variation patterns of the most reactive mineral phases of the system (plumbojarosite, anglesite, ferrihydrite and anorthite) in the simulation are presented in Fig. 6. Plumbojarosite kept dissolving during the simulation, but at different rates depending on the depth. This dissolution was quite slow and constant throughout the column, resulting in the dissolution after 6 months of 0.5 % of the quantity of plumbojarosite initially present (1.4 wt%). However, this dissolution increased in the first 10 cm of the column, reaching 4.7 % maximum

depletion after 6 months of simulation. The simulated ferrihydrite variation pattern was similar but in the opposite way. Indeed, only a precipitation reaction affected this mineral within the column at different rates depending on the depth. After 6 months, the quantity of ferrihydrite increased by 3.5 % between 10 and 70 cm depth, with a maximum 25 % increase at the top of the column (i.e. from 0.11 to 0.14 wt%). The amount of precipitated ferrihydrite would thus generate more available sorption sites for Pb and Zn retention. The anglesite variation pattern revealed two opposite behaviors according to the depth. In the simulation, anglesite was partially dissolved from 0 to 2.5 cm depth, while from 2.5 to 70 cm this mineral slightly precipitated; at 2.5 cm the quantity of anglesite remained temporally constant, with a value of 3.25
<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Temporal variations in Zn concentrations at D1 (a), D2 (b) and D3 (c), SO₃²⁻ at D1 (d), D2 (e) and D3 (f), and Ca at D1 (g), D2 (h) and D3 (i) depths. Red dots correspond to the experimental data while the lines represent simulated values.</figcaption></figure>

wt%. The dissolution rate was greater than the precipitation rate in the top 2.5 cm of the column. Indeed, after 6 months a maximum of 2.7 % of anglesite should be dissolved at the surface, while its more uniform precipitation within the column resulted in an increase of only 0.07 % in the anglesite quantity. Finally, the simulated anorthite variation pattern resulted in low uniform dissolution within the column, i.e. 1.4 % of the amount of anorthite initially present (0.1 wt%).

# 4. Discussion

## 4.1. Control of geochemical parameters

This study revealed several mechanisms controlling the Pb/Zn dynamics in the mesocosm. Plumbojarosite dissolution led to the release of Fe³⁺ ions in solution which were then incorporated by ferrihydrite during its precipitation, with the latter involving proton release. These concomitant dissolution and precipitation reactions highlighted by RTM resulted in maintenance of an acidic environment in the tailings. This was due to the regular watering, which resulted in a thermodynamic imbalance between the pore water and solid phase. This relationship between plumbojarosite and ferrihydrite was perfectly reflected by the RTM, with these two minerals showing opposite variation patterns. Ferrihydrite precipitation in response to jarosite dissolution has already been described in previous studies (Welch et al., 2008; Elwood Madden et al., 2012) and the resulting acidification of ~ pH 4 was observed (Smith et al., 2006). Anorthite also slightly affected the pH, as its dissolution reaction consumed protons. This may explain the slight

increase in pH with depth in the simulation.

The Pb concentration in solution was mainly controlled by the dissolution and precipitation of the Pb-bearing phases, i.e. anglesite and plumbojarosite. Dilution occurred during watering, thereby decreasing the Pb concentration in solution, which explained the marked dissolution of plumbojarosite at the top of the column due to the remoteness from thermodynamic equilibrium. This effect of higher mineral dissolution at the surface was previously observed by Pabst et al. (2017) when designing an RTM for a column filled with tailings. When the leaching water interacted with tailings, it became enriched in Pb, which slowed down the dissolution of plumbojarosite and anglesite. From 2.5 cm depth of, the Pb concentration in solution was high enough to be oversaturated with respect to anglesite, which began to precipitate. Several studies have highlighted these dynamics, i.e. Pb from plumbojarosite can be incorporated into anglesite (Smith et al., 2006; Lu and Wang, 2012). The initial 11 % drop that occurred in the first hours of the simulation could be directly related to anglesite precipitation, as the pore water was slightly oversaturated at the beginning of the simulation. The Pb concentration pattern then involved a smaller decrease before increasing and reaching an asymptote. The time lag required to reach this asymptote was due to the time required for the leaching water to reach D1, D2 then D3, thus delaying the establishment of a metastable state in the bottom of the column. The plateau level at a certain depth was affected by the behavior of anglesite and plumbojarosite, while the difference in Pb concentration between the different depths could be explained by the fact that the anglesite precipitation was greater than the plumbojarosite dissolution rate. This effect was noticeable at >2.5
<figure><img src="image_11.png" /><figcaption>**Fig. 6.** Depth-dependent variations in the quantity of mineral phases present in the tailings at 0 days, 20 days, 3 months and 6 months: (a) plumbojarosite, (b) ferrihydrite, (c) anglesite, and (d) anorthite.</figcaption></figure>

cm depth, where Pb tended to be incorporated into anglesite while its concentration in solution decreased. This effect was notably accentuated by the presence of a higher proportion of Pb in anglesite than in plumbojarosite. The solution pH also impacted the Pb concentration in solution, thus altering the sorption capacity of Pb by iron oxides. Once the column reached a metastable state, a higher pH would have led to a decrease in the Pb concentration in solution. The pH was controlled especially by the dissolution of anorthite. However, this effect remained minor and the Pb continued to be mainly controlled by dissolution of the Pb bearing phases. As the reaction rate of these processes was relatively fast, this decrease in Pb was mainly governed by thermodynamic phenomena.

The behavior of Zn with respect to pH differed from that of Pb. Indeed, Zn mobility was only controlled by its carrier phases, i.e. anglesite and plumbojarosite, whereas the amount of Zn adsorbed on iron oxide was negligible. In contrast to Pb, the amount of Zn present in

plumbojarosite was higher than in anglesite, so the amount of Zn released from plumbojarosite was greater than that incorporated by anglesite, resulting in an increase in the Zn concentration in solution with depth. This process was similar for $SO_4^{2-}$, whose concentration also increased with depth. Finally, the Ca variation pattern was only controlled by anorthite dissolution. Thermodynamic equilibrium could not be reached due to the slow dissolution kinetics of this mineral, so the solution became enriched in Ca throughout the column.

## 4.2. Scale change and process validity

A key issue in RTM is the consistency of processes across different spatial scales (Steefel et al., 2005, 2015). In this study, there was continuity in the results between the simulations carried out at a centimetric scale in the microcosms experiments (Mertz et al., 2021) and those carried out at a metric scale in the pilot experiment (Thouin et al., 2022).
The parameters describing the geochemical processes were identically transcribed between the two models. The only exception concerned the anorthite surface area, which was slightly modified compared to the previous model (0.011 instead of 0.014 m² g⁻¹) in order to obtain the best possible fit between the data and the model. This minor variation could be explained by a difference in solid phase compaction between the two experiments. As anorthite was the only mineral whose dissolution and precipitation reactions were highly controlled by its kinetics, the pore water was always undersaturated with respect to this mineral in the experimental conditions. A slight modification in its available surface could thus potentially influence the system.

The comparison between microcosm and column experiment simulations showed that the parameters used to describe the geochemical reactions had similar values. This suggests that the geochemical reaction intensity was comparable at both scales. The microcosm leaching experiments conducted with only 150 g of tailings (Mertz et al., 2021) thus generated a sufficient data-set to develop a biogeochemical model that helped determine the nature and intensity of the geochemical reactions controlling the fate of Pb and Zn in the Pontgibaud tailings. This hypothesis, based on a comparison of the simulations, suggests that the microcosm volume was greater than the representative elemental volume (REV) for the Pontgibaud tailings, thus curbing issues related to tailing heterogeneity. Nevertheless, some questions remain, notably on the REV size in other types of tailings. Indeed, the tailings were represented by a single material with a relatively simple mineralogy (four reactive mineral phases), mainly controlled by thermodynamic equilibrium. Tailings with a more complex mineralogical assemblage and/or a greater kinetic limitation in its reactivity could imply more complexity in the implementation of a scale change. However, this modeling methodology used to gain insight into the behavior of metallic contaminants in the Pontgibaud tailings could be applied to other old mining sites with a reasonable number of reactive mineral phases containing metallic contaminants. In a previous study carried out at the centimetric scale, transport was simulated via a simplified transport model based on successive mixtures of pore water and leaching water (Mertz et al., 2021). The experimental setup in microcosms did not enable advective flow measurement, thus hampering description of transport processes in the model. Monitoring of the instrumented column led to an accurate description of flow and solute transport. The parameters describing these processes were mainly determined on the basis of the study of Carsel and Parrish (1988), thereby enabling comparison of the hydrodynamic properties of the tailings to their granulometric properties. Some parameters were adjusted to be able to accurately describe the water and mass balances in the column, notably the saturation water content θₛ was adjusted to 0.38. This value was nevertheless in the range of those in the classification of Carsel and Parrish (1988), i.e. 0.43 ± 0.06 for materials with a sandy texture such as the studied tailings. Moreover, similar values for this parameter have been used in previous studies aimed at modeling water flow in sandy matrices (Looms et al., 2008; Zheng et al., 2017; Wang et al., 2018). The empirical parameter α was adjusted to 0.1 cm⁻¹, i.e. slightly lower than the value proposed by Carsel and Parrish (0.145 cm⁻¹ ± 0.029). Nevertheless, this was higher than the values used in studies on a similar medium (0.03 to 0.08 cm⁻¹) (Zheng et al., 2017; Wang et al., 2018). Finally, as no precise dispersivity could be determined due to the absence of a tracer test, a value of 8 cm was used in the simulations. This dispersivity value was in the same range as that used by Looms et al. (2008) in a sandy soil (5 cm). The present simulations revealed that different geochemical parameters (pH, Pb and Zn concentrations) eventually reached an asymptote. The dispersivity value only influenced the time needed to reach this plateau but not its magnitude. Consequently, the value configured in our simulation was a trade-off between a good representation of the data and a value in line with those used in previous studies for sandy matrices.

### 4.3. RTM contribution to mine site management

Over the past 30 years, RTM have emerged as powerful tools for exploring processes involved in the critical zone (Li et al., 2017), as reflected by the marked growth in the number of RTM studies since the 1990s (Deng et al., 2021). Authors of many studies focused on recent mining sites have developed RTM to understand and predict the effects of acid mine drainage on the mobility of metal contaminants (Walter et al., 1994; Wunderly et al., 1996; Amos et al., 2004; Brookfield et al., 2006; Mayer et al., 2015). The ability of RTM to relatively accurately describe the properties of a system offers promising prospects in terms of management of these mining sites. Indeed, RTM can be used as environmental engineering tools to assess the effectiveness of several remediation scenarios (Runkel and Kimball, 2002). However, few similar approaches have been undertaken in studies on former tailings with highly weathered primary minerals. The results of this study allowed assessment of the effects of a very fundamental phenomenon in terms of contaminant migration, i.e. the role of natural attenuation (Steefel et al., 2005). The decreased Pb concentration in solution with depth seemed at first sight to suggest a natural attenuation effect. Yet the intensity of this phenomenon remained relatively low and this attenuation of Pb was offset by an increase in the Zn concentration in deep horizons. Moreover, the system dynamics were relatively stable and this release of contaminants into the environment could likely extend over long periods. Indeed, although plumbojarosite dissolution was heterogeneous with depth, it is generally slow (~ 1 %/year). The results could hence be extrapolated to estimate the long-term fate of Pb and Zn fluxes from tailings. When taking the slowest dissolution that occurred throughout most of the column into account, a decrease of 0.015 wt %/year in the amount of plumbojarosite was noted. This would result in total dissolution of this mineral after 93 years, but this estimation was based on extrapolation of the simulations performed over the first 6 months. This total plumbojarosite dissolution would result in a pH increase due to the absence of Fe release, but not necessarily in a decrease in Pb or Zn concentrations. Indeed, continuous dissolution of plumbojarosite would lead to anglesite precipitation, thus increasing its already substantial total content in the tailings. Once plumbojarosite was completely dissolved, it was calculated that the total amount of anglesite would reach 3.7 wt%. In the absence of Pb release from plumbojarosite, anglesite should in turn dissolve and release significant amounts of Pb and Zn over long periods. In addition, a less acidic solution would have a relatively limited effect on metal mobility due to the small number of sorption sites available via Fe oxide. Hence, although it would be hard to accurately estimate, the release of metal contaminants from mine tailings appears to be sustainable over time. Consequently, natural attenuation of these contaminants cannot be relied on in remediation efforts.

### 5. Conclusion

Data recorded in a leaching experiment performed in a metric column filled with former mine tailings were used to design a reactive transport model that combines mechanistic descriptions of the different processes controlling the system reactivity, and a rigorous description of the flow and solute transport dynamics. The results suggest that the mobility of the most abundant metal contaminants (Pb, Zn) was mainly controlled by the dissolution and precipitation of their bearing phases, e. g. anglesite and plumbojarosite. The RTM designed in this study highlighted that the dissolution/precipitation dynamics between anglesite and plumbojarosite fluctuated with depth, resulting in variable Pb and Zn concentrations. It was shown that the pH was maintained at an acidic level (~ 4) due to the precipitation of ferrihydrite induced by Fe release via plumbojarosite dissolution. These different dynamics documented at the metric scale were similar to those observed at the centimetric scale. It could be assumed that the volume of leached tailings considered in the microcosms represents a volume greater than the REV, limiting the influence of the heterogeneities of the system. Finally, the variations in the
different mineral phases within the column over the 6 months of leaching allowed estimation of the long-term fate of the tailings. According to this RTM, plumbojarosite dissolution was likely to continue for about 93 years, while maintaining acidic conditions and a high Pb concentration. Thereafter, anglesite would in turn start to dissolve, also releasing large amounts of Pb into the solution. These estimations highlight the irrelevance of relying on natural attenuation in old mine tailings rehabilitation projects, even in the long term. One possible way to effectively immobilize these metals would involve the addition of inorganic and/or organic amendments. RTMs is thus a valuable tool for the management of these mining sites, at it enables assessment of temporal variations in the system reactivity.

### CRediT authorship contribution statement

**Samuel Mertz:** Writing – original draft, Formal analysis, Conceptualization. **Nicolas Devau:** Writing – original draft, Visualization, Validation, Supervision, Methodology. **Hugues Thouin:** Validation, Methodology, Conceptualization. **Fabienne Battaglia-Brunet:** Validation, Funding acquisition. **Marie-Paule Norini:** Visualization, Formal analysis. **Marc Crampon:** Visualization, Formal analysis. **Lydie Le Forestier:** Writing – original draft, Supervision, Funding acquisition.

### Declaration of competing interest

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

### Acknowledgements

This research was performed within the framework of the Phytoselect project funded by the Région Centre - Val de Loire (contract N°2016-00108485), and by the Labex Voltaire (ANR-10-LABX-100-01). S. Mertz benefited from a Ph.D. grant provided by the French Geological Survey (BRGM) and by the Région Centre - Val de Loire. The authors gratefully acknowledge the financial support provided to the PIVOTS project by the Région Centre - Val de Loire: ARD 2020 program, CPER 2015-2020, and the European Union, which invests in the Centre-Val de Loire via the European Regional Development Fund. Finally, the authors thank David Manley for English editing.

### Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.scitotenv.2024.178248.

### Data availability

Data will be made available on request.

### References

Acero, P., Ayora, C., Carrera, J., Saaltink, M.W., Olivella, S., 2009. Multiphase flow and reactive transport model in vadose tailings. Appl. Geochem. 24, 1238-1250. https://doi.org/10.1016/j.apgeochem.2009.03.008.

Akcil, A., Koldas, S., 2006. Acid mine drainage (AMD): causes, treatment and case studies. J. Clean. Prod. 14, 1139-1145. https://doi.org/10.1016/j.jclepro.2004.09.006.

Akhavan, A., Golchin, A., 2021. Estimation of arsenic leaching from Zn-Pb mine tailings under environmental conditions. J. Clean. Prod. 295, 126477. https://doi.org/10.1016/j.jclepro.2021.126477.

Amos, R.T., Mayer, K.U., Blowes, D.W., Ptacek, C.J., 2004. Reactive transport modeling of column experiments for the remediation of acid mine drainage. Environ. Sci. Technol. 38, 3131-3138. https://doi.org/10.1021/es0349608.

Anju, M., Banerjee, D.K., 2010. Comparison of two sequential extraction procedures for heavy metal partitioning in mine tailings. Chemosphere 78, 1393-1402. https://doi.org/10.1016/j.chemosphere.2009.12.064.

Atanacković, N., Dragićić, V., Stojković, J., Papić, P., Živanović, V., 2013. Hydrochemical characteristics of mine waters from abandoned mining sites in Serbia and their

impact on surface water quality. Environ. Sci. Pollut. Res. 20, 7615-7626. https://doi.org/10.1007/s11356-013-1959-4.

Bernhard, D., Reilly, I.I.J.F., 2020. Mineral Commodity Summaries 2020. US Geological Survey.

Blanc, Ph., Lassin, A., Piantone, P., Azaroual, M., Jacquemet, N., Fabbri, A., Gaucher, E. C., 2012. Thermoderm: a geochemical database focused on low temperature water/rock interactions and waste materials. Appl. Geochem. 27, 2107-2116. https://doi.org/10.1016/j.apgeochem.2012.06.002.

Blowes, D.W., Jambor, J.L., 1990. The pore-water geochemistry and the mineralogy of the vadose zone of sulfide tailings, Waite amulet, Quebec, Canada. Appl. Geochem. 5, 327-346. https://doi.org/10.1016/0883-2927(90)90008-S.

Bril, H., Zainoun, K., Puziewicz, J., Courtin-Nomade, A., Vanaecker, M., Bollinger, J.-C., 2008. Secondary phases from the alteration of a pile of zinc-smelting slag as indicators of environmental conditions: an example from Swietochowice, upper Silesia, Poland. Can. Mineral. 46, 1235-1248. https://doi.org/10.3749/canmin.46.5.1235.

Brookfield, A.E., Blowes, D.W., Mayer, K.U., 2006. Integration of field measurements and reactive transport modelling to evaluate contaminant transport at a sulfide mine tailings impoundment. J. Contam. Hydrol. 88, 1-22. https://doi.org/10.1016/j.jconhyd.2006.05.007.

Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24, 755-769. https://doi.org/10.1029/WR024i005p00755.

Cerdà, A., 1998. Soil aggregate stability under different Mediterranean vegetation types. CATENA 32, 73-86. https://doi.org/10.1016/S0341-8162(98)00041-1.

Chaney, K., Swift, R.S., 1984. The influence of organic matter on aggregate stability in some British soils. J. Soil Sci. 35, 223-230. https://doi.org/10.1111/j.1365-2389.1984.tb00278.x.

Courtin-Nomade, A., Waltzing, T., Evrard, C., Soubrand, M., Lenain, J.-F., Ducloux, E., Ghorbel, S., Grosbois, C., Bril, H., 2016. Arsenic and lead mobility: from tailing materials to the aqueous compartment. Appl. Geochem. 64, 10-21. https://doi.org/10.1016/j.apgeochem.2015.11.002.

Csavina, J., Field, J., Taylor, M.P., Gao, S., Landázuri, A., Betterton, E.A., Sáez, A.E., 2012. A review on the importance of metals and metalloids in atmospheric dust and aerosol from mining operations. Sci. Total Environ. 433, 58-73. https://doi.org/10.1016/j.scitotenv.2012.06.013.

Deng, H., Navarre-Sitchler, A., Heil, E., Peters, C., 2021. Addressing water and energy challenges with reactive transport modeling. Environ. Eng. Sci. 38, 109-114. https://doi.org/10.1089/ees.2021.00009.

Dove, P.M., Czank, C.A., 1995. Crystal chemical controls on the dissolution kinetics of the isostructural sulfates: Celestite, anglesite, and barite. Geochim. Cosmochim. Acta 59, 1907-1915. https://doi.org/10.1016/0016-7037(95)00116-6.

Dutrizac, J.E., 1984. The behavior of impurities during Jarosite precipitation. In: Bautista, R.G. (Ed.), Hydrometallurgical Process Fundamentals. Springer US, Boston, MA, pp. 125-169. https://doi.org/10.1007/978-1-4899-2274-8_6.

Dzombak, D.A., Morel, F.M.M., 1990. Surface Complexation Modeling: Hydrous Ferric Oxide. Wiley, New York.

Elwood Madden, M.E., Madden, A.S., Rimstidt, J.D., Zairah, S., Kendall, M.R., Miller, M. A., 2012. Jarosite dissolution rates and nanoscale mineralogy. Geochim. Cosmochim. Acta 91, 306-321. https://doi.org/10.1016/j.gca.2012.05.001.

Evangelou, V.P., Zhang, Y.L., 1995. A review: pyrite oxidation mechanisms and acid mine drainage prevention. Crit. Rev. Environ. Sci. Technol. 25, 141-199. https://doi.org/10.1080/10643389509388477.

Gemici, Ű., 2008. Evaluation of the water quality related to the acid mine drainage of an abandoned mercury mine (Alašehir, Turkey). Environ. Monit. Assess. 147, 93-106. https://doi.org/10.1007/s10661-007-0101-9.

Harris, J.R., 1988. The British iron industry, 1700-1850, Studies in economic and social history. Macmillan Education, Basingstoke.

Hayes, S.M., Webb, S.M., Bargar, J.R., O'Day, P.A., Maier, R.M., Chorover, J., 2012. Geochemical weathering increases Lead bioaccessibility in semi-arid mine tailings. Environ. Sci. Technol. 46, 5834-5841. https://doi.org/10.1021/es300603s.

Hiller, E., Tóth, R., Kučerová, G., Jurković, L., Sottník, P., Lalinská-Voleková, B., Vozár, J., 2016. Geochemistry of mine tailings from processing of siderite-cu ores and mobility of selected metals and metalloids evaluated by a pot leaching experiment at the Slovinky impoundment, eastern Slovakia. Mine Water Environ. 35, 447-461. https://doi.org/10.1007/s10230-016-0388-2.

Jurjovec, J., Blowes, D.W., Ptacek, C.J., Mayer, K.U., 2004. Multicomponent reactive transport modeling of acid neutralization reactions in mine tailings: modeling of acid neutralization reactions. Water Resour. Res. 40, W11202. https://doi.org/10.1029/2003WR002233.

Jurjovec, J., Ptacek, C.J., Blowes, D.W., 2002. Acid neutralization mechanisms and metal release in mine tailings: a laboratory column experiment. Geochim. Cosmochim. Acta 66, 1511-1523. https://doi.org/10.1016/S0016-7037(01)00874-2.

Lasaga, A.C., 1981. Transition state theory. Rev. Mineral. Geochem. 8, 135-168.

Lasaga, A.C., 1984. Chemical kinetics of water-rock interactions. J. Geophys. Res. 89, 4009-4025. https://doi.org/10.1029/JB089iB06p04009.

Lee, P.-K., Kang, M.-J., Jo, H.Y., Choi, S.-H., 2012. Sequential extraction and leaching characteristics of heavy metals in abandoned tungsten mine tailings sediments. Environ. Earth Sci. 66, 1909-1923. https://doi.org/10.1016/j.s12665-011-1415-2.

Li, L., Maher, K., Navarre-Sitchler, A., Druhan, J., Meile, C., Lawrence, C., Moore, J., Perdrial, J., Sullivan, P., Thompson, A., Jin, L., Bolton, E.W., Brantley, S.L., Dietrich, W.E., Mayer, K.U., Steefel, C.I., Valocchi, A., Zachara, J., Kocar, B., Mcintosh, J., Tutolo, B.M., Kumar, M., Sonnenthal, E., Bao, C., Beisman, J., 2017. Expanding the role of reactive transport models in critical zone processes. Earth Sci. Rev. 165, 280-301. https://doi.org/10.1016/j.earscirev.2016.09.001.
Looms, M.C., Jensen, K.H., Binley, A., Nielsen, L., 2008. Monitoring unsaturated flow and transport using cross-borehole geophysical methods. Vadose Zone J. 7, 227–237. https://doi.org/10.2136/vzj2006.0129.

Lu, X., Wang, H., 2012. Microbial oxidation of sulfide tailings and the environmental consequences. Elements 8, 119–124. https://doi.org/10.2113/gselements.8.2.119.

Mayer, K.U., Alt-Epping, P., Jacques, D., Arora, B., Steefel, C.I., 2015. Benchmark problems for reactive transport modeling of the generation and attenuation of acid rock drainage. Comput. Geosci. 19, 599–611. https://doi.org/10.1007/s10596-015-9476-9.

McCarty, D.K., Moore, J.N., Marcus, W.A., 1998. Mineralogy and trace element association in an acid mine drainage iron oxide precipitate; comparison of selective extractions. Appl. Geochem. 13, 165–176. https://doi.org/10.1016/S0883-2927(97)00067-X.

Mertz, S., Le Forestier, L., Bataillard, P., Devau, N., 2021. Leaching of trace metals (Pb) from contaminated tailings amended with iron oxides and manure: new insight from a modelling approach. Chem. Geol. 579, 120356. https://doi.org/10.1016/j.chemgeo.2021.120356.

Moncur, M.C., Ptacek, C.J., Blowes, D.W., Jambor, J.L., 2005. Release, transport and attenuation of metals from an old tailings impoundment. Appl. Geochem. 20, 639–659. https://doi.org/10.1016/j.apgeochem.2004.09.019.

Morbidelli, R., Saltalippi, C., Flammini, A., Govindaraju, R.S., 2018. Role of slope on infiltration: a review. J. Hydrol. 557, 878–886. https://doi.org/10.1016/j.jhydrol.2018.01.019.

Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12, 513–522. https://doi.org/10.1029/WR012i003p00513.

Muniruzzaman, M., Karlsson, T., Kauppila, P.M., 2018. Modelling tools for the prediction of drainage quality from mine wastes. Geological Survey of Finland 408, 27–42. https://doi.org/10.30440/bt408.2.

Ouangrawa, M., Molson, J., Aubertin, M., Bussière, B., Zagury, G.J., 2009. Reactive transport modelling of mine tailings columns with capillarity-induced high water saturation for preventing sulfide oxidation. Appl. Geochem. 24, 1312–1323. https://doi.org/10.1016/j.apgeochem.2009.04.005.

Pabst, T., Molson, J., Aubertin, M., Bussière, B., 2017. Reactive transport modelling of the hydro-geochemical behaviour of partially oxidized acid-generating mine tailings with a monolayer cover. Appl. Geochem. 78, 219–233. https://doi.org/10.1016/j.apgeochem.2017.01.003.

Palandri, J.L., Kharaka, Y.K., 2004. A Compilation of Rate Parameters of Water-Mineral Interaction Kinetics for Application to Geochemical Modeling (No. 2004–1068).

Parkhurst, D.L., Appelo, C.A.J., 2013. Description of input and examples for PHREEQC version 3: a computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations (report no. 6-A43). In: techniques and methods. Reston, VA. doi: https://doi.org/10.3133/tm6A43.

Plumlee, G., Morman, S., 2011. Mine wastes and human health. Elements 7, 399–404. https://doi.org/10.2113/gselements.7.6.399.

Runkel, R.L., Kimball, B.A., 2002. Evaluating remedial alternatives for an acid mine drainage stream: application of a reactive transport model. Environ. Sci. Technol. 36, 1093–1101. https://doi.org/10.1021/es0109794.

Sand, W., Jozsa, P.-G., Kovacs, Z.-M., Sasāran, N., Schippers, A., 2007. Long-term evaluation of acid rock drainage mitigation measures in large lysimeters. J. Geochem. Explor. 92, 205–211. https://doi.org/10.1016/j.gexplo.2006.08.006.

Shaw, S.C., Groat, L.A., Jambor, J.L., Blowes, D.W., Hanton-Fong, C.J., Stuparyk, R.A., 1998. Mineralogical study of base metal tailings with various sulfide contents, oxidized in laboratory columns and field lysimeters. Environ. Geol. 33, 209–217. https://doi.org/10.1007/s002540050239.

Simate, G.S., Ndlovu, S., 2014. Acid mine drainage: challenges and opportunities. J. Environ. Chem. Eng. 2, 1785–1803. https://doi.org/10.1016/j.jece.2014.07.021.

Šimúnek, J., Genuchten, M.Th., Sejna, M., 2008. Development and applications of the HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7, 587–600. https://doi.org/10.2136/vzj2007.0077.

Smith, A.M.L., Dubbin, W.E., Wright, K., Hudson-Edwards, K.A., 2006. Dissolution of lead- and lead-arsenic-jarosites at pH 2 and 8 and 20 °C: insights from batch experiments. Chem. Geol. 229, 344–361. https://doi.org/10.1016/j.chemgeo.2005.11.006.

Steefel, C., Depaolo, D., Lichtner, P., 2005. Reactive transport modeling: an essential tool and a new research approach for the earth sciences. Earth Planet. Sci. Lett. 240, 539–558. https://doi.org/10.1016/j.epsl.2005.09.017.

Steefel, C.I., Beckingham, L.E., Landrot, G., 2015. Micro-continuum approaches for modeling pore-scale geochemical processes. Rev. Mineral. Geochem. 80, 217–246. https://doi.org/10.2138/rmg.2015.80.07.

Thouin, H., Battaglia-Brunet, F., Gautret, P., Le Forestier, L., Breeze, D., Séby, F., Norini, M.P., Dupraz, S., 2017. Effect of water table variations and input of natural organic matter on the cycles of C and N, and mobility of as, Zn and Cu from a soil impacted by the burning of chemical warfare agents: a mesocosm study. Sci. Total Environ. 595, 279–293. https://doi.org/10.1016/j.scitotenv.2017.03.218.

Thouin, H., Battaglia-Brunet, F., Norini, M.P., Le Forestier, L., Charron, M., Dupraz, S., Gautret, P., 2018. Influence of environmental changes on the biogeochemistry of arsenic in a soil polluted by the destruction of chemical weapons: a mesocosm study. Sci. Total Environ. 627, 216–226. https://doi.org/10.1016/j.scitotenv.2018.01.158.

Thouin, H., Norini, M.P., Battaglia-Brunet, F., Gautret, P., Crampon, M., Le Forestier, L., 2022. Temporal evolution of surface and sub-surface geochemistry and microbial communities of Pb-rich mine tailings during phytostabilization: a one-year pilot-scale study. J. Environ. Manag. 318, 115538. https://doi.org/10.1016/j.jenvman.2022.115538.

Thouin, H., Norini, M.P., Le Forestier, L., Gautret, P., Motelica-Heino, M., Breeze, D., Gassaud, C., Battaglia-Brunet, F., 2019. Microcosm-scale biogeochemical stabilization of Pb, as, Ba and Zn in mine tailings amended with manure and ochre. Appl. Geochem. 111, 104438. https://doi.org/10.1016/j.apgeochem.2019.104438.

van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898. https://doi.org/10.2136/sssaj1980.03615995004400050002x.

Walter, A.L., Frind, E.O., Blowes, D.W., Ptacek, C.J., Molson, J.W., 1994. Modeling of multicomponent reactive transport in groundwater: 2. Metal mobility in aquifers impacted by acidic mine tailings discharge. Water Resour. Res. 30, 3149–3158. https://doi.org/10.1029/94WR00954.

Wang, X., Li, Y., Wang, Y., Liu, C., 2018. Performance of HYDRUS-1D for simulating water movement in water-repellent soils. Can. J. Soil Sci. 98, 407–420. https://doi.org/10.1139/cjss-2017-0116.

Welch, S.A., Kirste, D., Christy, A.G., Beavis, F.R., Beavis, S.G., 2008. Jarosite dissolution II—reaction kinetics, stoichiometry and acid flux. Chem. Geol. 254, 73–86. https://doi.org/10.1016/j.chemgeo.2008.06.010.

Wunderly, M.D., Blowes, D.W., Frind, E.O., Ptacek, C.J., 1996. Sulfide mineral oxidation and subsequent reactive transport of oxidation products in mine tailings impoundments: a numerical model. Water Resour. Res. 32, 3173–3187. https://doi.org/10.1029/96WR02105.

Zheng, C., Lu, Y., Guo, X., Li, H., Sai, J., Liu, X., 2017. Application of HYDRUS-1D model for research on irrigation infiltration characteristics in arid oasis of Northwest China. Environ. Earth Sci. 76, 785. https://doi.org/10.1007/s12665-017-7151-2.