Galton-Watson Tree

Thierry E Huillet

THE HEIGHT OF THE LATEST COMMON ANCESTOR OF TWO RANDOMLY CHOSEN LEAVES FROM A (SUB-)CRITICAL

Keywords: (sub-)critical Galton-Watson process, Birth and death, Pair coalescence time for two leaves

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.

Background: Ancestral and coalescence processes have received a lot of attention in the context of Moran branching processes conditioned on having fixed population size over the non-overlapping generations (in the sense of [START_REF] Karlin | Direct product branching processes and related Markov chains[END_REF]- [START_REF] Karlin | Direct product branching processes and related Markov chains I: Calculations of rates of approach to homozygosity[END_REF]), starting from the early work of [START_REF] Kingman | The coalescent[END_REF]. See for example [START_REF] Möhle | A classification of coalescent processes for haploid exchangeable population models[END_REF], [START_REF] Huillet | On the extended Moran model and its relation to coalescents with multiple collisions[END_REF] and [START_REF] Huillet | Asymptotics of Symmetric Compound Poisson Population[END_REF] and the references therein. There was also some recent interest on similar problems in the context of the standard unconditional Galton-Watson branching processes, either (sub)-critical or supercritical. In both cases, the problem is to characterize coalescence times distributions, such as the time to most recent common ancestor of a random sample of arbitrary size drawn from the current population, including full size. See [START_REF] Zubkov | Limiting distributions of the distance to the closest common ancestor[END_REF], [START_REF] Lambert | Coalescence times for the branching process[END_REF], [START_REF] Chen | Smaller population size at the MRCA time for stationary branching processes[END_REF], [START_REF] Athreya | Coalescence in critical and subcritical Galton-Watson branching processes[END_REF], [START_REF] Athreya | Coalescence in the recent past in rapidly growing populations[END_REF], [START_REF] Athreya | Ancestor Problem for Branching Trees[END_REF], [START_REF] Le | Coalescence times for the Bienaymé-Galton-Watson process[END_REF] and [START_REF] Grosjean | On the genealogy and coalescence times of Bienaymé-Galton-Watson branching processes[END_REF] for example. Specifically, the question of whether coalescence time, both for pairs of tips and for the whole population, occur in the recent past, in the distant past or in-between has been addressed. In this short note, we consider the modest and related problem of computing the distribution of the height of the latest common ancestor of two randomly chosen leaves from a complete (sub-)critical Galton-Watson branching tree. Leaves play an important role in the tree of leaves of [START_REF] Rizzolo | Scaling limits of Markov branching trees and Galton-Watson trees conditioned on the number of vertices with out-degree in a given set[END_REF] and the Brownian continuum random tree of [START_REF] Aldous | A critical branching process model for biodiversity[END_REF].

Branching processes are simple and fundamental models of non-linear birth and death population dynamics, [START_REF] Harris | The theory of branching processes[END_REF]. In a discrete-time Galton-Watson (GW) process, some founder initially gives birth to a random number M of offspring at the next generation, each daughter proceeding similarly independently of its sisters, and so on till all daughters are found sterile, [START_REF] Harris | The theory of branching processes[END_REF]. Let then

(1) E z M = φ (z) , z ∈ [0, 1]
be the branching probability generating function (pgf) of this GW process with φ (0) = 0 (explosion possible). We let µ = E (M ) = φ (1) < ∞ be the mean number of offspring per capita. We also let π m = P (M = m), the probability mass function of M . Starting from a single founder, the birth and death process is iterated indefinitely, ending up in a population with a total random amount N of descendants (the asymptotic total progeny of the founder at equilibrium i.e. after time was iterated indefinitely). A productive individual in the GW tree will produce M > 0 offspring at each step, including itself, meaning that a productive individual passes to the next generation and generates M -1 new daughters in such a reproduction event. If M = 1, the reproduction event is reduced to itself in a self-regeneration no-event (Equivalently, a productive individual, as an internal node of the tree, dies out at the reproduction event while giving birth to M > 0 offspring, the net production of newborn individuals there being M -1). As a result, the only individuals that matter after the whole lifetime of the tree was exhausted are the leaves (corresponding to nodes at which the death event M = 0 takes place, so with outdegree 0). The leaves are thus the sterile individuals contributing to the population extinction. Clearly then, denoting N m the number of nodes of the complete tree with outdegree m, the number of leaves N 0 with outdegree 0 obey ( 2)

N 0 = 1 + m≥1 (m -1) N m ,
together with of course N = m≥0 N m . The latter displayed formula indicates that the terminal number of leaves corresponds to the total number of newborns over the whole population lifetime.

Take a (sub)-critical GW branching process with µ ≤ 1, so going extinct with probability 1 in finite time. Its complete tree (i.e. after its extinction) will present finitely many nodes N and leaves N 0 a.s.. If this tree has at least two leaves, picking two distinct leaves at random, we ask for the height H (as measured from the root) of their latest common ancestor. The height H, as a coalescence time between pairs, is a global measure of the proximity between any two leaves of the complete GW tree.

If the founder individual at the root generates M ≥ 2 offspring, the complete tree is a concatenation of M first-generation sub-trees which are iid copies of the complete tree. Upon conditioning on the first branching event M at the root, with H a copy of H, we have the characterization in distribution ( 3)

H d =        ∞ if M = 0 H + 1 if M = 1 0 if M ≥ 2 and the 2 leaves / ∈ same sub-tree H + 1 if M ≥ 2 and the 2 leaves ∈ same sub-tree.
Here, H was set to ∞ if the tree is made of only one leaf.

Given M = m ≥ 2, we let

P m = P (2 leaves ∈ same sub-tree | M = m) .
With ϕ (z) = E z H the pgf of H, we thus have from ( 3)

ϕ (z) =   π 1 zϕ (z) + m≥2 π m [(1 -P m ) + P m zϕ (z)]   so that, with P = m≥2 π m P m ϕ (z) = 1 -(π 0 + π 1 + P ) 1 -z (π 1 + P ) . Note ϕ (1) = 1 -(π 0 + π 1 + P ) 1 -(π 1 + P ) = P (H < ∞) < 1 and (4) E z H | H < ∞ = ϕ (z) ϕ (1) = 1 -(π 1 + P ) 1 -z (π 1 + P ) .
We are left with the problem of computing P m and P . Our main result is Theorem 1. With Φ 0 (z) the pgf of the total number of leaves N 0 in a complete GW tree generated by φ given by the solution of Eq. ( 5),

P = m≥2 π m P m = (1 -π 1 ) (1 -Φ 0 (0)) -π 0 1 0 log Φ 0 (z) Φ 0 (0) dz.
Given H < ∞, from Eq. ( 4), H is geometrically distributed with mean

m = E (H | H < ∞) = π 1 + P 1 -(π 1 + P ) . Proof: Let (n) 2 = n (n -1)
. Given M = m ≥ 2, the probability that the two distinct leaves belong to the sub-tree m ∈ {1, ..., m} is

E     N (m ) 0 2 m m =1 N (m ) 0 2     ,
where N (m ) 0

; m ∈ {1, ..., m} are the iid number of leaves of each of the m subtrees, each distributed like N 0 . Note that this probability is 0 for a sub-tree with N (m ) 0 = 1 (i.e. a sub-tree reduced to a single leaf). So P m = mE

(N 0 ) 2 m m =1 N (m ) 0 2
.

We now have (see [START_REF] Grosjean | On the genealogy and coalescence times of Bienaymé-Galton-Watson branching processes[END_REF], Lemma 9)

P m = m Γ (2) 1 0 (1 -z) Φ 0 (z) m-1 Φ 0 (z) dz, where (5)
Φ 0 (z) = π 0 (z -1) + φ (Φ 0 (z)) , Φ 0 (0) = 0 and Φ 0 (1) = 1
is the functional equation giving the pgf of the total number of leaves N 0 in the tree generated by φ, (see [START_REF] Drmota | Random Trees: An Interplay between Combinatorics and Probability[END_REF], p. 84; for example). Note

Φ 0 (z) = π 0 + φ (Φ 0 (z)) Φ 0 (z) , leading to Φ 0 (0) = P N 0 = 1 = π 0 1 -π 1 and Φ 0 (1) = EN 0 = π 0 1 -µ ≤ ∞.
(with EN 0 → ∞ as µ → 1 -, while approaching the critical regime). We thus have m≥2

