<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRGEOS,Unicode,biblatex,published]{cedram}
\addbibresource{crgeos20241013.bib}
\newenvironment{noXML}{}{}
\let\citep\parencite
\let\citet\textcite
\def\xcitealp#1#2{\citeauthor{#1}, \citelink{#1}{#2}}
\def\citealt#1#2{\citeauthor{#1}, \citelink{#1}{#2}}
\newcommand*{\citelink}[2]{\hyperlink{cite.\therefsection @#1}{#2}}
\def\defcitealias#1#2{}
\let\citetalias\textcite
\let\citepalias\parencite
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.92\linewidth}{#1}}
\skip\footins=13pt
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\ndash{\text{--}}
\def\hyphen{\text{-}}
\def\myRef{\ref}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\def\0{\phantom{0}}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\RequirePackage{etoolbox}
\usepackage{upgreek}
\def\jobid{crgeos20241013}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\def\xmorerows#1#2{#2}
\newcounter{runlevel}
\usepackage{multirow}
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\let\ubreak\break
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\def\back#1{}
\newcommand*{\hyperlinkcite}[1]{\hyper@link{cite}{cite.#1}}
\DOI{10.5802/crgeos.286}
\datereceived{2024-12-05}
\daterevised{2025-02-13}
\datererevised{2025-03-12}
\dateaccepted{2025-03-17}
\ItHasTeXPublished
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc}
\def\CDRdoi{10.5802/crgeos.286}
\def\citelink#1#2{\citeyear{#1}}
\def\xcitealp#1#2{\citealp{#1}}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
\def\hyperlinkcite#1#2{\citeyear{#1}}
}
\makeatother

\usepackage{upgreek}

\dateposted{2025-06-16}
\begin{document}

\begin{noXML}

\TopicFR{Hydrologie, hydrog\'ologie}
\TopicEN{Hydrology, hydrogeology}

\CDRsetmeta{articletype}{research-article}

\title{Reactive transport modelling to analyse groundwater flow path in
the young oceanic hotspot volcano of Piton des Neiges, Reunion Island
(France)}

\alttitle{Mod\'{e}lisation du transport r\'{e}actif pour analyser
l'\'{e}coulement souterrain d'un jeune volcan de point chaud
oc\'{e}anique, le Piton des Neiges, \^{I}le de La R\'{e}union (France)}

\author{\firstname{Jean-Lambert} \lastname{Join}\CDRorcid{0000-0002-2321-5341}\IsCorresp}
\address{Universit\'{e} de La R\'{e}union, Laboratoire G\'{e}oSciences
R\'{e}union, 97744 Saint-Denis, Reunion, France}
\address{Universit\'{e} Paris-Cit\'{e}, Institut de Physique du Globe
de Paris, CNRS, UMR 7154, 75005 Paris, France}
\email[J.-L. Join]{jean-lambert.join@univ-reunion.fr}

\author{\firstname{Julien} \lastname{Bonnier}}
\address{Office de l'Eau de la R\'{e}union, 33 Impasse Henri
Guillaumet, 97411 Saint-Paul, R\'{e}union, France}

\author{\firstname{Alain} \lastname{Cl\'{e}ment}}
\address{Institut Terre et Environnement Strasbourg (ITES),
Universit\'{e} de Strasbourg et CNRS, 5 rue Descartes, 67084
Strasbourg, France}

\author{\firstname{Bertrand} \lastname{Fritz}\CDRorcid{0000-0003-4896-1700}}
\addressSameAs{4}{Institut Terre et Environnement Strasbourg (ITES),
Universit\'{e} de Strasbourg et CNRS, 5 rue Descartes, 67084
Strasbourg, France}

\shortrunauthors

\keywords{\kwd{Groundwater}\kwd{Volcanic terrains}\kwd{Reunion
island}\kwd{Hydrochemical modelling}\kwd{Erorun critical zone
observatory}\kwd{Kirmat}\kwd{Vadose and phreatic weathering processes}}

\altkeywords{\kwd{Eaux souterraines}\kwd{Terrains
volcaniques}\kwd{\^{I}le de La R\'{e}union}\kwd{Mod\'{e}lisation
hydrochimique}\kwd{Observatoire de zone critique
Erorun}\kwd{Kirmat}\kwd{Processus d'alt\'{e}ration dans les zones
vadoses et phr\'{e}atiques}}

\thanks{University of Reunion Island, CNRS/INSU France, Office de
l'Eau R\'{e}union}

\begin{abstract} 
This study examines groundwater flow paths and geochemical processes in
the Piton des Neiges volcano on Reunion Island. Using data from the
Salazie Amont tunnel, bored at a depth of 1000 m under the
\textit{plan\`{e}ze} of Plaine des Foug\`{e}res, we characterise
groundwater composition and weathering processes in volcanic terrains.
The study area features high permeability in shallow lava flows,
enabling deep water infiltration, but also demonstrates a decrease in
permeability with increasing depth due to mineral precipitation of
secondary phases, mainly clays.

Reactive transport in the vadose and phreatic zones was modelled with
KIRMAT, a geochemical code simulating water--rock interactions. Initial
rainwater undergoes chemical changes while percolating through volcanic
rock, increasing in pH and mineral content. In the phreatic zone,
simulation results suggested magmatic CO$_{2}$ input, impacting water
pH and mineral composition. The model captured observed mineral trends
effectively, predicting significant secondary mineral formation at
depth, particularly clay minerals. Results illustrate how
mineralisation reduces permeability in volcanic formations over time,
progressively creating an impermeable barrier in the volcano's core.
This research provides insights into the long-term evolution of
groundwater systems in shield volcanoes, contributing to the
understanding of aquifer dynamics and material flow in volcanic
environments.
\end{abstract} 

\begin{altabstract} 
Nous \'{e}tudions la circulation des eaux souterraines dans le volcan
oc\'{e}anique du Piton des Neiges, \`{a} La R\'{e}union, ennous
appuyant sur la mod\'{e}lisation du transport r\'{e}actif pour
comprendre l'\'{e}volution chimique des eaux. Les venues d'eau
souterraines \'{e}chantillonn\'{e}es au cours du creusement du tunnel
{\og} Salazie Amont {\fg}, r\'{e}alis\'{e} \`{a} 1000 m de profondeur
sous la plan\`{e}ze de La Plaine des Foug\`{e}res, nous permettent
d'\'{e}tudier l'acquisition de la min\'{e}ralisation des eaux
d'infiltration et les processus d'alt\'{e}ration min\'{e}ralogique au
sein du massif volcanique. Deux zones hydrog\'{e}ologiques sont
distingu\'{e}es : la zone vadose, avec une infiltration verticale
rapide des eaux de pluie, et la zone phr\'{e}atique avec une
\'{e}volution chimique tr\`{e}s diff\'{e}rente marqu\'{e}e par un
enrichissement net en sodium aux d\'{e}pens du magn\'{e}sium et
calcium.

En utilisant des donn\'{e}es d'analyse et le mod\`{e}le KIRMAT, nous
simulons les transformations chimiques des eaux dans ces zones. Les
r\'{e}sultats montrent une augmentation progressive de la
min\'{e}ralisation et du pH, li\'{e}e \`{a} l'alt\'{e}ration des
basaltes. Dans la zone phr\'{e}atique, l'hypoth\`{e}se d'un apport en
CO$_{2}$ magmatique est test\'{e}e pour expliquer la variabilit\'{e}
des \'{e}chantillons d'eau. Ce mod\`{e}le permet aussi de quantifier
les min\'{e}raux secondaires, principalement des argiles,
pr\'{e}cipitant en profondeur, contribuant \`{a} la r\'{e}duction de la
perm\'{e}abilit\'{e}. Finalement ces travaux mettent en lumi\`{e}re la
dynamique d'imperm\'{e}abilisation du volcan bouclier.
\end{altabstract} 

\maketitle

\twocolumngrid

\end{noXML}

\defcitealias{Sanqueretal2010}{ibid.}

\section{Introduction} \label{sec1}\vspace*{-.3pc}

\looseness=-1
Young or active basaltic oceanic volcanoes form high islands mainly
consisting of stacks of lava flows with high hydraulic conductivities
that promote deep infiltration of surface waters, first into shallow
aquifers and later deep into the basaltic bedrock basal aquifer.
Consequently, and despite very high rainfall rates, surface water
resources are scarce. Moreover, prospecting inland groundwaters is a
challenge due to the extreme depth (hundreds of metres) of the water
table. The analysis of deep volcanic aquifers is therefore an important
recent field of research. Numerous studies which have been conducted on
other islands and volcanic archipelagos emphasise the difficulty in
defining a general representation of the circulation and water--rock
interaction conditions of deep groundwater in such terrains
\citep{CruzSilva2001, Custodio1989, Ecker1976, Izquierdo2014,
Izukaetal2010, LauMink2006, Macdonaldetal1983, PetersonMoore1987}.

These various analyses highlight how difficult it is to define a
general representation of the deep groundwater circulation conditions
in such media. 

In Reunion Island, \citet{Joinetal1997,
Joinetal2005,Joinetal2016} have
defined a general framework to explain the continuity of the water
table close to the sea level in the coastal area towards the higher
altitude water table inland. This involves lower permeability at the
core of the island. \citet{Bretetal2003} described the presence of
secondary mineralisation (clays and zeolites) that should contribute to
reducing initial volcanic hydraulic conductivity and contributed to
this conceptual model, which remained to be improved by direct
investigations. Pursuing that approach, we benefited from the Salazie
Amont tunnel project, a water diversion gallery drilled inside the
Piton des Neiges volcano at an average depth of 1000~m. Bored in the
early 2000s, this tunnel provided a unique opportunity to observe the
groundwater flow path and to analyse the water composition in the
deepest part of the volcano: the springs occurring at various levels in
the massif's cliff walls were surveyed and sampled. The chemical
contents of this spring water are used as testimony to the chemical
processes which are active during gradual infiltration.\looseness=-1

Using flow and chemical data from the tunnel's seepage waters, the
purpose of this paper is to characterise the groundwater flow paths and
to document the weathering processes resulting in the genesis of clays
as secondary minerals of alteration processes within the volcanic
terrains down to the tunnel level, in coherence with the significant
evolution of the water's chemical composition. Hypotheses for
groundwater patterns in volcanic terrains are improved by chemical
modelling. The results are interpreted with regard to the structural
and geological knowledge of the flank of this volcano.

\section{Study area} \label{sec2}
\subsection{Site description} \label{sec2.1}

Located in the south-western Indian Ocean
(21{\textdegree}S/55{\textdegree}E), east of Madagascar, Reunion Island
has the climate of a high humid tropical island, mainly defined by a
hot rainfall season from December to April and a fairly temperate dry
season from May to December. Some of the rainfall values recorded are
among the highest in the world (up to 14,000~mm${\cdot}$yr$^{-1}$ and
270~mm${\cdot}$h$^{-1}$). From a geomorphological perspective, Reunion
fits typically in to the category of high volcanic islands. It is
formed of two shield volcanoes (Figure~\ref{fig1}). Covering the
northwestern two thirds of the island, Piton des Neiges (3069~m) has
been extinct for 20,000~years \citep{GillotNativel1982}, while the
still-active volcano of Piton de la Fournaise has erupted approximately
once a year since 1998. Both present high relief marked by deep
erosional \textit{cirques} with steep cliff walls known as
\textit{remparts}.

\begin{figure*}[p!]
\includegraphics{fig01}
\caption{\label{fig1}Location map and sampling. Sampling points are
identified by a code and a symbol according to their type and location.
Plxx ${=}$ precipitations, Frxxx ${=}$ boreholes, Scxx ${=}$ springs, Rvxx
${=}$ stream or river, Tuxx ${=}$ tunnel.}
\end{figure*}

Located in the northern part of Piton des Neiges, the study area
consists of the Rivi\`{e}re des Pluies watershed, which comprises two
preserved flanks of Piton des Neiges surrounding the
proto-\textit{cirque} named Grand \'{E}boulis, drained by an entrenched
torrential river. This exceptional landscape, marked by a one thousand
meter-high cliff \textit{rempart}, provides a unique opportunity to
observe the interior of a volcano from the summit of the cliff
corresponding to the preserved structural surface, known as a
\textit{plan\`{e}ze}, down to the oldest terrains of Piton des Neiges,
around 1000~m below. At the cliff bottom, a water tunnel has been
drilled, revealing significant water inflows with physico-chemical and
isotopic characteristics at depths rarely investigated in such volcanic
environments.\vspace*{-.35pc}

\subsection{Geology and mineralogy} \label{sec2.2}

The lithology of the Rivi\`{e}re des Pluies watershed is linked to the
volcanic formations of Piton des Neiges \citep{UptonWadsworth1972,
BillardVincent1974, Kluska1997, Salvanyetal2012}. The authors agree in
defining two very different major periods of activity, reflected in
both the petrology and geomorphology of the massif. Its typical shield
volcano period 
\citep[phases~1 and 2 according to][]{BillardVincent1974},
from 2.1 to 0.43~Ma, consisted
mainly of olivine basalts and oceanic lava flows. After a period of
volcanic quiescence from 0.430 to~0.30~Ma, dominated by alteration and
erosion of the massif, volcanic activity shifted to a
stratovolcano-building stage 
\citep[phases~3 and~4 in][]{BillardVincent1974}, with the emission of
differentiated feldspar-rich lavas (0.34 to 0.029~Ma). Olivine
phenocrysts predominate in the shield volcano period, while
plagioclases becomes the dominant stage in the stratovolcano series
\citep{Valer2016}. Plagioclases are clearly more calcic in the
differentiated lava series, reaching An88 composition
\citep{Smietana2011} in particularly plagioclase-rich lavas.

Secondary mineral phases have been identified both in the shield
volcano series and soils in the differentiated stratovolcano series.
They confirm efficient weathering processes from surface to depth.
Smectites and occasionally vermiculites were observed. Microscopically,
these clays appear as a fine band (${&lt;}$10~${\upmu}$m) lining the
vacuoles in the basaltic lava. The oldest oceanites sampled in the
hydraulic gallery contain numerous zeolites, which are massive, fibrous
or lamellar and tend to fill the voids in the vacuoles. These include
phillipsite, chabazite, thomsonite and mesolite 
\citep{Bret2002,Bretetal2003}.\vspace*{-.35pc}

\subsection{Hydrogeological setting} \label{sec2.3}
\subsubsection{Water balance}

At the top of the \textit{plan\`{e}ze} above the tunnel, an
experimental watershed of the Ravine M\`{e}re Canal is part of the
Rivi\`{e}re des Pluies watershed covering an area of 0.65~km$^{2}$
upstream from the hydrometric station. The site was instrumented with
two meteorological stations, enabling water-balance components to be
calculated. Average annual rainfall is 4207~mm, of which 20\% is
evapotranspiration, 38\% runoff, for a deduced infiltration of 42\%,
i.e., 1768~mm or 56l/s per km$^{2}$ \citep{Bonnier2010}.

\subsubsection{Groundwater conceptual model}

The stacking of lava flows and associated scoriaceous products presents
high hydraulic conductivities favoring rapid and deep infiltration, as
described for active shield volcanoes. However, these hydrodynamic
characteristics display spatial variability (e.g., variations in lava
flow thickness, or presence of dikes or pyroclastic ash products), and
there is still little information available for estimating hydrodynamic
properties on the scale of a massif.

In the high elevation part of the Piton des Neiges volcano, groundwater
flow is evidenced by springs that mark aquifer outlets and appear in
deeply incised valleys. A hydrogeological conceptual model has been
proposed \citep{Joinetal2016, Joinetal2005}, describing three types of
aquifers. (1)~Superficial groundwaters within the perhydrated andosols
occur in the cloud forest area at high elevation (above 1300~m~Ngr).
(2)~Perched groundwaters refer to the intercepted part of deep
infiltration due to local impervious horizons such as palaeosols and
paleovalleys. Springs of that type occur mainly in the cliff walls of
erosion canyons. (3)~Finally, a deep saturated zone extends above the
oldest volcanic units, and defines the phreatic zone at the origin of
the permanent flow of the entrenched valley of Rivi\`{e}re des Pluies
located approximately 1000~m below the cloud forest area. On the basis
of this diagram, we can assume that the groundwater flows sampled in
the bored tunnel at the bottom of the valley allow us to characterise
the basal aquifer.

The conceptual model can be reduced to two layers: a vadose zone,
comprising superficial and perched groundwater, and a phreatic zone
corresponding to the base water table. The vadose zone in Reunion
Island is of the exceptional thickness of several hundred metres, which
is due to the extreme permeability of the most recent volcanic
formations. Vertical infiltration predominates, with temporary surface
or perched water tables. The phreatic zone, on the other hand,
represents a saturated zone associated with the base water intercepted
by the tunnel and tapped by drilling in coastal areas.

During the tunnel drilling operations, fluid pressure measurements of
more than 30~bars were recorded at the bottom of the phreatic zone in a
borehole drilled horizontally as the tunnel progressed. These
measurements enabled us to locate the top of the water table at an
altitude of around ${+}$900~m~Ngr. As the average ground level above is
around ${+}$1500~m~Ngr, we can deduce a thickness of 600~m for the vadose
zone.

Several water samples in the catchment area and the tunnel were dated
using the tritium isotope method. The results enabled us to refine this
conceptual model. They provide residence time for \mbox{water} taken
from the vadose zone (surface and perched water tables) and the
phreatic zone (tunnel and borehole). \citet{Sanqueretal2010} suggest
ages ranging from 75 years at the top of the phreatic zone to 150 years
for the oldest water samples dated in the tunnel.

\section{Materials and methods} \label{sec3}
\subsection{Water sampling and analysis} \label{sec3.1}

An initial campaign of sampling and chemical analysis was carried out
during the excavation of the water transfer tunnel (2006--2014). Waters
from the tunnel, springs and boreholes corresponding to the vadose and
phreatic zones were sampled simultaneously. This initial dataset was
complemented by chemical analyses carried out as part of the
hydrological monitoring of the Rivi\`{e}re des Pluies watershed within
France's research infrastructure of \mbox{OZCAR} (Critical Zone Observatories:
Research and Application) and its observatory ERORUN (Observatory of
Erosion in Reunion) dedicated to observing the water cycle and critical
zone \citep{Toulieretal2024}.

From the complete dataset, 419 analyses were selected on the basis of a
quality criterion defined by a percent charge-balance error (\%CBE) of
less than 5\% (Supplementary Table~1).

The samples analysed were grouped into five categories according to the
type of sampling carried out, distinguishing between borehole waters
(Frxxx) taken from the water table downstream from the catchment area,
rainwaters (PlPF1340) taken from a rain gauge located at 1340~m
altitude in the cloud forest area, surface water taken from the
Rivi\`{e}re des Pluies or its tributaries on the massif (Rvxxxx),
tunnel waters (Tuxxx) taken from various depths in the water transfer
tunnel, and finally spring waters (Scxxx) taken from various points on
the massif and corresponding to superficial or perched aquifers.

\subsection{Hydrochemical modelling} \label{sec3.2}
\subsubsection{The modelling strategy} \label{sec3.2.1}

Reunion Island provides a most interesting field for studying
water--rock interaction processes, with extremely high rainfall on
basaltic rocks known for being very reactive to alteration in the
weathering process. We thought that the modelling approach might be
useful in order to better understand how the composition of the
infiltrated rainwater undergoes major change, to reach the composition
of well waters. We know from the analytical data on samples (cf.\ 
Supplementary Table~1) that the resulting solutions are still diluted
in subsurface aquifers first, but become more concentrated in deep
aquifers as shown by samples collected in the tunnel, with pH values as
high as 9 to~10.

To approach the problem we first applied a modelling strategy with two
major steps according to the hydrogeological conceptual model defined
in Figure~\ref{fig2}: (1)~a \textit{vadose} episode with the first
reaction of rainwater with basaltic rocks at near-surface CO$_{2}$
conditions, (2)~a subsequent \textit{phreatic} episode with proceeding
alteration under variable $p$CO$_{2}$, no longer controlled by the
atmosphere, with a possible sharp increase in the pH values due to the
loss of the pH buffering effect of CO$_{2}$. In our modelling approach
the results of the first simulation led us to define a more complex
modelling process, as described later, with a third episode at depth
with the probable influence of some deep CO$_{2}$.

\begin{figure*}
\includegraphics{fig02}
{\vspace*{.25pc}}
\caption{\label{fig2}Hydrogeological settings. Hydrogeological
conceptual model of inland aquifers and water pathways from rainwater
to the tunnel waters. Representation of the vadose zone
(${=}$superficial and perched groundwater) and phreatic zone (${=}$basal
groundwater).}
\end{figure*}

\subsubsection{Modelling the reactive transport} \label{sec3.2.2}

We used the geochemical code KIRMAT (KInetics of Reaction and Mass
Transfer) first described by \citet{Gerardetal1998a, Gerardetal1998b},
when it was created for simulating water--rock interaction processes in
an open 1-D system with the circulation of fluid, in order to predict
the cumulative effect of alteration with time in open cells. The
simulated primary minerals may dissolve and secondary precipitated
minerals are predicted. Among the latter, the clay minerals are
considered as solid solutions with variable compositions
\citep{TardyFritz1981}. 

The KIRMAT modelling is a one-dimensional model particularly adapted to
the simulation of the consequence of alteration in a fractured pathway.
Each calculation is made at a given temperature in the range from
0~{\textdegree}C to 300~{\textdegree}C. Input data describe the initial
composition of the reactive cells, including rock and aqueous solution.
The aqueous solution is considered to saturate porosity. Its chemical
composition is described in terms of major and minor chemical elements
with indication of Eh and pH values. The $p$CO$_{2}$ partial pressure may
be used as a fixed parameter. The rock characteristics include total
porosity, permeability (cf.\ Section~\ref{sec3.2.3}), reactive surface
area, and mineral composition with kinetic dissolution rate constants
of minerals (cf.\ Section~\ref{sec3.2.4}).

Primary minerals are dissolved as a function of dissolution rates
calculated along the process following Equation~(\ref{eq1}):
{\begin{equation}\label{eq1}
r_{\mathrm{d}} = k_{\mathrm{d}}S_{\mathrm{effm}}\ 
an\mathrm{H}^{+} \ 
[1-((Q_{\mathrm{m}}/K_{\mathrm{m}})na)nb]
\end{equation}}\unskip
where $r_{\mathrm{d}}$ is the dissolution rate of a mineral m
(mol${\cdot}$yr$^{-1}{\cdot}$kg$^{-1}{\cdot}$w); $k_{\mathrm{d}}$ is
the dissolution rate constant (mol${\cdot}$m$^{-2}{\cdot}$yr$^{-1}$;
Table~\ref{tab1}); $S_{\mathrm{effm}}$ is the reactive surface area of
the mineral (m$^2{\cdot}$kg$^{-1}{\cdot}$w); $a$H$^{+}$ is the activity
of protons; $n$ is an experimental power, where $na$ and $nb$ are
experimental exponents depending on the pH of the solution;
$Q_{\mathrm{m}}$ is the ion activity product of the mineral for a given
composition of the aqueous solution; and $K_{\mathrm{m}}$ is the
equilibrium constant for the hydrolysis reaction of the mineral at the
given temperature and pressure. Kinetic reaction rates $k_{\mathrm{d}}$
are given in Table~\ref{tab1} with three different values depending on
the pH range in the solution.

\begin{table*}%tab1
\caption{\label{tab1}Dissolution kinetic rates $k_{\mathrm{d}}$ (moles/kg/yr)
from \citep{Made1991}}
\tabcolsep=8pt
\begin{tabular}{lccc}
\thead
&amp; Kda &amp; Kdn &amp; Kdb\\
\endthead
\multicolumn{4}{l}{$T=19$~{\textdegree}C}\\
Forsterite &amp; $1.907 \times 10^{-0}$ &amp; $1.875\times 10^{-3}$ &amp; $7.447\times 10^{-6}$\\
Fayalite &amp; $2.374\times 10^{-1}$ &amp; $2.377\times 10^{-4}$ &amp; $9.376\times 10^{-7}$\\
Anorthite &amp; $2.112\times 10^{-2}$ &amp; $6.762\times 10^{-4}$ &amp; $2.711\times 10^{-6}$\\
Albite &amp; $2.124\times 10^{-2}$ &amp; $7.244\times 10^{-4}$ &amp; $2.844\times 10^{-6}$\\
Mugearite glass &amp; $1.355\times 10^{-3}$ &amp; $1.352\times 10^{-5}$ &amp; $5.408\times 10^{-8}$\vspace*{.35pc}\\
\multicolumn{4}{l}{$T=21$~{\textdegree}C}\\
Forsterite &amp; $2.109\times 10^{-0}$ &amp; $2.042\times 10^{-3}$ &amp; $8.128\times 10^{-6}$\\
Fayalite &amp; $2.656\times 10^{-1}$ &amp; $2.600\times 10^{-4}$ &amp; $1.023\times 10^{-6}$\\
Anorthite &amp; $2.394\times 10^{-2}$ &amp; $7.800\times 10^{-4}$ &amp; \0$3.00\times 10^{-6}$\\
Albite &amp; $2.438\times 10^{-2}$ &amp; $8.128\times 10^{-4}$ &amp; $3.162\times 10^{-6}$\\
Aphyric glass &amp; $1.596\times 10^{-3}$ &amp; $1.622\times 10^{-5}$ &amp; $6.456\times 10^{-8}$\vspace*{.35pc}\\
pH limits for Kds &amp; pH ${&lt;}$ 6 &amp; 6 ${&lt;}$ pH ${&lt;}$ 7 &amp; pH ${&gt;}$ 7
\botline
\end{tabular}
{\vspace*{-.45pc}}
\end{table*}

For validation of this kinetic approach, the \mbox{KIRMAT} code has
been progressively improved and applied to different types of
alteration systems: weathering processes \citep{Ackereretal2018,
Schaffhauseretal2014}, geothermal studies \citep{Lucasetal2017,
Lucasetal2020}, nuclear waste storage \citep{Martyetal2010,
Ngoetal2014} and recently the simulation of possible hydrogen
generation in geothermal systems \citep{Murrayetal2020}.

The Code KIRMAT does not use the reactive surface area values as
generally derived from BET measurements which lead to very high values
of surface areas (about tens of m$^{2}$/g of rock). These values
correspond to all forms of porosity including microporosity almost, not
connected to circulated fluids. In KIRMAT we refer to the rock surface
in contact with a given amount of circulating fluid in macroporosity
and we estimate the total surface ($S_{\mathrm{t}}$) area of minerals
in contact with the fluid in small fractures or cracks using a
geometrical model of fissures. The surface area corresponding to each
rock-forming mineral ($S_{\mathrm{eff}}$), used in the kinetic rate law
for dissolution, is then calculated using the volumic fraction (Fv) of
this mineral in the basaltic rock ($S_{\mathrm{eff}} = S_{\mathrm{t}}
\cdot F_{\mathrm{v}}$). In this study we used two different values for
the VGW modelling runs ($S_{\mathrm{t}} =  10~{\mathrm{m}}^2$/kg
H$_2$O) and PGW runs ($S_{\mathrm{t}} = 100~\mathrm{m}^2$/kg H$_{2}$O)
respectively.

The formation of secondary minerals produced in the process is
simulated at thermodynamic \mbox{equilibrium}. 

The dissolution equations for both primary and secondary minerals are
presented in Table~2 in the Supplementary Data with the solubility
products.

The vertical evolution of the water chemistry in the volcanic massif
has been modelled using a succession of boxes intersecting two
successive different lithological zones (Figure~\ref{fig3}): 
\begin{enumerate}[(i)]
\item First box representing the 600~m-thick vadose zone made up of
differentiated mugearite volcanic flows (30 simulated cells of 20~m in
length). 
\item Second box representing the 300~m-thick phreatic zone located in
basaltic volcanic flows corresponding to the ancient phase of the
oceanite series (15 simulated cells of 20~m in length).
\end{enumerate}

\begin{figure*}
\includegraphics{fig03}
\caption{\label{fig3}Chemical evolution of major ion components.
(a)~Piper diagram for all sampled water; (b)~Schoeller diagram for
representative waters of rainwaters (PLPFxx), spring waters (SCxx) and
tunnel waters (Tuxx).}
\end{figure*}

The circulation of the fluid flow between the successive cells is
defined by a given effective velocity (Darcy's velocity/porosity).

A feedback effect of the calculated dissolution/precipitation processes
on porosity and permeability in the reactive cells is included
\citep{Lucasetal2017}. Concerning the control of CO$_2$ along the
simulated process, KIRMAT offers two options as indicated in Table~7:
using a selected fixed value for a given calculation run (see runs PGW2
a and b, and VGW) or a given initial value allowed to evolve as a
response to the chemical evolution of the system (see runs PGW1~a
and~b).

For the purposes of this study, KIRMAT has been adapted to store not
only the calculations but also all the chemical data necessary to apply
the code DIAGRAMMES \citep{Simler2021} in order to produce the Piper
and Schoeller diagrams for the solutions.

The kinetics of dissolution of primary minerals in the reactive cells
is controlled by kinetic reaction rates given in Table~\ref{tab1}
\citep{Made1991} and in Complementary Material (Annexe Table~2).

The parameters governing fluid flow and solute transport differ between
the vadose and phreatic zones, and have therefore been modelled
sequentially along two distinct pathways to reflect the contrasting
petrological characteristics.\ The vadose zone, primarily composed of
mugearite and plagioclase-rich lavas, differs from the phreatic zone,
which lies within the basaltic shield volcano phase. Additionally,
$p$CO$_{2}$ is assumed to remain close to \mbox{atmospheric} levels in
the vadose zone, whereas it may decrease in the phreatic zone.

\subsubsection{Estimated hydrodynamical parameters} \label{sec3.2.3}

In order to simulate the mineralisation processes from infiltration of
rainwater (Figure~\ref{fig2}: point PL1340) through 600~m towards a
300-m thick phreatic zone (Figure~\ref{fig2}: point Tu3000), we use
hydrodynamic properties that could represent the global structure of
the volcano defined by two main volcanic units (a shield volcano
substratum overlaid by the recent stratovolcano). These two units have
been assimilated to hydrogeological layers (Figure~\ref{fig2}). Their
permeability and porosity parameters are defined on the basis of
previous hydrogeological modelling work carried out in Reunion Island.
Values of hydraulic conductivity and drainage porosity range
respectively from $1\times 10^{-9}$ to $1.5\times 10^{-4}$~m/s and
0.1\% to 15\% \citep{Belleetal2018, Folio2001, Marechaletal2014}.

Using KIRMAT, effective velocity is of major importance. For that
purpose, we used carbon-14 and tritium dating of the groundwaters
sampled from the streams, springs and tunnel waters within the
Rivi\`{e}re des Pluies watershed \citep{Sanqueretal2010}.

\paragraph{The vadose zone}

The time taken for the water to go through the vadose zone is assessed
at 77 years using isotopic time dating \citepalias{Sanqueretal2010} of the
top of the basal groundwater (from 75 to 82 years), i.e., an effective
velocity of $600/77=7.79$~m/yr or a Darcy velocity of 0.779~m/yr
for a porosity of 10\%.

\paragraph{The phreatic zone}

Dating of the tunnel waters has provided ages ranging from 77 years to
a maximum of 147 years. The youngest waters correspond to significant
water inflows in draining sectors marked by fractures and numerous
dikes, and we can assume that the water sampled in these sectors
characterises the top of the water table. In contrast, the oldest water
is found in the least drained areas and corresponds, conversely, to
sampling close to the wall of the water table, whose maximum thickness
has been estimated at around 300~m. These observations enable us to
estimate the vertical transfer time between the roof and the water
table wall as $147-77 = 70$ years. From this we can deduce the
effective velocity for the vertical component of the flow in the
phreatic zone, i.e., a Darcy velocity of 0.429~m/yr considering a
porosity of 10\% in the phreatic zone (300~m and 70~y). These values
are consistent with the hydraulic conductivity estimates used to
simulate the tunnel's drainage capacity \citep{Marechaletal2014}.

\subsubsection{Geochemical conditions and parameters} \label{sec3.2.4}
\paragraph{Initial composition of rocks}

To characterise the average mineralogy of the volcanic terrains in
these two environments, we used the mineralogical compositions
collected and analysed by \citet{Hoareauetal2005} as part of his
experimental work on leaching volcanic rocks from Reunion Island. We
choose two aphyric-textured samples whose composition was deemed
representative of the two modelized volcanic series (Table~\ref{tab2}).

\begin{table}%tab2
\caption{\label{tab2}Mineralogical composition of the basalts and
volcanic glasses used in the simulations (in moles)}
\tabcolsep9pt
\begin{tabular}{ccc}
\thead
Simulation &amp; Vadose &amp; Phreatic\\
\endthead
Lava &amp; Mugearite &amp; Aphyric basalt\\
An/Ab &amp; 1.5 &amp; 1.5\\
Forsterite &amp; 0.0643 &amp; 0.1444\\
Fayalite &amp; 0.1112 &amp; 0.1283\\
Albite &amp; 0.1469 &amp; 0.1136\\
Anorthite &amp; 0.2204 &amp; 0.1703\\
Glass &amp; 0.4572 &amp; 0.4434\vspace*{.3pc}\\
Total &amp; 1 &amp; 1\vspace*{.6pc}\\
\multicolumn{3}{l}{{Composition of the glass (atoms per 2O)}}\\
Si &amp; 0.8892 &amp; 0.805\0\\
K &amp; 0.1542 &amp; 0.0558\\
Na &amp; 0.2398 &amp; 0.0487\\
Ca &amp; 0.0247 &amp; 0.3778
\botline
\end{tabular}
\tabnote{$^*$The composition of the glass has been calculated from the
total composition of the basalt (in \% oxides) and the An/Ab ratio ${=}$
1.5. The solubility product K of the glasses has been calculated from
the solubility products of oxides using an ideal solid solution mode.}
\end{table}

\paragraph{Initial composition of waters (Table~\myRef{tab3})}

The water entering the vadose zone has the composition of typical
rainwater collected on the volcano (sample reference PlPF1300 of 17
July 2019).

\begin{table*}%tab3
\caption{\label{tab3}Major ion chemistry (mg/L) of entering water in
the modelling process for the vadose (No~3) and phreatic pathways
(No~459) (a)~cations, (b)~anions}
\begin{tabular}{ccccccccccc}
\tbody
\multicolumn{11}{l}{(a)}\\
$n$ &amp; Reference &amp; Name &amp; Date &amp; $T$~({\textdegree}C) &amp; pH &amp; c25~{\textdegree}C &amp; Ca &amp; Mg &amp; Na &amp; K\\
3 &amp; PlPF1300 &amp; PlPF1300 &amp; 17 July 2019 &amp; 23.9 &amp; 6.3 &amp; 15.8 &amp; \00.38 &amp; 0.21 &amp; \01.53 &amp; 0.15\\
459 &amp; 38714 &amp; Phreatic pathway &amp;&amp; 21 &amp; 8.2 &amp; 183 &amp; 16.40 &amp; 5.65 &amp; 11.90 &amp; 0.25\vspace*{.3pc}\\
\multicolumn{11}{l}{(b)}\\
$n$ &amp; Reference &amp; Name &amp; Date &amp; HCO$_{3}$ &amp; CO$_{3}$ &amp; Cl &amp; SO$_{4}$ &amp; NO$_{3}$ &amp; SiO$_{2}$ &amp;\\
3 &amp; PlPF1300 &amp; PlPF1300 &amp; 17 July 2019 &amp; \01.23 &amp; 0 &amp; 2,54 &amp; 0.65 &amp; 0 &amp; ${&lt;}$DL &amp;\\
459 &amp; 38714 &amp; Phreatic pathway &amp;&amp; 97.7\0 &amp; 2.70 &amp; 4.62 &amp; 7.07 &amp; 0 &amp; 24 &amp;
\botline
\end{tabular}
\end{table*}

The water entering the phreatic zone has the composition of the
simulated water leaving the overlying vadose zone.

When the initial composition of the water in rock (infiltrated
rainwater) was defined, the results of preliminary modelling indicated
that its precise composition does not exert a significant impact on the
final concentration profile during simulation times that are longer
than the water residence time in the column (77 years in the vadose
zone and 70 years in the phreatic zone).

\paragraph{Temperature}

The mean temperature of spring water sampled in the vadose zone was
19.14~{\textdegree}C, with a standard deviation of 1.43. Typical
phreatic water samples from the tunnel had a mean temperature of
21.26~{\textdegree}C, with a standard deviation of 2.37. Accordingly,
we chose a temperature of 19~{\textdegree}C to simulate the vadose
zone, and 21~{\textdegree}C for the phreatic zone.
\vspace*{-.4pc}

\section{Results and discussion} \label{sec4}
\vspace*{-.15pc}
\subsection{Concentrations of major elements} \label{sec4.1}
\vspace*{-.15pc}
\subsubsection{Chemical characteristics of the aqueous{\hfill\break}
solutions sampled in the field} \label{sec4.1.1}

The evolution of the chemical composition of rainwater from Plaine des
Foug\`{e}res to the transfer tunnel clearly shows two distinctly
different pathways in the vadose zone and the phreatic zone. The first
pathway, named Vadose Geochemical Way (VGW), concerns the chemical
evolution of rainwater as it moves towards the perched groundwater and
the top of the water table (springs and boreholes). The second, named
Phreatic Geochemical Way (PGW), concerns the marked evolution of
samples taken at different depths in the tunnel towards water at the
base of the water table. All the water samples taken are represented in
a Piper diagram (Figure~\ref{fig3}a).

The VGW pathway corresponds to an increase in concentration of the four
major cations (Ca, Mg, Na, K) together with alkalinity
($\mathrm{HCO}_{3}^{-}$, $2\mathrm{CO}_{3}^{2-}$) as a direct
consequence of alteration processes of the mineral components of the
basalts. The PGW solutions, from the deepest part of the phreatic zone
in the tunnel itself, become very alkaline (Mg$^{++}$, Na$^{+}$,
$\mathrm{HCO}_{3}^{-}$, $\mathrm{CO}_{3}^{2-}$).

We have selected nine representative samples from each type of water
(rainwater, borehole water, tunnel spring water) in order to detail the
global mineralisation with a Schoeller diagram (Figure~\ref{fig3}b).

This diagram clearly discriminates the rainwater group (PlPF1300) from
spring waters (ScApol and ScCH) and groundwaters sampled at different
depths in the tunnel (TU4925, Tu2935, Tu3279). It highlights the
general increase in mineralisation between rainwater and groundwater.
It also shows the significant chemical evolution between the vadose
zone and the phreatic zone, marked mainly by an increase in Na at the
expense of Ca and Mg.

It should be noted that this geochemical trend goes together with a
continuous increase in pH, reaching 8 to 9 for VGW and 9 to 10 for PGW.
This general evolution in pH is typical of alteration processes at low
temperature (around 20~{\textdegree}C) when the gas CO$_{2}$ is
available for providing protons:
{\begin{equation}\label{eq1a}
\mathrm{CO}_{2} + \mathrm{H}_{2}\mathrm{O} = \mathrm{H}^{+} +
\mathrm{HCO}_{3}^{-}
\end{equation}}\unskip
The alteration process takes advantage of these protons for the
dissolution of silicates in a general scheme of alteration:
{\fontsize{9.7}{12}\selectfont $$
\mathrm{silicates} + \mathrm{H}^{+} + \mathrm{H}_{2}\mathrm{O} +
\mathrm{CO}_{2} = \mathrm{cations} + \mathrm{silica} +
\mathrm{alkalinity}
$$}\unskip
The equilibrium condition for Equation~(\ref{eq1a}) relates the
activity of H$^{+}$ to the $p$CO$_{2}$:
{\begin{eqnarray}
\log \mathrm{K}(\mathrm{CO}_{2}) &amp;=&amp;
\log(\mathrm{H}^{+}) + \log
(\mathrm{HCO}_{3}^{-})
\nonumber\\&amp;&amp; 
-\, \log (p\mathrm{CO}_{2}) - \log (\mathrm{H}_{2}\mathrm{O})
\label{eq2}
\end{eqnarray}}\unskip
as $\mathrm{pH} = -\mathrm{log}(\mathrm{H}^{+})$
{\begin{equation}\label{eq3}
\mathrm{pH} = \log(\mathrm{HCO}_{3}^{-}) - \log(p\mathrm{CO}_{2})
-\log(\mathrm{H}_{2}\mathrm{O}) -\log\mathrm{K} (\mathrm{CO}_{2})
\end{equation}}\unskip
with (${\ldots}$) referring to activities.

Equation~(\ref{eq3}) shows clearly that increased alkalinity goes
together with increased pH at a given $p$CO$_{2}$ and also that the pH
will increase even further if the $p$CO$_{2}$ diminishes as expected with
depth in the phreatic zone.

\subsubsection{Observed secondary mineralogy sampled in the tunnel}
\label{sec4.1.2}

Rock samples collected during drilling operations were analysed, with
four samples examined using a Bruker X-ray diffractometer to
characterise the secondary mineralogy. In the deepest section of the
tunnel, between PM 3065 and 3209, clay minerals were predominant,
specifically identified as smectite-type montmorillonite across all
four samples. These findings align with observations in other tunnel
sections, indicating that the abundance of zeolites seen in outcrops of
older volcanic formations diminishes towards the volcano's core, where
clay minerals are more prevalent \citep{Bretetal2003}.

\subsection{Hydrochemical functioning} \label{sec4.2}
\subsubsection{Simulation results} \label{sec4.2.1}

The quality test of the composition of water samples collected in the
field was conducted using the numerical code KIRMAT. In this necessary
first step in any modelling approach to water--rock interactions, the
calculations determined the chemical speciation of all chemical
elements in the aqueous solutions. As mentioned in
Section~\ref{sec3.1}, this makes it possible to begin by checking the
electrical neutrality of the solution as an important quality index of
the analytical data set. For the discussion and comparison with
simulated values, we only selected representative samples shown in the
Schoeller{\break} diagram.

\begin{figure*}
\includegraphics{fig04}
\caption{\label{fig4}Chemical evolution of major ion components for
simulated versus sampled waters described with a Piper diagram and a
Chadha diagram (all samples) and a Schoeller diagram (selected samples
and simulated waters). (a)~Piper diagram for all sampled waters;
(b)~Schoeller diagram for selected waters. Red crosses correspond to
increments of calculation of simulated vadose pathway (Increment
38,000) and phreatic pathway (Increment 4500); (c)~Chadha diagram for
all sampled waters.}
\vspace*{-.05pc}
\end{figure*}

Figure~\ref{fig4} presents the results of our simulation based on the
proposed modelling approach.\ In the initial simulation stage, starting
with rainwater (PlPF1300) interacting with mugearite rocks, we
observed the evolution of the water composition along the vadose
zone under constant $p$CO$_{2}$ conditions. Following this, the final
water composition from the vadose zone simulation served as the input
for a second simulation step, which models basalt weathering under
phreatic conditions with variable $p$CO$_{2}$. The Piper diagram
(Figure~\ref{fig4}a) shows that the simulated water composition closely
aligns with the measured data, a result further corroborated by the
Chadha diagram \citep{Chadha1999} (Figure~\ref{fig4}b). The Schoeller
diagram confirms that the KIRMAT model effectively predicts the
formation of spring water in the vadose zone and tunnel water in the
phreatic\break zone.

\looseness=-1
In contrast, in the simulation of the phreatic zone with variable
$p$CO$_{2}$, the less evolved samples from the tunnel are not as well
simulated. The rapid increase in pH is accompanied by a marked drop in
Mg concentration, which becomes underestimated at the start of the
simulation compared to the lowest mineralized waters collected in the
tunnel. Sensitivity tests revealed that the system is strongly
influenced by the partial pressure of CO$_{2}$. Assuming complete
consumption of dissolved CO$_{2}$ in the phreatic zone would result in
a total loss of Mg and very high pH values (${\gg}$10), not observed in
the most basic \mbox{tunnel} waters.

The issue of CO$_{2}$ control (or not) in the alteration process is a
geochemical key point, as it will regulate/control the pH. In the upper
part of the volcanic massif (VGW), the $p$CO$_{2}$ may remain not far
from the atmospheric value ($10^{-3.5}$) and this will be considered in
our simulation conditions. However, the alteration process will
progressively cut this relation with near-surface conditions as
solutions move and infiltrate the deeper part of the volcanic{\break}
massif. 

It is clear that deep-waters collected in the phreatic zone (PGW) show
a high pH value (9~to~10).\ This is consistent with a \'{e} in
$p$CO$_{2}$ far from atmospheric conditions. 

However, the possible saturation with respect to calcite may introduce
the well-known pH buffer of calcite. 

We used the modelling approach to further understand what happens at
depth with CO$_{2}$. To better represent the diversity of the waters
sampled, we propose the hypothesis of a magmatic CO$_{2}$ input
\citep{Benardetal2023a},
thus maintaining a minimal level of $p$CO$_{2}$ in the fractures.

We suggest modelling a second circulation pathway in order to account
for fractured zones subject to potential magmatic CO$_{2}$ diffusion.
Two models of water evolution in the phreatic zone result in the
\mbox{following}:

Firstly, a PGW1 path, corresponding to the least draining zones of the
aquifer (2.5\% to 10\% porosity), where $p$CO$_{2}$ is consumed from its
initial value in the vadose zone and drops to very low values, reaching
as low as $10^{-7}$. With this option, a very low value of $p$CO$_{2}$
($10^{-7}$) is reached after under 70 years (54 years) because we
limited the decrease to that value. This evolution of the $p$CO$_{2}$
goes together with a substantial increase in the pH, which reaches
values of between 9 and 10, and even more if the simulation is forced
to 70 years, which was ultimately unrealistic if we compared the final
solutions of the simulations to the composition of samples collected in
the\break tunnel.

For that reason, we tested an alternative PGW2 pathway, representing
groundwater circulation in draining fractures (10\% or 2.5\% porosity)
connected to a CO$_{2}$ supply of deep origin, enabling the initial
$p$CO$_{2}$ of the vadose zone to be maintained. Such a CO$_2$ supply
helps to control the $p$CO$_2$ in equilibrium with the circulation of
solutions at depth at a level not very different from the near surface
conditions (${-}4 &lt; \log p\mathrm{CO}_2 &lt; {-}3$) instead of decreasing
$p$CO$_2$ to very low values like that obtained in the modelling with
variable CO$_2$ option ($\log p\mathrm{CO}_2 = {-}7$) which goes together
with a very high pH value up to 10 and more. This may happen in some
zones of the massif not concerned by this deep CO$_2$ effect but is not
realistic for most of the circulation pathways concerned by our study
as shown by the data for waters samples collected in the tunnel
(Table~\ref{tab4}). This is the reason why we tested both options for
CO$_2$ (fixed under deep CO$_2$ influence or variable.

\begin{table*}%tab4
\caption{\label{tab4}pH and $p$CO$_{2}$ in the samples from the tunnel\vspace*{3pt}}
\tabcolsep3pt\fontsize{9.3}{12}\selectfont
\begin{tabular}{ccccccccc}
\tbody
Sample no &amp; 300 &amp; 302 &amp; 307 &amp; 308 &amp; 312 &amp; 320 &amp; 325 &amp; 327\\
Location &amp; TU 8751 &amp; TU 4932 &amp; TU 9075 &amp; TU 1238 &amp; TU 2935 &amp; TU 5087 &amp; TU 3000 &amp; TU 1238\\
Date &amp; \parbox[t]{1.4cm}{\centering 09~June 2017} &amp; 
\parbox[t]{1.7cm}{\centering 10~October 2010} &amp; 
\parbox[t]{1.1cm}{\centering 08 May 2021} &amp; 
\parbox[t]{1.7cm}{\centering 04 October 2019} &amp; 
\parbox[t]{1.4cm}{\centering 08 March 2019} &amp; 
\parbox[t]{1.2cm}{\centering 10 April 2007} &amp; 
\parbox[t]{1.7cm}{\centering 05 October 2021} &amp; 
\parbox[t]{1.6cm}{\centering 05 October 2021}\vspace*{2.5pt}\\
$t$~({\textdegree}C) &amp; 18.9 &amp; 25.1 &amp; 20.2 &amp; 19.7 &amp; 23,9 &amp; 24.1 &amp; 24.5 &amp; 19.7\\
Ionic balance &amp; 3\% &amp; 1\% &amp; 3\% &amp; 2\% &amp; 2\% &amp; 0\% &amp; 2\% &amp; 4\%\\
pH &amp; 8.48 &amp; 9.5 &amp; 8 &amp; 8.6 &amp; 9.68 &amp; 8.6 &amp; 10.02 &amp; 8.85\\
log $p$CO$_{2}$ &amp; ${-}$3.763 &amp; ${-}$4.965 &amp; ${-}$3.135 &amp; ${-}$3.912 &amp; ${-}$5.145 &amp; ${-}$3.983 &amp; ${-}$5.525 &amp; ${-}$4.186\vspace*{1pt}
\botline
\end{tabular}
{\vspace*{.6pc}}
\end{table*}

The relation between pH and $p$CO$_{2}$ in the samples from the tunnel is
documented in Table~\ref{tab4} for eight samples covering the pH range
8 to~10. Their $p$CO$_{2}$ equilibrium varies from $10^{-3.14}$ to
$10^{-5.5}$, which supports our hypothesis of the presence of some deep
source of CO$_{2}$ in the system, as samples with $p$CO$_{2}$ values as
low as $10^{-7}$ were never found. With the two options PGW1 and PGW2,
we simulated a kind of minimum and maximum effect of that possible
CO$_{2}$ diffusion at depth.

The results of the PGW1 and PGW2 simulations not only enable us to
better describe the variable chemical composition of the various waters
collected in the tunnel, but also suggest the hypothesis of \mbox{continuous}
diffusion of deep magmatic CO$_{2}$ associated with double porosity in
fractures of volcano-tectonic origin. This potential contribution of
deep magmatic CO$_{2}$ is well documented in the two massifs of Reunion
Island, based on observations of the soils and thermo-mineral waters of
the two volcanic edifices \citep{Liuzzoetal2015, Boudoireetal2023,
Benardetal2023b, Benardetal2024}.

\subsubsection{Evolution of minerals in the simulated zones}
\label{sec4.2.2}

\begin{table*}%tab5
\caption{\label{tab5}Composition of clay produced in the simulation}
\begin{tabular}{lcccccccc}
\thead
Simulation &amp; Porosity (\%) &amp; Si &amp; Al &amp; Fe (III) &amp; Mg &amp; O &amp; (OH) &amp; Ca\\
\endthead
VGW &amp; 10 &amp; 3.371 &amp; 2.25 &amp; 0.126 &amp; 0.38 &amp; 10 &amp; 2 &amp; 0.315\\
PGW1a &amp; 10 &amp; 3.01 &amp; 1.656 &amp; 0.188 &amp; 1.72 &amp; 10 &amp; 2 &amp; 0.495\\
PGW1b &amp; 2.50 &amp; 3.01 &amp; 1.656 &amp; 0.188 &amp; 1.72 &amp; 10 &amp; 2 &amp; 0.495\\
PGW2a &amp; 10 &amp; 3.407 &amp; 1.76 &amp; 0.146 &amp; 1.031 &amp; 10 &amp; 2 &amp; 0.297\\
PGW2b &amp; 2.50 &amp; 3.407 &amp; 1.76 &amp; 0.146 &amp; 1.031 &amp; 10 &amp; 2 &amp; 0.297
\botline
\end{tabular}
{\vspace*{.2pc}}
\end{table*}

Each primary mineral in the rock is dissolved acccording to its kinetic
dissolution rate. Of the secondary minerals predicted by the
simulation, the majority are due to a clay phase whose detailed
composition in the different simulations is given in
\mbox{Table}~\ref{tab5}. These compositions were predicted in the last
cell of the simulated domains for the maximum duration. Both phases are
simulated as ideal solid solutions. The kinetics of dissolution of
primary minerals in the reactive cells is controlled by kinetic
reaction rates given in Table~\ref{tab1}.

As shown in Table~\ref{tab5} for run PGW1 or PGW2, the simulations made
with two different values of porosity did not show differences in terms
of clay composition in the simulated cells after the selected duration:
the chemical evolution of the solutions is directly related to the
interaction duration and not to porosity. However, in the next section
we will see that the consequence in terms of amounts of minerals
produced can be observed in Table~\ref{tab6} total mass balance in the
reservoir is related to the amount of solution reacting with the
basalt, and therefore related to the porosity for a given flux of
water. Decreasing the porosity from 10\% to 2.5\% decreases the volumic
balance by a \mbox{factor} of~4.

\begin{table*}%tab6
\caption{\label{tab6}Minerals dissolved or precipitated (cm$^{3}$/kg
H$_{2}$O)}
\begin{tabular}{lccccc}
\thead
Alteration pathways &amp; PGW1a &amp; PGW1b &amp; PGW2a &amp; PGW2b &amp; VGW\\
\endthead
Temperature &amp; 21~{\textdegree}C &amp; 21~{\textdegree}C &amp; 21~{\textdegree}C &amp; 21~{\textdegree}C &amp; 19~{\textdegree}C\\
U ${=}$ DARCY global (KIRMAT) (m/y) &amp; 0.4286 &amp; 0.107 &amp; 0.4286 &amp; 0.107 &amp; 0.78\\ 
p ${=}$ porosity $-(0&lt;\mathrm{p}&lt;1)$ &amp; 0.1 &amp; 0.025 &amp; 0.1 &amp; 0.025 &amp; 0.1\\
incr ${=}$ final increment of calculation &amp; 10,786 &amp; 10,787 &amp; 2918 &amp; 2919 &amp; 38,714\\
$t$~${=}$ duration (y) &amp; 54.84 &amp; 54.77 &amp; 70 &amp; 70 &amp; 77\\
Final pH &amp; 10.91 &amp; 10.91 &amp; 8.57 &amp; 8.57 &amp; 8.27\\
Final log $p$CO$_{2}$ &amp; ${-}$7 &amp; ${-}$7 &amp; ${-}$3.26 &amp; ${-}$3.26 &amp; ${-}$3.33{\vspace*{6pt}}\\
\multicolumn{6}{l}{Dissolution (cm$^{3}$/kg H$_{2}$O)} \\
Forsterite &amp; 0.2403 &amp; 0.240 &amp; 0.3067 &amp; 0.3067 &amp; 0.01299\\
Fayalite &amp; 0.0285 &amp; 0.0284 &amp; 0.0363 &amp; 0.0363 &amp; 0.00298\\
Anorthite &amp; 0.5536 &amp; 0.553 &amp; 0.7016 &amp; 0.7016 &amp; 0.08493\\
Aphiric glass &amp; 0.0018 &amp; 0.0018 &amp; 0.0022 &amp; 0.0022 &amp; 0.000225\\
High albite &amp; 0.3719 &amp; 0.371 &amp; 0.4897 &amp; 0.4897 &amp; 0.05945\\
\textbf{Total dissolved} &amp; \textbf{1.1921} &amp; \textbf{1.1942} &amp; \textbf{1.5365} &amp; \textbf{1.5365} &amp; \textbf{0.1626}{\vspace*{6pt}}\\
\multicolumn{6}{l}{Precipitation (cm$^{3}$/kg H$_{2}$O)}\\
Clay solid solution end-members: &amp;&amp;&amp;&amp;&amp;\\
Pyrophyllite &amp; $2.79\times 10^{-4}$ &amp; $2.75\times 10^{-4}$ &amp; 0.064 &amp; 0.064 &amp; 0.03326\\
Talc &amp; $8.23\times 10^{-3}$ &amp; $8.21\times 10^{-3}$ &amp; 0.3975 &amp; 0.3975 &amp; 0.1578\\
Fe3 Pyrophyllite &amp; $2.60\times 10^{-7}$ &amp; $2.6\times 10^{-7}$ &amp; $3.27\times 10^{-5}$ &amp; $3.27\times 10^{-5}$ &amp; $1.21\times 10^{-5}$\\
Ca Muscovite &amp; $2.79\times 10^{-1}$ &amp; $2.79\times 10^{-1}$ &amp; 0.6578 &amp; 0.6578 &amp; 0.07253\\
Ca Phlogopite &amp; $5.07\times 10^{-1}$ &amp; $5.06\times 10^{-1}$ &amp; 0.2516 &amp; 0.2516 &amp; 0.002\\
Ca-Ferrimuscovite &amp; $8.46\times 10^{-2}$ &amp; $8.45\times 10^{-2}$ &amp; 0.1079 &amp; 0.1079 &amp; 0.00885\\
Total clay solid solution &amp; 0.8791 &amp; 0.8663 &amp; 1.4799 &amp; 1.4789 &amp; 0.1323\\
Calcite &amp; 0.0271 &amp; 0.0283 &amp; 0.1324 &amp; 0.1324 &amp; 0.00017{\vspace*{5pt}}\\

Phillipsite-K &amp; 0.0024 &amp; 0.0024 &amp; 0.0016 &amp; 0.0016 &amp; 0\\
Laumontite &amp; 0.2676 &amp; 0.2649 &amp; 0 &amp; 0 &amp; 0{\vspace*{5pt}}\\

Minerals other than clay &amp; 0.2971 &amp; 0.2956 &amp; 0.1324 &amp; 0.1324 &amp; 0.00017\\
\textbf{Total precipitated} &amp; \textbf{1.1762} &amp; \textbf{1.1654} &amp; \textbf{1.6113} &amp; \textbf{1.61131} &amp; \textbf{0.1325}{\vspace*{5pt}}\\

\textbf{Volumic balance of minerals} &amp; \textbf{${-}$0.0199} &amp; \textbf{${-}$0.0288} &amp; \textbf{0.0748} &amp; \textbf{0.0748} &amp; \textbf{${-}$0.0307}
\botline
\end{tabular}
\end{table*}

\subsubsection{Extrapolating simulation results} \label{sec4.2.3}

In KIRMAT modelling, all the calculations are made as a reference to
1~kg H$_{2}$O of solution. In order to predict the evolution of a given
reservoir in the basaltic massif, one must consider the porosity of the
rock to extrapolate the simulated mass balances to the scale of the
reservoir under consideration. This is shown in Table~\ref{tab6} for
PGW1 and PGW2 with porosity 10\% (a) and 2.5\% (b).

It is interesting to note that in all the simulations the volumic
balance between dissolved minerals and precipitated minerals remains
near equilibrium after 147 years, which corresponds to an evolution in
porosity of less than 1\%. The tendency is to increase the porosity in
the vadose zone slightly and in the first phreatic pathway
(PGW1$_{\mathrm{a},\mathrm{b}}$) and to reduce the porosity slightly in
the final one. A significant evolution of the porosity would only be
predicted for the long geological term. 

However, these evolutions correspond to significant volumes of
dissolved and produced minerals after only 147 years. A volume $V$ of a
given dissolved or precipitated mineral, calculated by KIRMAT in
cm$^{3}$/kg H$_{2}$O, can be extrapolated at the site scale
(approximately 1~km$^{3}$ of basaltic terrains), taking into account
the porosity~$n$. Therefore, the volumic balance of minerals
($V$[cm$^{3}$/kg H$_{2}$O]) corresponds in 1~km$^{3}$ to a total
dissolved or precipitated amount of mineral: $V_{\mathrm{tot}}
(\mathrm{m}^{3}) = 10^{6}~\mathrm{V}{\cdot}\mathrm{n}$.

It appears that the predicted amounts of minerals precipitated are
about 10 times higher in the phreatic phase than in the vadose phase.
In 1~km$^{3}$ of basaltic terrains, clay minerals are the main
secondary minerals and amount to 16,000~m$^{3}$ in the vadose phase and
120,000 to 160,000~m$^{3}$ in the phreatic pathway.

\section{Conclusion} \label{sec5}

{\vspace*{-2pt}}

Since 2008, first as part of the monitoring of the excavation of the
Salazie Amont tunnel and then as part of the ERORUN research
infrastructure, sampling of rainwater, springs, river boreholes and
tunnel water has been conducted from the rainwater infiltration zone at
1300~m altitude down to the deepest level of the phreatic zone reached
by the tunnel (alt.~350~m). By selecting more than 400 water analyses,
we were able to identify an initial evolution in water mineralisation
in the vadose zone, which is perfectly distinct from that observed in
the various samples taken in the phreatic zone.

To parameterize the KIRMAT water/rock interaction model, we made
realistic assumptions about the geochemical composition of the rocks
and the partial pressure of CO$_{2}$ in the vadose and phreatic zones.
In~the phreatic zone, we tested the possibility of magmatic CO$_{2}$
degassing, with results that enable us to better frame the variability
of sampled waters when they were taken from fractures or dikes.

In this context, the KIRMAT simulation tool enabled us to reconstitute
the dual evolution of rainwater flowing through 1000~m of volcanic
formations. The simulated mineralogical balances enable us to
characterise the precipitation of neo-formed clay minerals in the
phreatic zone.

Ultimately, these results make it possible to clarify the flow of
materials in underground volcanic flows and for the first time to
characterise the dynamics of progressive impermeabilization within a
shield volcano. These dynamics partly explain the decrease in
permeability of volcanic formations on a spatial and temporal scale,
which diminishes towards the inner zones of the massif and according to
the age of the volcanism.

{\vspace*{-2pt}}

\section*{Declaration of interests}

{\vspace*{-2pt}}

The authors do not work for, advise, own shares in, or receive funds
from any organisation that could benefit from this article, and have
declared no affiliations other than their research organisations.

{\vspace*{-2pt}}

\section*{Acknowledgments}

{\vspace*{-2pt}}

The ERORUN observatory funded by OSU-Reunion is a part of the \mbox{OZCAR}
research infrastructure (\url{http://www.ozcar-ri.org}), which is
financially supported by the University of Reunion Island and the
CNRS/INSU France. We thank Y.~Benoit, M.~David, G.~Lebeau, S.~Ferreira
and P.~Stamenoff for providing technical laboratory assistance.
Geochemical data and dating conducted through the Salazie Amont tunnel
project were funded by the Office de l'Eau R\'{e}union.

The authors wish to acknowledge the contribution of both reviewers for
their critical reading of the manuscript and their many useful
remarks.

We are grateful to E.~Nicolini for our helpful discussions and also the
Office de l'Eau R\'{e}union, the CVH-DEAL, M\'{e}t\'{e}o-France and the
D\'{e}partement de La R\'{e}union for participating in the ERORUN
observatory work by providing part of the Supplementary Data included
in this data paper.

C.~Northam is thanked in particular for her careful proofreading and
improvement of the English language. 

\back{}

\section*{Supplementary data}

Supporting information for this article is available on the journal's
website under \printDOI\ or from the author.

\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{crgeos-286-suppl.pdf}}
\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{crgeos-286-suppl.xlsx}}

\printbibliography
\refinput{crgeos20241013-reference.tex}

\end{document}
</pre></body></html>