π m P m = 1 0 (1 -z) m≥2 mπ m Φ 0 (z) m-1 Φ 0 (z) dz = 1 0 (1 -z) φ (Φ 0 (z)) -π 1 Φ 0 (z) dz = 1 0 (1 -z) (1 -π 1 ) Φ 0 (z) -π 0 Φ 0 (z) Φ 0 (z) dz = (1 -π 1 ) 1 0 (1 -z) dΦ 0 (z) -π 0 1 0 (1 -z) d log Φ 0 (z) .
By integration by parts

1 0 (1 -z) dΦ 0 (z) = [(1 -z) Φ 0 (z)] 1 z=0 + 1 0 Φ 0 (z) dz = 1 -Φ 0 (0) = 1 -(π 0 + π 1 ) 1 -π 1 1 0 (1 -z) d log Φ 0 (z) = [(1 -z) log Φ 0 (z)] 1 z=0 + 1 0 log Φ 0 (z) dz = 1 0 log Φ 0 (z) Φ 0 (0) dz Thus P = m≥2 π m P m = (1 -π 1 ) (1 -Φ 0 (0)) -π 0 1 0 log Φ 0 (z) Φ 0 (0) dz.
Example (binary tree): consider the branching mechanism φ (z

) = π 0 + π 1 z + π 2 z 2 with π 0 + π 1 + π 2 = 1 and µ = 2π 2 + π 1 = 1 + π 2 -π 0 ≤ 1 (else π 2 ≤ π 0 ). We have Φ 0 (z) = 1 -π 1 -(1 -π 1 ) 2 -4π 0 π 2 z 2π 2 Φ 0 (z) Φ 0 (0) = 1 -4 π 0 π 2 (1 -π 1 ) 2 z -1/2
.

Therefore, if we first restrict our attention to µ < 1 (subcritical case)

1 0 log Φ 0 (z) Φ 0 (0) dz = - (1 -π 1 ) 2 8π 0 π 2 1 1-4 π 0 π 2 (1-π 1 ) 2 log z dz = - (1 -π 1 ) 2 8π 0 π 2 [z log z -z ] 1 z =1-4 π 0 π 2 (1-π 1 ) 2 = (1 -π 1 ) 2 8π 0 π 2 1 -4 π 0 π 2 (1 -π 1 ) 2 log 1 -4 π 0 π 2 (1 -π 1 ) 2 + 1 2 = 1 2 + (π 0 -π 2 ) 2 4π 0 π 2 log π 0 -π 2 1 -π 1 .
As a result,

P = m≥2 π m P m = 1 -(π 0 + π 1 ) - π 0 2 - (π 0 -π 2 ) 2 4π 2 log π 0 -π 2 1 -π 1 . Q : = P 1 -(π 0 + π 1 ) = 1 - π0 2 + (π0-π2) 2 4π2 log π0-π2 1-π1 1 -(π 0 + π 1 ) = 1 - π 0 2π 2 + π 0 -π 2 2π 2 2 log π 0 -π 2 π 0 + π 2 = π 0 -π 2 2π 2 2 log π 0 + π 2 π 0 -π 2 - π 0 2π 2 -1 . Note 0 < Q < 1/2 and Q only depends on the ratio π 0 /π 2 . Indeed, if π 0 = aπ 2 with a > 1, Q = a-1 2 2 log a+1 a-1 -a-2
2 which is convex decreasing from .5 to 0 as a varies from 1 to ∞. The range a → ∞ corresponds to π 2 → 0.

In the critical case where µ → 1 -, = π 0 -π 2 → 0 + and (observing -2 log → 0 as → 0 + ) the above expression of Q can be prolonged by continuity, leading to Q → 1 2 .

In any case (see Figure 1 for a plot of P versus a and π 2 ),

0 < P < (1 -(π 0 + π 1 )) /2 = π 2 /2 and 0 < (1 -(π 0 + π 1 )) /2 < 1 -(π 0 + π 1 + P ) < 1 -(π 0 + π 1 ) < 1 0 < (1 -(π 1 -π 0 )) /2 < 1 -(π 1 + P ) < 1 -π 1 < 1, with m := π 1 + P 1 -(π 1 + P ) = π 1 + Qπ 2 1 -π 1 -Qπ 2 .
-If π 1 is large, close to 1, π 2 is small together with Q and m is close to π 1 / (1 -π 1 ), so large.

-On the other extreme hand, if

π 1 = 0, Q = - 2π 0 -1 2 (1 -π 0 ) 2 log (2π 0 -1) - π 0 2 (1 -π 0 ) -1 m = Q (1 -π 0 ) / (1 -Q (1 -π 0 ))
where π 0 ≥ 1/2. The mean height m = E (H) is a decreasing function of π 0 taking values between 1/3 (in the critical case when π 0 = 1/2) and 0, tending to 0 as π 0 → 1.

For a general branching mechanism φ, Φ 0 is not available in closed form as in the above binary example, neither is the integral appearing in the expression of P . However a series expansion of P is available. Indeed, Eq. ( 5) Φ 0 (z) = π 0 (z -1) + φ (Φ 0 (z)) , Φ 0 (0) = 0 and Φ 0 (1) = 1 is also

f (Φ 0 (z)) = z, where f (z) = 1 - 1 π 0 (φ (z) -z) with f (0) = 0. We further have Φ 0 (z) = π 0 / 1 -φ (Φ 0 (z)) , Φ 0 (0) = π 0 / (1 -π 1 ), so log Φ 0 (z) Φ 0 (0) = Ψ (Φ 0 (z))
where Ψ (z) = log (1 -π 1 ) / 1 -φ (z) . By Lagrange inversion formula (see [START_REF] Comtet | Analyse Combinatoire. Tomes 1 et 2[END_REF], Part I, p. 160; for example)

[z n ] Ψ (Φ 0 (z)) = 1 n z n-1 Ψ (z) f (z) z -n = 1 n z n-1 φ (z) 1 -φ (z) f (z) z -n
.

The computation of .

This gives from (6) the series expansion of P and so, from Eq. ( 4), the general law of H | H < ∞.

Concluding remarks:

The question of finding the height of the latest common ancestor of N 0 > i > 2 randomly chosen leaves from a (sub-)critical GW tree is left open so far as the conditioning on the first branching event M at the root that guarantees that the complete GW tree has at least i leaves is not so straightforward. Neither is also the problem of determining to which combination of sub-trees these selected leaves belong.

( 6 )

 6 P = (1 -(π 0 + π 1 )) -

Figure 1 .

 1 Figure 1. P versus a and π 2 .

Acknowledgments:

The author acknowledges partial support from the "Chaire Modélisation mathématique et biodiversité" of Veolia-Ecole Polytechnique-MNHN-Fondation X. He also acknowledges support from the labex MME-DII Center of Excellence (Modèles mathématiques et économiques de la dynamique, de l'incertitude et des interactions, ANR-11-LABX-0023-01 project).


Soumaia M'nassri

Yann Lucas

Gerhard Schäfer
email: schafer@unistra.fr

Lotfi Dridi

Rajouene Majdoub

Coupled hydrogeochemical modelling using KIRMAT to assess water-rock interaction in a saline aquifer in central-eastern Tunisia

Keywords: groundwater salinization, rock-water interaction, reactive mass transport, Sidi El Hani aquifer, numerical model

To quantify water-rock interactions and forecast their spatiotemporal evolution in the shallow aquifer of Sidi El Hani located in central-eastern Tunisia, the numerical code KIRMAT was used. Based on available geological, mineralogical, thermodynamic, kinetic and hydrodynamic data, the model was built and calibrated on two transects of the aquifer oriented along the flow lines of groundwater. The numerical results show that the interactions between water and rock are responsible for groundwater salinity, which is strongly dominated by chlorides, sodium and sulphates. The variations of dissolved species in space and time are the result of halite dissolution, anhydrite and gypsum dissolution and precipitation, calcite dissolution and dolomite precipitation. The precipitation of dolomite can locally induce a reduction in permeability and further enhance an increase in leaching because the residence time between rock and water will be amplified. Three numerical sensitivity studies have been conducted to evaluate the influence of water-rock interactions on groundwater salinization.

The study highlights the finding that the modelled fate of dissolved species in the unconfined aquifer is rather sensitive to the Darcy velocity but not very sensitive to the initial volume fractions of minerals in the two transects and the composition of groundwater injected at the upstream boundary.

Introduction

In arid and semi-arid climate regions, groundwater represents the primary source of water for agricultural, industrial and domestic uses [START_REF] Li | Assessement of groundwater quality for irrigation purpose and identification of hydro-geochemical evolution mechanisms in Pengyang County, China[END_REF]. Unfortunately, since the last few 2 decades, these resources have been affected by salinization problems, which occur in many aquifers around the world (De Montaty et al., 2008;[START_REF] Kumar | Assessment of hydrochemistry and groundwater quality in the coastal area of South Chennai, India[END_REF], as well as in the Mediterranean areas [START_REF] Farid | Origin of the groundwater salinity and geochemical processes in detrial and carbonate aquifers: case of Chougafiya basin (Central Tunisia)[END_REF][START_REF] Eissa | Groundwater recharge and salinization in the arid coastal plain aquifer of the Wadi Watir delta, Sinai, Egypt[END_REF]. In these regions, the reasons for the degradation of the groundwater quality are mainly seawater intrusion, climate change, irrigation return flow, uncontrolled fertilizer applications, dissolution of evaporites and upwelling of saline waters from deep layers [START_REF] Milnes | Process-based groundwater salinisation risk assessment methodology: Application to the Akrotiri aquifer (Southern Cyprus)[END_REF][START_REF] Salem | Hydrochemical and isotope evidence of groundwater contamination of cultivated fields of semi-arid environments in Tunisia[END_REF][START_REF] Li | Accelerate research on land creation[END_REF]Giambastiani et al., 2014;[START_REF] Han | Chemical and isotopic constraints on evolution of groundwater salinization in the coastal plain aquifer of Laizhou bay, China[END_REF][START_REF] Najib | Salinization process and coastal groundwater in Chaouia, Morroco[END_REF]. The interaction of these processes can lead to a varied composition of dissolved species in the aquifer.

Therefore, eliminating processes that alter groundwater quality is crucial for powerful water management and protection and necessary to ensure the sustainability of precious groundwater resources in arid regions [START_REF] Sikdar | Geochemical evolution of groundwater in the Quaternary aquifer of Calcutta and Howrah, India[END_REF]Moran-Ramírez et al., 2016).

In recent years, a strong effort has been undertaken to thoroughly develop hydrogeochemical tools to assess groundwater quality, such as statistical and geostatistical approaches [START_REF] Alberto | Pattern recognition techniques for the evaluation of spatial and temporal variations in water quality. A case study: Suqua River Basin (Cordoba-Argentina)[END_REF][START_REF] Belkhiri | Application of multivariate statistical methods and inverse geochemical modeling for characterisation of groundwater-A case study : Ain Azel plain (Algeria)[END_REF], the empirical measurement of major chemical indicators, the design of hydrogeochemical diagrams [START_REF] Edmunds | The natural (baseline) quality of groundwater: aUK pilot study[END_REF]Andre et al., 2005) and water analysis based on isotopic techniques [START_REF] Bouchaou | Application of multiple isotopic and geochemical tracers for investigation of recharge, salinization and residence time of water in the Souss-Massa aquifer, southwest of Morocco[END_REF]. These tools provide a better understanding of the interference of physical and chemical processes that affect water chemistry and the residence time of the major dissolved species involved in the system. However, each approach has its own limitations due to the multiplicity of geochemical processes between the different phases, which adds complexity to the analysis. Additionally, the mass fluxes in the aquifer and the chemical composition of the rock may strongly vary in time and space [START_REF] Xiao | Solute gochemsitry and its sourcs of the groundwaters in the Qinghai Lake catchment, NW china[END_REF][START_REF] Peikam | Application of inverse geochemical modeling for predicting surface water chemistry in Ekbatan watershed, Hamedan, Western Iran[END_REF]. One way to overcome these problems is the use of numerical reactive transport models [START_REF] Beaucaire | Groundwater characterisation and modelling of water-rock interaction in an argillaceous formation (Tournemire, France)[END_REF][START_REF] Lucas | Geochemical tracing and hydrogeochemical modeling of water-rock interaction during salinization of alluvial groundwater (Upper Rhine Valley, France)[END_REF][START_REF] Montcoudiol | A Conceptual model for groundwater flow and geochemical evolution in the southern Outaouais region, Québec Region, Québec, Canada[END_REF][START_REF] Jing | Hydro-geochemical for the evolution of groundwater quality in Yinchuan plain, China[END_REF][START_REF] Nasri | Modeling the hydrogeochemical evolution of brine in saline systems: Case study of the sabkha of Oum El Khialate in South East Tunisia[END_REF].

The modelling approach is an effective tool to study and predict the degradation of groundwater quality [START_REF] Lopez-Chicano | Factors which determinate the hydrogeochemistry behavior of karstic springs: a case study from the Betic Cordilleras, Spain[END_REF][START_REF] André | Using geochemical data and modelling to nhance the understanding of groundwater flow in a regional deep aquifer, Aquitaine Basin, south-west of France[END_REF]Nilde et al., 2008;[START_REF] Li | Geochemical modeling of groundwater in southern plain area of Pengyang County, Ningxia, China[END_REF]. In the past, a large number of geochemical programmes have been developed. Among these, some of them attempt to determine the final thermodynamic equilibrium of the geochemical system [START_REF] Helgeson | Evaluation of irreversible reactions in geochemical processes involving minerals ans aqueous solutions. I. Thermodynamic relations[END_REF], while others take into account the phase mole transfer throughout kinetic reactions [START_REF] Beig | Albite dissolution kinetics as a function of distance from equilibrium: implications for natural feldspar weathering[END_REF][START_REF] Burch | Free energy dependence of albite dissolution kinetics at 80°C and pH 8.8[END_REF]. Such models are useful because they provide a thorough understanding of the aquifer system. However, their use remains somewhat limited for modelling mineral digenesis coupled with dynamic processes of mass transport (Nourtier-Mazauric et al., 2005;[START_REF] Nourtier-Mazauric | Modélisation géochimique et numérique des interactions entre des solutions solides et une solution aqueuse[END_REF]. Therefore, more advanced numerical tools need to be implemented for reactive transport models [START_REF] Gérard | Modélisation géochimique thermodynamique et cinétique avec prise en compte des phénomènes de transport de masse en milieux poreux saturé[END_REF][START_REF] Guo | The geochemistry simulation of groundwater chemical environment evolution in Hebei plain[END_REF][START_REF] Eary | Geochemical controls on groundwater composition at the Cripple Creek Mining District, Colorado[END_REF][START_REF] Li | Hydro-geochemical modeling on saturation of minerals in groundwater. A case study at Niangzigun, Northern China[END_REF][START_REF] Lakshmanan | Major ion chemistry and identification of hydro-geochemical processes of groundwater in a part of Kancheepuran district, Tamil Nadu, India[END_REF]Appelo et al., 2008;[START_REF] Belkhiri | Application of multivariate statistical methods and inverse geochemical modeling for characterisation of groundwater-A case study : Ain Azel plain (Algeria)[END_REF][START_REF] Steefel | Reactive transport codes for subsurface environmental simulation[END_REF][START_REF] Laabidi | Reactive Henry problem: effect of calcite dissolution on seawater intrusion[END_REF]. They may be used to quantify the geochemical processes and predict possible changes in the composition of both water and rock and the fate of dissolved species in the aquifer.

Therefore, advanced quantitative studies using numerical hydrogeochemical models are needed to understand and study the horizontal extent of the saturated zone that may be influenced by geologic alteration.

This paper aims to quantify the main processes of salinization of the groundwater in the basin of Sidi El Hani. A coupled hydrogeochemical model was used to study and predict the chemical changes that occur as water evolves along a flow path. Therefore, we used the multi-component reactive transport code KIRMAT (Kinetic Reaction and Mass Transport) [START_REF] Gérard | Introduction des phénomènes de transport dans le modèle thermo-cinétique KINDIS : le modèle KIRMAT[END_REF][START_REF] Gérard | Numerical validation of an Eurlian hydrochemical code using a 1-D multisolute mass transport system involving heterogeous kineticallycontroled reactions[END_REF] that combines geochemical reactions and 1D mass transport equations.

The main objectives are to characterize and identify the possible geochemical processes acting in the groundwater and to predict the evolution of the chemical composition. To meet these objectives, a multi-component reactive transport model KIRMAT was applied in this study. The Sidi El Hani aquifer, where the salinity exceeds 6 g.L -1 to reach 8 g.L -1 , is therefore a very suitable site to evaluate the influence of geochemical process on groundwater salinization. In the first step, the response of the numerical model is evaluated by the rootmean-square-error (RMSE) between the computed and observed concentrations of each monitoring point. Second, once the numerical model is calibrated, a sensitivity analysis is performed to identify the model input parameters that may have an influence on the reactive mass transport in the aquifer. This allows us to understand and interpret the numerical results of the hydrochemical model. The sensitivity study is conducted for three specific cases: i) variation of the volume fraction of halite; ii) variation of the chemical composition of incoming water and (iii) variation of the Darcy velocity.

Material and methods

Study area, geological and hydrogeological setting

The investigated area is located in the centre-eastern part of Tunisia, between 35°21'-35°35'N and 10°10'-10°25'E, covering approximately 346 km 2 . This area extends between the mountain chains of El Guessat and Ktitir, which barely reach more than 100 m of altitude in the western part, and the lowlands of a continental sabkha named the Sidi El Hani in the eastern and Cherita in the southern regions. The largest and permanent hydrological network in this region is Wadi Chertia, which runs from Sabkha Cherita to Sabkha Sidi El Hani (Fig. 1). The plain lies within a continental arid to semi-arid climate region. The annual average rainfall is below 270 mm. The monthly average temperature varies between 9°C and 30°C in January and August, respectively. The evapotranspiration is estimated to be 1500 mm per year. As in many arid and semi-arid regions of the world, the surface waters are scarce. The aquifer is considered to be one of the main water resources of the agricultural region. The plain has been intensively cultivated for the last decade, and 75% of the area is occupied by trees, crops and cereals [START_REF] Majdoub | Etude de la dynamique de l'eau et du transfert des sels dans un sol sablo-limoneux du sahel Tunisien[END_REF]. However, the groundwater is no longer suitable for agricultural uses because the total concentration of dissolved salts (chlorides, sodium, sulphates) generally exceeds the limits values for irrigation water (M'nassri et al., 2016). The Ouled Chamekh plain is Mio-Plio-Quaternary in age, comprising the Holocene (Qh), the Pleistocene (Qp) and the Mio-Pleistocene (Qmp) (Fig. 2). The Holocene formation is largely present in the study area as clays, sands, silts and evaporites [START_REF] Tagorti | Geochemical controls of groundwater upwelling in saline envirnments : Case study the discharge playa of Sidi El Hani (Sahel, Tunisia)[END_REF]. It represents the superficial layer. The Pleistocene deposits can be divided into three layers. The first layer is generally below 35 m and is composed of sand, halite and carbonate deposits.

The second layer is limestone and is between 35 and 45 m thick. The third layer is hard limestone, soft white limestone and silts and is 45 to 50 m thick. The Mio-Pleistocene is composed of sand with an interlayered 50 m-thick limestone. In addition, the study area lies within a collapsed zone with a synclinal form and cuttings of saline depressions along NW-SE orientations [START_REF] Essefi | Geodynamic framework of saline systems in Eastern Tunisia: Saline Depressions inherited from the Triassic intrusion and/or the Messinian Salinity Crisis[END_REF]. Furthermore, the Ouled Chamekh plain is sliced by a fault network that controls the Mio-Plio-Quaternary deposition, including the Ktitir, Hajeb Layoun and Sidi El Hani faults. These faults are inherited from the Cenozoic and Mesozoic formations (Ghribi, 2010;[START_REF] Essefi | Hydrocarbons migration through groundwater convergence towards saline depression. A case study Sidi El Hani, discharge playa, Tunisian Sahel[END_REF]. The Sidi El Hani fault causes subsidence of the Sabkha Sidi El Hani along the eastern boundary of the study area. According to the mineralogical composition, the study area has quartz (60-70%), carbonate minerals and evaporites (15-20%) and clay (5-10%). The clays are mainly composed of kaolinite, illite and muscovite. The depth of the water table in this shallow aquifer ranges from 3 m at the eastern part of the basin up to 26 m in the centre. The groundwater flow direction is oriented from the west to the Sabka Sidi El Hani in the north-east and to the Sabkha Cherita in the south (Fig. 3). The hydraulic heads vary between 80 and 35 m above sea level. Irregularities in the water table contours indicate high variability in the hydraulic conductivity of the alluvial aquifer from the west to the south and from the west to the east due to the heterogeneity of the sediments. In addition, the transmissivity varies between 2 10 -4 and 5 10 -3 m 2 .s -1 , and the porosity varies between 3.5 and 4% [START_REF] Dridi | Modélisation hydrogéologique des écoulements souterrains au niveau de la région d'OuledChamekh[END_REF].

Numerical reactive transport model

Governing equations of reactive mass transport

Most coupled reactive transport models rely on the resolution, simultaneous or not, of the partial differential equation in a saturated porous medium and the equation resulting from geochemical phenomena, which is added to the transport equation as a sink / source term.

The classical hypotheses of instantaneous homogeneous reactions with equal diffusion coefficients of all present aqueous species that only consider the kinetics of heterogeneous reactions make it possible to express the mass balance equation in terms of generalized (or total) concentrations Ψ (in mole per water mass or volume) of primary species. Under these conditions, the mass balance equation of reactive transport in a one-dimensional porous medium is written as (Lichtner, 1988;[START_REF] Gérard | Numerical validation of an Eurlian hydrochemical code using a 1-D multisolute mass transport system involving heterogeous kineticallycontroled reactions[END_REF]:

ϕΨ = ϕD -U -∑ α ϕ V j = 1, … . , N$ (1)
and

ϕ V = v r = 1, … . , M$, (2)
where Equation (1) refers to N aqueous primary species and Equation (2) refers to M primary species of reacting minerals. In these equations, ϕ denotes the porosity of the porous medium, ϕ and α denote the volume fraction and the stoichiometric reaction coefficients, respectively, of the r th mineral with molar volume V ( , v represents the reaction rates of the irreversible reaction of minerals and fluids equivalent to the rate of precipitation or dissolution of reacting minerals r per unit of the rock and fluid system (by convention, v is positive for precipitation and negative for dissolution reactions), D denotes the hydrodynamic dispersion coefficient, and U denotes the Darcy velocity. The generalized concentration Ψ is defined according to the expression:

Ψ = C + ∑ α + C + + , with C + = K + γ + ∏ γ C 3 4 5 , ( 3
)
where 6 refers to the concentration of the j th primary species, and the sum runs over all aqueous secondary species with concentration C + related to the concentrations of the primary species through the mass action equation. The quantity α + denotes the molar stoichiometric coefficient of species j in secondary species i, γ is the activity coefficient of the aqueous species, and K + denotes the equilibrium constant. Using (1) and ( 2), one obtains a system of (N +M) coupled nonlinear partial differential equations.

Description of the numerical model KIRMAT

The hydrochemical model developed in this study is based on the numerical code KIRMAT (KInetic Reaction and Mass Transport) (Gerard, 1996), which combines geochemical reactions and 1D mass transport equations. KIRMAT was developed from the thermokinetic section of the geochemical code KINDIS by integration of the reaction terms in an advectivedispersive scheme of mass transport in a saturated porous medium [START_REF] Madé | Modelling mineral/solution interactions: The thermodynamic and kinetic code KINDIS[END_REF]. The purpose of KIRMAT is to contribute to the understanding of groundwater salinization [START_REF] Lucas | Geochemical tracing and hydrogeochemical modeling of water-rock interaction during salinization of alluvial groundwater (Upper Rhine Valley, France)[END_REF], weathering processes [START_REF] Lucas | Hydrogeochemical modeling (KIRMAT) of spring and deep borehole water compositions in the small granitic Ringelbach catchment (Vosges Mountains, France)[END_REF][START_REF] Ackerer | Monitoring and reactive-transport modeling of the spatial and temporal variations of the Strengbach spring hydrochemistry[END_REF] and hydrothermal alteration [START_REF] Ngo | Modeling the impact of temperature on the saturation state and behavior of minerals in the Soultz-sous-Forêts geothermal system[END_REF]. It can also be used to simulate the fate of clay minerals under specific environmental conditions, such as the storage of nuclear waste (Montes et al., 2005a;2005b;[START_REF] Marty | Influence of reaction kinetics and mesh refinement on the numerical modelling of concrete/clay interactions[END_REF][START_REF] Marty | Modeling the long term alteration of the engineered bentonite barrier in an underground radioactive waste repository[END_REF]2006;[START_REF] Ngo | Coupled transport-reaction modelinf of the long-term interaction between iron, bentonite and Callovo-Oxfordian claystone in radioactive confinement systems[END_REF].

The geochemical modelling of the rock-water interaction via KIRMAT requires a thorough knowledge of input parameters, such as the groundwater chemical composition, the mineralogical properties of the aquifer, and the thermodynamic and kinetic data of the reaction involved. The workflow in KIRMAT is organized in three steps. In the first step, the ionic speciation (initial equilibrium) and the saturation index are calculated. During the second step, there is a perturbation of this initial partial equilibrium state by irreversible water-rock interaction due to dissolution or precipitation of minerals (reaction term vr, see Eq. 2), which results in either supplying or removing elements from the aqueous solution.

Finally, the (N+M) partial differential equations ( 1) and ( 2) are solved simultaneously with a one-step algorithm. The mass transport equations are solved using a finite difference approximation in an explicit scheme.

Saturation index and ionic speciation

The saturation index (SI) of the reactions occurring between the aqueous species and minerals is defined by Equation 4:

SI + = log < 4 = 4 , ( 4
)
where Qi denotes the ion activity product and Ki denotes the thermodynamic equilibrium For each species and mineral in the studied system, the ion activity products are obtained by applying the corresponding activity coefficients:

Q ? = ∏ [A ? ] 3+ = ∏ [m ? γ ? ] 3+ +∈E +∈E , ( 5
)
where Ai denotes the activity of species i, mi denotes the molality, γi denotes the activity coefficient, E represents the basic species assemblage of the study system and αij is the stoichiometric coefficient of the specie i in the complex j. In our current study, the activity coefficient γi has been fitted with the Debye-Hückel model ( 6), which is largely used in geochemical models [START_REF] Debye | Zur Theorie der Elektrolyte[END_REF]:

log γi = F G 4 √I JK L M √I + 6N, (6)
where A, B and C denote the Debye-Hückel constants, si, I and Zi denote the size, the ionic strength of the solution and the electric charge, respectively, of the specie i.

Water-rock interaction

In KIRMAT, the interaction between water and rock (see Eq. 2) involves dissolution and precipitation, which are controlled by many factors, including the mineral surface, the solution composition, the temperature and the degree of saturation of the aqueous solution with respect to the studied mineral. In our study, only the dissolution of minerals was modelled by irreversible kinetics. The kinetic of rock dissolution rd (in mole per water mass) is quantified by [START_REF] Gérard | Modélisation géochimique thermodynamique et cinétique avec prise en compte des phénomènes de transport de masse en milieux poreux saturé[END_REF]:

r O ≝ v = k O, S RSS α T U 1 - < V = V $, (7)
where kd,r denotes the dissolution rate constant (mole m -2 year -1 ) of the reactive mineral r, Sr eff stands for the reactive surface area of mineral r (m 2 .kgH2O -1 ), W X Y denotes the proton activity

where n depends on the pH of the solution, Qr is the ion activity product of reactive mineral r and Kr denotes the thermodynamic equilibrium constant of the hydrolysis reaction of mineral r at a given temperature and pressure.

Finally, the dynamic system is discretized over a time step Δt until the time t+Δt, which is equal to the simulation time requested by the user. A self-adjusting numerical time step was used during the simulations. Furthermore, the KIRMAT code can also describe the feedback effect of the chemical and mineralogical evolution of the permeability at any node of the mesh due to the dissolution and precipitation reaction, which can modify the pore geometry [START_REF] Laabidi | Nonstationary porosity evolution in mixing zone in coastal carbonate aquifer using an alternative modeling approach[END_REF]. The intrinsic permeability k (m 2 ) is updated after each time step as follows (Gerard, 1996):

Z = 6 [ [ ϕ \ ] ^ ]$ _ ] `, (8)
where C0, c, ϕ and S denote the cementing factor, an experimental constant, and the porosity of the porous medium and the grid cell surface in contact with the adjacent cell (m 2 ), respectively.

Model domain

According to the flow direction of the groundwater, two transects T1 and T2 are chosen to model the behaviour of the aquifer. The first transect is located from the lower (inlet) boundary near El Guessat Mountain to the Sabkha Sidi El Hani. The second transect is oriented to the Sabkha Cherita (Fig. 3).

The lengths of transect T1 and T2 are 16 km and 14 km, respectively. The two transects are discretized with a horizontal spacing of 300 m-long meshes. Based on the area and the geological characteristics (see Fig. 2), T1 contains four zones, whereas T2 is composed of three zones. The mineralogical compositions involved in the modelled transects are halite, calcite, gypsum, anhydrite, dolomite, illite and quartz. All of the zones of each transect are assumed to have the same mineralogy. The only exceptions are the first and second zones of T1 and T2, where halite is not present because the ion ratios of Na + to Cl -of the groundwater samples are lower than 1. This relationship suggests that the concentrations of both ions are not due to halite dissolution.

The simulation covered 50 and 100 year periods. These periods were chosen to explain the chemical composition of groundwater in 2015, which was used as the initial state. The two modelled transects are presented in Figure 3. In our study, the mineralogical compositions are considered homogenous in each transect zone due to the non-availability of the data. In addition, the geologic map showed significant lithological variation. All numerical simulations are performed at a constant temperature of 25°C because the average of the temperature measured in the different samples of groundwater is of the order of 25°C.

Initial and boundary conditions of the model

The initial conditions (IC) of our model are based on the average of the chemical composition of all observation wells located along each streamline. Concentrations of the upstream boundary (UB) were applied as fixed-concentrations corresponding to a 1 st order type, or Dirichlet boundary. The measured concentrations of the first observation well located at the upstream of the modelled section were used as the prescribed concentrations. The downstream boundary (DB) was set as a free flux boundary.

The physical and chemical characteristics of the injected fluid in T1 and T2 are reported in Table 1. The measurement of physical parameters such as pH and temperature (T) were performed in the field, whereas the analyses of chemical compositions of the injection fluid in both T1 and T2 were performed in the laboratory in accordance with the standard methods of the American Association of Public Health [START_REF] Apha | Standard method for the examination of water and wastewater, 19th ed[END_REF]. Calcium (Ca 2+ ) and magnesium (Mg 2+ ) concentrations were determined by the complexometric method. Sodium (Na + ) and potassium (K + ) concentrations were measured by atomic absorption spectrometry. Sulphate (SO4 2-) concentrations were determined by the gravimetric method, and chloride (Cl -) and bicarbonate (HCO3 -) concentrations were quantified by titration. Validation of the chemical analysis results was performed by verification of the error percentage, which was judged acceptable if the value is less or equal 5% [START_REF] Appelo | Geochemistry, groundwater and pollution[END_REF].

Model parameters

Mineralogical composition and reactive surface

To study the geochemical processes between the rock and groundwater in the Sidi El Hani aquifer, it is crucial to define the mineralogical composition of the aquifer. Hence, a borehole profile was performed in the study area with a depth of 30 m. There were 27 core samples collected along the profile. Rock samples were dried before undergoing mineralogical analysis. X-ray diffraction analysis of the cuttings highlights that the Sidi El Hani reservoir is heterogeneous. The reservoir contains quartz, gypsum, halite, anhydrite, dolomite, and calcite with 56%, 10%, 4%, 3%, 3% and 14% volume fractions, respectively (Table 2). In addition, the X-ray diffraction analysis underlined that the reservoir contains illite in a volume fraction of approximately 10%. Because illite has a low cation exchange capacity compared to the other clay minerals [START_REF] Sparks | Ion exchange processes[END_REF], sorption effects are neglected in our study. Heterogeneous reactions modelled with KIRMAT thus only include dissolution and precipitation processes between rock and water. The reactive surface of the minerals that is involved in both dissolution and precipitation reactions is often difficult to estimate because of insufficient knowledge of some natural factors. Note that this is one of the most sensitive parameters of the modelling of geochemical systems [START_REF] Wilson | Thermodynamic and fully-coupled reactive transport models of a steel-bentonite interface[END_REF]. Its estimation is often questionable due to certain ambiguities, such as the roughness of the surface and the changing of the surface area throughout the reaction progress (Savage et al., 2002;[START_REF] Molins | An investigation of the effect of pore scale flow on average geochemical reaction rates using direct numerical simulation[END_REF]. Methods for measuring reactive surfaces have been discussed in numerous studies [START_REF] Gérard | Modélisation géochimique thermodynamique et cinétique avec prise en compte des phénomènes de transport de masse en milieux poreux saturé[END_REF][START_REF] Kieffer | An experimental study of the reactive surface area of the Fontainebleau sandstone as a function of porosity, permeability and fluid flow rate[END_REF][START_REF] Kieffer | An experimental study of the reactive surface area of the Fontainebleau sandstone as a function of porosity, permeability and fluid flow rate[END_REF][START_REF] Marty | Modélisation couplée (transport-réaction) des interactions fluids-argiles et de leurs effets en retour sur les proprieties physiques de barriers ouvragées en bentonite[END_REF][START_REF] Ngo | Coupled transport-reaction modelinf of the long-term interaction between iron, bentonite and Callovo-Oxfordian claystone in radioactive confinement systems[END_REF]. In our study, we assume that the reactive surface area of the involved minerals remains constant during the numerical simulation. The reactive surface Sr eff (m 2 kg -1 H2O) was calculated as follows:

S RSS = a ( _ bcde f bcde , ( 9
)
where Vr (m 3 .kg -1 ) denotes the volume of the reactive mineral per kilogram of water, and Scube (m 2 ) and Vcube (m 3 ) stand for the surface and the volume of the cube, respectively. Vr was calculated by

V = F ]$ * [[[/j ] , (10)
where Fr and k represent the volume fraction and the volumetric mass (kg.m -3 ) of the mineral, respectively.

It is assumed that all minerals have a cubic shape. Therefore, the long sides should be 10 -3 m.

However, long sides of clay minerals are assumed to be 1 µm because they are formed with a smaller size than the other minerals [START_REF] Meunier | Why are the clay minerals small?[END_REF]. The porosity of the aquifer was estimated to be 0.04. The calculated reactive surfaces of the minerals are shown in Table 3.

Thermodynamic and kinetic parameters

The thermodynamic parameters of the minerals used in our numerical studies, at a temperature of 25°C, are summarized in Table 4 for six primary minerals and one secondary mineral. [START_REF] Gérard | Numerical validation of an Eurlian hydrochemical code using a 1-D multisolute mass transport system involving heterogeous kineticallycontroled reactions[END_REF] Note that kinetic data are less available in the literature than thermodynamic data, mainly for the precipitation of reactive minerals. Considering Equation 7, the dissolution rates of primary or secondary minerals were determined from data sets proposed by [START_REF] Palandari | A compilation of rate parameters of water-mineral kinetics for application to geochemical modeling[END_REF], which are based on kinetic values at 25°C and an activation energy term (Table 5). To model dissolution kinetics at any temperature, KIRMAT code uses the following expression for mineral dissolution controlled by surface reactions [START_REF] Madé | Modelling mineral/solution interactions: The thermodynamic and kinetic code KINDIS[END_REF]:

l m,? = n o pqq 1 - r L s L t u u u v Z `w xy?m exp } ~• €bL• ' ƒ -`"…. w † W xy?m Yxy?m +Z `w ‡pˆ‰(Š‹ OE•Ž } ~• •ec•'•' ' ƒ -`……. w † +Z `w "ŠKp OE•Ž } ~• "••e ' ƒ - `"…. w † W "ŠKp Y"ŠKp - - - - ˜, (11)
where rd,i represents the dissolution rate controlled by the surface reaction (mol year -1 kg - 1 H2O), and Seff is the reactive surface area for a given mineral i. Qi and Ki denote the activity product and the thermodynamic equilibrium constant at the given temperature, respectively.

k25 Acid , k25 Neutral and k25 Base represent the dissolution rate constant at 25°C and pH=0, for acidic, neutral and basic conditions, respectively. nAcid and nBase are the reactions order with respect to the pH of the solution, R the gas constant (J.K -1 mol -1 ) and T the temperature (K).

Hydrodynamic parameters

The hydrodynamic parameters of the Sidi El Hani aquifer were derived from the hydrogeological model created by [START_REF] Dridi | Modélisation hydrogéologique des écoulements souterrains au niveau de la région d'OuledChamekh[END_REF] using Visual Modflow. A Darcy velocity of approximately 91 and 152 m.year -1 for T1 and T2, respectively, were chosen according to hydraulic conductivity and hydraulic gradient measurements obtained from the contour map of hydraulic heads of the Sidi El Hani aquifer. In addition, the numerical value of the porosity for both transects (T1 and T2) was estimated to be 4% [START_REF] Dridi | Modélisation hydrogéologique des écoulements souterrains au niveau de la région d'OuledChamekh[END_REF].

However, the estimated values of hydraulic conductivity were approximately 30×10 -5 and 9×10 -5 m s -1 for T1 and T2, respectively.

Model calibration and sensitivity study

The calibration of the model was carried out by minimizing the differences between observed and simulated groundwater concentrations. The comparison was made with concentration data measured in March and April 2015 in numerous observation wells located along both transects. Using the reactive surface areas given in Table 3 and the initial dissolution rate constants (Table 5), the preliminary simulation led to an overestimation or underestimation of the chemical element concentrations. Hence, it was decided to adjust the dissolution rates of minerals, which are often different than those estimated by many authors [START_REF] Kim | Plagioclase weathering in the groundwater system of a sandy, silicate aquifer[END_REF]White and Brantely, 2003;[START_REF] Maher | The mineral dissolution rate conundrum: Insights from reactive transport modeling of Uisotopes using U isotopes and pore fluid chemistry in marine sediments[END_REF]. In the case of carbonates, we used the dissolution rate constants for basic conditions. The adjusted dissolution rates obtained for gypsum, anhydrite, calcite, dolomite and halite are shown in Table 6. They are in the range proposed by [START_REF] Palandari | A compilation of rate parameters of water-mineral kinetics for application to geochemical modeling[END_REF]. To obtain a good match between the computed and observed chemical composition, a large number of trial-and-error runs were performed. The response of the model was evaluated by the root-mean-square-error (RMSE) between computed and observed concentrations of each monitoring point (Eq. 12):

RMSE = › Y ∑ oe ? -6 ? $ `, (12)
where n denotes for the number of observations, and Oi and Ci are the measured and calculated groundwater concentrations of species i (mol kgH2O -1 ).

Sensitivity analysis was performed individually to identify the model input parameters that may have an influence on the reactive mass transport in the aquifer. This allows us to understand and interpret the numerical results of the hydrochemical model. The sensitivity study was conducted for three specific cases: i) variation of the volume fraction of halite ; ii) variation of the chemical composition of incoming water and (iii) variation of the Darcy velocity.

Results and discussion

Observed concentrations of major elements

The spatial variation of dissolved species measured in the groundwater samples along the two transects T1 and T2 is shown in Figure 4. The concentration of the major elements are Cl (0.02-0.07 mol.mol.kg H2O -1 ), Na (0.03-0.07 mol.kg H2O -1 ), SO4 (0.019-0.020 mol.kg H2O -1 ), Ca (0.007-0.023 mol.kg H2O -1 ), Mg (0.010-0.015 mol mol.kg H2O -1 ), HCO3 (0.003-0.006 mol.kg H2O -1 ), and K (0.001 mol mol.kg H2O -1 ) with the ionic sequence Na>Ca>Mg>K and Cl>SO4>HCO3. All of the water samples show a large deficit of HCO3 and K. However, a great excess of Cl and Na has been shown in zones 3 and 4 of T1 and in zone 3 of T2 compared to the first and second zones. In these areas groundwater is characterized by a high salinity with a concentration of total dissolved solids (TDS) above 6 g L -1 . In fact, the history of the salinity of the given aquifer system stems from the geomorphology of the area during the Quaternary. Essefi el al. (2013) showed that the Sidi El Hani-Cherita system had a salty heritage in the Tunisian Sahel area. Furthermore, they suggested two hypotheses for the origin of the salt. The first is linked to the halokinesis activities of the Triassic evaporitic sedimentation in the subsurface. The second states that the Messinian salinity crisis has left huge quantities of salt and salt water in the Sidi El Hani aquifer. This saline subsurface is until now influencing the behaviour of the area as a whole. Through groundwater convergence, the salts are accumulated within depression zones. Recently, M'nassri et al.

(2016 and 2018) investigated the salinization problem of the unconfined aquifer of Sidi El

Hani. Their studies showed that the high groundwater salinization is principally due to the rock-water interaction, including chemical and physical reactions, such as dissolution and precipitation. These reactions probably altered not only the chemical compositions of the rocks and the quality of the groundwater, but they also contributed to the formation of new minerals along the flow pathway (Edmunds et al., 2002;[START_REF] Rajmohan | Identification and volution of hydrogeochemical processes in the groundwater envirnment in area the Palar and Cheyyar River Basins, Southern India[END_REF]Meredith et al., 2012;[START_REF] Pacheco | Modeling rock weathering in small watersheds[END_REF]. Additionally, measured SO4 and Ca concentrations are increased, except in the water samples of zone 1. These two ions appear to be inherited from the same rock. The variation of the chemical elements in the different zones seems to be determined by the nature of the geology.

Simulated concentrations of dissolved species along the transects at different times

The calculated dominant aqueous compositions, including sodium, calcium, magnesium, sulphate and chloride ions, are presented along the two transects in Figure 5 and are compared to the measured concentrations. There is globally a good agreement between calculated and measured concentrations after 100 years.

It is found that the concentrations of Cl and Na increase progressively with travel distance.

They reach values of 0.08 mol kg -1 H2O and 0.09 mol kg -1 H2O in T1 and value of 0.055 mol kg -1 H2O and 0.062 mol kg -1 H2O in T2, respectively. The increased concentrations are observed in cells containing the mineral halite. Therefore, the interaction between halite and groundwater has a strong influence on the spatiotemporal evolution of dissolved sodium chloride. Assuming the dissolution of halite, the geochemical process can be expressed by halite (Eq. 13), gypsum (Eq. 14) and anhydrite (Eq. 15).

Halite:

NaCl ž$ → Na J ¡$ + Cl ¡$ (13)
The dissolved concentrations calculated in T1 and T2 after 30, 50 and 100 years show a slight decrease over time. This can be explained by a decrease in mineral mass relative to its reactive surface [START_REF] Kieffer | An experimental study of the reactive surface area of the Fontainebleau sandstone as a function of porosity, permeability and fluid flow rate[END_REF][START_REF] Colon | Experimental investigation of the effect of dissolution on sandstone permeability, porosity, and reactive surface area[END_REF][START_REF] Prikryl | Mineral dissolution in porous media: An experimental and modeling study kinetics, porosity and surface area evolution[END_REF] during the dissolution process of halite. In addition, in comparing the two transects, the Cl and Na concentrations observed in T2 are lower, due to the presence of halite only in the third zone, which is very short compared to the overall length of zones 2, 3 and 4 in transect T1.

However, the calculated SO4 and Ca concentrations generally increase with travel distance (Fig. 5). The most logical source of these ions is the dissolution of gypsum and anhydrite according to equations ( 14) and ( 15) (Jones et al., 1999;[START_REF] Wu | Hydrogeochemical characteristics of surface water and groundwater in the karst basin, southwest China[END_REF]. The gypsum and anhydrite minerals both release one mole of SO4 and one mole of Ca. However, their temporal evolution is different. The SO4 and Ca concentrations calculated at T1 and T2 do not vary between 50 and 100 years. The Ca release is very likely due to the precipitation of calcite and/or dolomite:

Gypsum: CaSO £ 2 H `O$ ž$ → Ca `J ¡$ + SO £ ` ¡$ + 2 H `O, (14)
Anhydrite:

CaSO £ ž$ + 2 H `O → Ca `J ¡$ + SO £ ` ¡$ + 2 H `O. (15)
Furthermore, the spatial distribution of the calculated Mg concentration can be characterized as follows: a decrease in concentration from the first zone to the fourth zone of T1, while the concentrations are more or less the same along transect T2. The decrease in concentrations as a function of travel distance is likely due to the precipitation of dolomite, which can occur along the flow path of T1.

Numerical modelling of the saturation indexes

The dissolution and precipitation of carbonate minerals and evaporites are determined by the saturation level of the water solution. The results of the simulation of these processes are shown in Figure 6. Dissolution of halite, anhydrite, gypsum and calcite minerals is very high in the first zones of each transect.

The computed saturation index (SI) of halite is lower than zero everywhere, which shows that the salt rock is not saturated. However, the saturation index increases at the downstream boundary because of the washing out of the dissolved quantity of halite along the flow path.

The SI values of anhydrite are negative or close to zero, which suggests that this mineral is in equilibrium with water, especially in the downstream zones of the modelled transects.

However, the saturation index of gypsum increases gradually until it reaches 0.23 and 0.21 in T1 and T2, respectively. The positive values indicate that a saturated state of this mineral is achieved where the ions accumulate in the water solution.

The saturation indexes computed for calcite and dolomite are in the range of -0.14 and -0.06 and 0.08 and 0.20 for T1 and in the range of -0.16-and -0.08 and 0.12 and 0.2 for T2, respectively. This clearly reflects the conversion from under-saturation of calcite to the oversaturation of dolomite. The lixiviation of the calcite releases the moles of Ca, which are then involved in the precipitation of dolomite. Calcite would thus be scarce, as water could not reach the necessary saturation. One may conclude that the most important reactions that can occur in the groundwater of the Sidi Hani aquifer are the dissolution of halite, anhydrite and calcite, the dissolution and precipitation of gypsum and the precipitation of dolomite. Water may instantaneously react with halite, anhydrite and gypsum and then differently increase the concentrations of Cl, Na, and SO4 in the water, which leads to an increase in the total dissolved solids. In addition, in zones where anhydrite is present, its dissolution could provoke precipitation of gypsum.

Spatiotemporal evolution of permeability

The final hydraulic conductivity profile in both modelled transects T1 and T2 is shown in Figure 7. There is a change in permeability after 100 years of interaction between the water and minerals observed in the two transects. They varied between 2.9×10 -4 and 3.7×10 -4 m s -1 and between 9×10 -5 and 11.4×10 -5 m y -1 for T1 and T2, respectively. Globally, there is an increase of the hydraulic conductivity. In the first zones of each transect this increase is due to the high dissolution of anhydrite and calcite minerals, then a decrease occurs due to dolomite precipitation. Further, downstream of zone 2, the dissolution of halite leads to an increase of hydraulic conductivity.

Model calibration

Different numerical runs were performed to quantify the numerical values of the reactive surface and the kinetic rate of mineral dissolutions and to test their influence on the chemical composition of the groundwater. The best fit of the measured chemical compositions was obtained for the adjusted dissolution rates (Table 6). For the comparison, we used the numerical results after 100 years. The model provides a good estimate of chemical elements for the tow modelled transects T1 and T2 with a RMSE value of approximately 0.01 for all of the chemical species. Note that the measured concentrations of Cl, Na, Ca, SO4 and Mg concentrations for T1 vary between 0.03 and 0.07, 0.04 and 0.07, 0.007 and 0.023, 0.019 and 0.029, and 0.015 and 0.011 mol kg -1 H2O, respectively.

For transect T2, the measured concentrations of Cl, Na, Ca, SO4 and Mg vary between 0.03 and 0.04, 0.04 and 0.05, 0.007 and 0.020, 0.020 and 0.025, and 0.015 and 0.012 mol kg -1 H2O, respectively. The numerical results were close to the field data. The obtained correlation coefficients show an acceptable correlation between measured and simulated chemical concentrations for T1 and T2 (Table 7).

Sensitivity studies

Water-rock interaction is the major factor responsible for groundwater salinity. Its influence was first analysed using two numerical studies that analysed (1) the variation of volume sensitivity of the simulated results of both cases was tested by varying their initial optimised values and taking as reference the previous numerical results obtained at t=100 years. The results of these two sensitivity studies are presented and compared to the reference results in Figure 8.

In the first sensitivity study, we increased (case 1a) and decreased (case 1b) the volume fraction of minerals up to 10%. It was shown that the dissolved ions, such as chlorides, calcium, magnesium and sulphates, are not sensitive to this parameter. If sodium is a little bit more sensitive, its response to the volume fraction variation remains lower than 10%. In the second sensitivity case, the chemical composition of the injected water flux was replaced by the chemical composition of the first observation well located upstream of each transect. Here also, the calculated dissolved concentrations of major elements do not differ significantly from those of the reference case and highlight the finding that the composition of the water at the upstream boundary is not a sensitive model parameter. As the results of case 3a and case 3b present the same spatial distribution of dissolved major elements one may conclude that the residence time is not a key parameter of influence.

However, by comparing the results of cases 3c and 3d, especially in the case of transect T1, there are significant differences between modelled concentrations of Na and Cl downstream of zone 2. Furthermore, it is clearly shown that all modelled dissolved major elements exhibit large differences in concentration when increasing (case 3b) or decreasing (case 3c) the reference Darcy velocity by a factor of 2. This underlines that the modelled fate of dissolved species is rather sensitive to the Darcy velocity. But the concentrations of the reference case are located between these extreme cases.

Conclusions

The Ouled Chamekh plain located in central-eastern Tunisia is a general example of a semiarid region where the problem of groundwater salinization is becoming more and more severe. The shallow unconfined aquifer of Sidi El Hani is the main supply source of water.

Although a large amount of data has been obtained in previous studies, these studies did not highlight the main factors responsible for the origin of the groundwater salinization in Sidi El

Hani. Therefore, further research was needed to identify the causes of groundwater salinization at this site in a more integrated manner for application to the sustainable management of groundwater resources.

The present study addresses the quantification and prediction of the chemical changes of dissolved species in the groundwater, as well as the main geochemical processes that evolve in the aquifer along a flow path. The rock-water interactions and mass transport along two transects of the aquifer were simulated using the numerical hydrogeological code KIRMAT.

The calibrated reactive transport model underscores the observation along the two transects:

the concentrations of dissolved species such as chlorides, sodium, sulphates and calcium are highly increased from the middle of each transect compared to the concentrations at the inflow edge. However, in accordance with the field data, the calculated concentrations of magnesium decrease along the flow path. These changes are clearly the result of the waterrock interaction expressed by dissolution and precipitation. One may conclude that the groundwater composition along the first 3 km of each transect is more specifically controlled by the dissolution of halite, anhydrite, gypsum, calcite and dolomite. These processes are also active in the remaining distance to the outflow edge of both transects, while gypsum changes from under-saturated to over-saturated due to the accumulation of SO4 ions. In addition, the groundwater in these zones becomes saturated with respect to anhydrite. The results of the modelled feedback effect of reactive mineral surfaces highlighted that the permeability of the porous medium is globally increased.

The present study provides useful insights into the origin of salinization and the source of the major elements in the groundwater in the Sidi El Hani basin: it specifies the origin of each chemical element and the current state of its flux in the aquifer. The numerical studies using the dynamic reactive transport model KIRMAT complete the actual knowledge of the functioning of the aquifer and provide answers concerning the transport process involved.

Figure 1 :

 1 Figure 1: Location map of the study area.

Figure 2 :

 2 Figure 2: Geological map and cross-section of the study area.

Figure 3 :

 3 Figure 3: Hydrochemical conceptual model: a) contour map of hydraulic heads of the Sidi El Hani aquifer with flow directions illustrated by red arrows (from Dridi et al., 2013) and b) spatial discretization of the T1 and T2 transects.

  constant. If the water sample is completely saturated with the dissolved mineral, SI equals zero. Positive values of SI indicate super-saturation under which the mineral would tend to precipitate; negative values indicate under-saturation and mineral dissolution.

Figure 4 :

 4 Figure 4: Distribution of Cl, Na, SO4, Ca, Mg and HCO3 ions measured in the groundwater along (a) transect 1 (T1) and (b) transect 2 (T2).

Figure 5 :

 5 Figure 5: Spatiotemporal evolution of selected dissolved species along transects T1 and T2:

Figure 6 :

 6 Figure 6: Numerically calculated saturation index (SI) as a function of travel distance and

Figure 7 :

 7 Figure 7: Numerically calculated hydraulic conductivity distribution in (a) transect 1 and (b) transect 2.

Figure 8 :

 8 Figure 8: Sensitivity study of the influence of the volume fraction of minerals (cases 1a and 1b) and the chemical composition of injected water flux (case 2) on modelled dissolved major elements of the ground water: (a) Cl, (b) Na, (c) SO4, (d) Ca, and (e) Mg

Figure 9 :

 9 Figure 9: Sensitivity study of the influence of the Darcy velocity ( cases 3a, 3b, 3c and 3d) on modelled dissolved major elements of the ground water: (a) Cl, (b) Na, (c) SO4, (d) Ca, and (e) Mg.

Table 1 :

 1 Physical and chemical compositions of water chosen at the upstream boundary (UB)

	and downstream boundary (DB) of T1 and T2, compared to the initial conditions (IC) in
			2015.
	Parameter		T1			T2
		UB	IC	DB	UB	IC	DB
	pH	7.0	8.4		7.0	8.0
	T(°C)	25.0	25.0		25.0	25.0
	Element		Concentration (mmol/kgH2O)
			Free flux			Free flux
	Ca	7	16		7	15
	Mg	15	12		15	13
	Na	40	59		40	39
	K	1	1		1	1
	Cl	30	46		30	31
	S	30	30		30	24
	C	3	3		3	5

Table 2 :

 2 Mineralogical composition and corresponding volume fraction of the Sidi El Hani reservoir.

	Mineral	Structural formula	Volume fraction (%)
	Quartz	SiO2	56
	Gypsum	Ca SO4.2H2O	10
	Halite	NaCl	4
	Anhydrite	CaSO4	3
	Dolomite	Ca Mg (HCO3)2	3
	Calcite	CaCO3	14
	Illite	Si3.5Al2.3Mg0.25O10(OH)2K0.6	10

Table 3 :

 3 Length and area of reactive surfaces of the minerals.

	Mineral	Long sides of the mineral	Reactive surface area
		(m)	(m 2 .kg H2O -1 )
	Quartz	10 -3	8.566
	Gypsum	10 -3	1.000
	Halite	10 -3	0.856
	Anhydrite	10 -3	0.428
	Dolomite	10 -3	0.571
	Calcite	10 -3	0.713
	Illite	10 -6	0.576

Table 4 :

 4 Parameters of the modelled primary and secondary minerals.

	Mineral	Structural formula	Log K (25°C)
	Illite	Si3.5Al2.3Mg0.25O10(OH)2K0.6	-43.192
	Gypsum	CaSO4.2H2O	-4.630
	Anhydrite	CaSO4	-4.390
	Halite	NaCl	1.570
	Dolomite	CaMg(CO3)2	-18.28
	Calcite	CaCO3	-8.370
	Quartz	SiO2	-3.999

Table 5 :

 5 Kinetic constants of dissolution reactions of primary and secondary minerals at 25°C.

	Mineral	Acidic conditions	Neutral conditions	Basic conditions
		a log k	b E	n c	a log k	b E	n a log k	b E	c n
	Illite	-11,31	65,4	0,77 -13,18	22,20		-17,05 17,90 -0,47
	Gypsum	-	-	-	-3,19	14,30		-	-	-
	Anhydrite	-	-	-	-6,50	31,30 -	-	-	-
	Halite	-	-	-	-0,21	7,40	-	-	-	-
	Dolomite	-3,19	36,1	0,50	-7,53	52,20 --5,11 34,80	0,50
	Calcite	-0,30	14,40	1,00	-5,81	23,50 --3,48 35,40	1,00

a: rate constant k computed at 25°C; b: Arrhenius activation energy E (kJ mole -1 ); c: reaction order n with respect to H + .

Table 6 :

 6 Initial and adjusted reactive dissolution rates.

	Mineral	Dissolution rate
		Log k (

mol years -1 kgH2O -1 ) Adjusted dissolution rate Log k (mol years -1 kgH2O -1 )

	Gypsum	-2.79	-1.92
	Anhydrite	-3.19	-1.30
	Halite	-0.21	-1.69
	Calcite	-3.00	-1.92
	Dolomite	-3.11	-3.15

Table 7 :

 7 Correlation coefficients.

		Transect 1	Transect 2
	Cl	0.86	0.77
	Na	0.70	0.77
	Ca	0.64	0.99
	Mg	0.81	0.83
	SO4	0.93	0.63

Acknowledgements

The authors sincerely thank the CRDA of the Mahdia for giving free access to the plain. We gratefully thank the farmers for their help during field work.