\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRGEOS,Unicode,Thematic]{cedram}
\usepackage[square,comma,authoryear]{natbib}
\renewcommand\bibsection{\section*{\refname}}
\newenvironment{noXML}{}{}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{}%\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.94\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\ndash{\text{--}}
\def\hyphen{\text{-}}
\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{crgeos20240525}
%\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{}
\def\rmDelta{\Delta}
\DOI{10.5802/crgeos.277}
\datereceived{2024-06-28}
\daterevised{2024-10-21}
\dateaccepted{2024-10-29}
\ItHasTeXPublished
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc}
\def\CDRdoi{10.5802/crgeos.277}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
}
\makeatother

\usepackage{upgreek}

\begin{document}

\begin{noXML}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Geophysics for an in-depth understanding of the Critical Zone}
%\def\cdr@specialtitle@french{La g\'eophysique pour comprendre la zone critique en profondeur}
%\makeatother

\CDRsetmeta{articletype}{research-article}

\title{Defining the sensitivity of cosmic ray muons to groundwater
storage changes}

\alttitle{Caract\'erisation de la sensibilit\'e des muons cosmiques aux
variations dans le stockage des eaux souterraines}

\author{\firstname{Mat{\'\i}as} \lastname{Tramontini}\CDRorcid{0000-0003-4635-5010}\IsCorresp}
\address{CONICET - Facultad de Ciencias Astron\'omicas y Geof\'isicas,
Universidad Nacional de La Plata, Argentina}
\email[M. Tramontini]{mtramontini@fcaglp.unlp.edu.ar}

\author{\firstname{Marina} \lastname{Rosas-Carbajal}\CDRorcid{0000-0002-5393-0389}}
\address{Universit\'e de Paris, Institut de Physique du Globe de Paris,
CNRS, UMR 7154, F-75238 Paris, France}

\author{\firstname{Fabio Iv\'an} \lastname{Zyserman}\CDRorcid{0000-0002-2534-6521}}
\addressSameAs{1}{CONICET - Facultad de Ciencias Astron\'omicas y
Geof\'isicas, Universidad Nacional de La Plata, Argentina}

\author{\firstname{Laurent} \lastname{~Longuevergne}\CDRorcid{0000-0003-3169-743X}}
\address{Geosciences Rennes -- UMR 6118, Univ. Rennes, CNRS, 35000
Rennes, France}

\author{\firstname{Christophe} \lastname{Nussbaum}\CDRorcid{0000-0002-8763-2320}}
\address{Swiss Geological Survey at Swisstopo, Seftigenstrasse 264,
CH-3084 Wabern, Switzerland}

\author{\firstname{Jacques} \lastname{Marteau}\CDRorcid{0000-0003-1941-6842}}
\address{Institut de Physique des 2 Infinis de Lyon, UMR 5822,
CNRS-IN2P3, Universit\'e de Lyon, Universit\'e Claude Bernard Lyon 1,
France}

\shortrunauthors

\keywords{\kwd{Muography}\kwd{Groundwater variations}\kwd{Hydrogeophysics}}

\altkeywords{\kwd{Muographie}\kwd{Eaux souterraines}\kwd{Hydrog\'eophysique}}

\thanks{Swisstopo}

\begin{abstract} 
Quantifying groundwater storage variations is key for sustainable water
resource management, estimating droughts and climate change impacts,
and the environmental protection of the Critical Zone. Despite the many
geophysical methods traditionally used to indirectly infer groundwater
storage, there is still a gap in the spatial and temporal scales
covered. Furthermore, the groundwater mass quantification is usually
done through proxies as few methods characterize a physical property
directly linked to the volumetric water content. In this work we employ
for the first time, the muography method to investigate groundwater
variations at the scale of 100 m. We study a mountainous karstic
aquifer in the Mont Terri Underground Rock Laboratory (URL),
Switzerland. We use a muon detector to compute a muon image of average
density, and to study muon rate variations over time. After correcting
these variations for atmospheric effects, we interpret the muon time
series in terms of groundwater storage changes. Comparisons with river
stream flow, and volumetric soil water content data, show similar
variation patterns to the groundwater estimates from our muon
measurements. Since we only observe significant changes in a particular
region of the whole volume scanned by our detector, we analyze in
detail the role of the geometrical setup in the sensitivity of
muography to groundwater variations. We find that the geometrical
relationship between the trajectories of detected muons and the
topography has a significant effect in the amplitude of the detected
variations. This finding points to a potential muon-based groundwater
storage monitoring with optimized experimental setups to amplify the
groundwater storage signal.
\end{abstract}

\begin{altabstract} 
La quantification des variations du stockage des eaux souterraines est
essentielle pour la gestion durable des ressources en eau, l'estimation
des s\'echeresses et des impacts du changement climatique, et la
protection environnementale de la zone critique. Malgr\'e les
nombreuses m\'ethodes g\'eophysiques traditionnellement utilis\'ees
pour d\'eduire indirectement le stockage des eaux souterraines, il
existe encore des lacunes dans les \'echelles spatiales et temporelles
couvertes. En outre, la quantification de la masse d'eau souterraine
est g\'en\'eralement r\'ealis\'ee \`a l'aide d'approximations, car peu
de m\'ethodes caract\'erisent une propri\'et\'e physique directement
li\'ee \`a la teneur en eau volum\'etrique. Dans ce travail, nous
utilisons pour la premi\`ere fois la muographie pour \'etudier les
variations des eaux souterraines sur une \'echelle de 100 m. Nous
\'etudions un aquif\`ere karstique montagneux dans le laboratoire
souterrain du Mont Terri (URL), en Suisse. Nous utilisons un
d\'etecteur de muons pour calculer une image muonique de la densit\'e
moyenne et pour \'etudier les variations du taux de muons dans le
temps. Apr\`es avoir corrig\'e ces variations des effets
atmosph\'eriques, nous interpr\'etons les s\'eries temporelles de muons
en termes de changements dans le stockage des eaux souterraines. Les
comparaisons avec le d\'ebit des cours d'eau et les donn\'ees
volum\'etriques sur la teneur en eau du sol sont corr\'el\'es aux
mod\`eles de variations des eaux souterraines estim\'ees par nos
mesures muographiques. \'Etant donn\'e que nous n'observons des
changements significatifs que dans une r\'egion particuli\`ere du
volume total balay\'e par notre d\'etecteur, nous analysons en d\'etail
le r\^ole de la configuration g\'eom\'etrique dans la sensibilit\'e de
la muographie aux variations des eaux souterraines. Nous constatons que
la relation g\'eom\'etrique entre les trajectoires des muons
d\'etect\'es et la topographie a un effet significatif sur l'amplitude
des variations d\'etect\'ees. Cette d\'ecouverte indique qu'il est
possible de surveiller le stockage des eaux souterraines \`a l'aide de
muons en optimisant les configurations exp\'erimentales afin
d'amplifier le signal de stockage des eaux souterraines.
\end{altabstract}

\maketitle

\twocolumngrid

\end{noXML}

\clearpage

\section{Introduction}

The external surface of the continents, extending from the outer limits
of the vegetation down to and including the zone of groundwater, has
been called ``the Critical Zone'' to highlight its crucial role for the
``availability of life-sustaining resources''
\citep{gaillardet2018ozcar}. While surface freshwater storage
components (streams, lakes, snow, glaciers) are easily observable
components of our landscapes, a major part of continental water
resources resides below ground in groundwater systems and is broadly
inaccessible to direct observations. Groundwater systems are hosted in
thick and complex heterogeneous aquifers \citep{gleeson2014glimpse}.
Water flows slowly in geological layers over depths up to kilometers,
locally intercepting the surface where it interacts with rivers or
wetlands, and therefore atmosphere and ocean
\citep{alley2002flow}.\looseness=-1

Geophysical methods are key for characterizing and monitoring the
subsurface part of hydrological systems, at spatial and temporal scales
ranging from the pore scale to kilometers, and from seconds to years
\citep{hermans2023advancing}. These are indirect methods that infer
different physical properties of the subsurface, such as electrical
resistivity \citep{mcgarr2021geophysical}, seismic velocity
\citep{blazevic2020time}, density \citep{crossley2013measurement}, 
spontaneous potential \citep{grobbe2021multi}, and hydrogen content of
pore fluid \citep{schmidt2020quantifying}, which are related to
volumetric water content. Satellite-based methods such as
Interferometry of Synthetic Aperture Radar (InSAR), and the Gravity
Recovery and Climate Experiment (GRACE) are crucial tools that can be
used to quantify groundwater storage changes across large geographical
areas. InSAR provides precise surface deformation measurements with a
precision on the order of a few mm/year across areas ranging from
100~km$^2$ to 5000~km$^2$ \citep{castellazzi2016assessing}, but its
application is limited to aquifers that undergo measurable deformation
due to poroelastic responses to changes in the water table. GRACE
measures the Earth's gravity field with high sensitivity
\citep{schmidt2008hydrological}, allowing it to detect mass changes
equivalent to a water disk about 1 centimeter thick at the land
surface, spanning a few hundred kilometers or more
\citep{longuevergne2010grace}. These techniques are effective for
monitoring groundwater fluctuations at regional and basin
scales{\break} \citep{leblanc2009basin, syed2008analysis}, though their
spatial resolution is limited and not suitable for detailed local
studies, typically on the order of hundreds of meters.

\looseness=-1
Cosmic-ray geophysical methods offer promising alternatives for soil
moisture estimation and groundwater assessment at intermediate spatial
scales.\ Cosmic-ray neutron sensing has proven effective in estimating
soil moisture content over areas spanning several hundred meters,
providing near-surface volumetric soil moisture measurements
\citep{Evans2016, Khli2015}.\ Similarly, the attenuation of underground
cosmic-ray electron flux has been proven to be sensitive to variations
in soil water content following rainfall events \citep{Taketa2022}.
However, these methods are limited to shallow depths, typically at the
centimeter scale, offering insights primarily into the uppermost layers
of soil. At this scale, there is a gap of geophysical methods that can
accurately quantify groundwater variations, that is, that can directly
relate the physical property measured to the volumetric water content
in a specific area.

Muography is a rather recent geophysical method that is directly
sensitive to the density of geological bodies. It makes use of
cosmic-ray muons, which are sub-atomic particles generated from the
interaction of cosmic rays with the Earth's atmosphere. Due to their
penetrating nature, muons can reach deep underground, making them
suitable for non-invasive imaging of geological bodies
\citep{bonechi2020atmospheric}. Muography consists of measuring the
attenuation of the natural muon flux arriving to the surface of the
Earth after it traverses matter. By detecting the flux of muons from
different directions, it is possible to reconstruct images of density
contrasts in geological bodies. A key consideration of muography is
that the detector must be positioned downhill or beneath the region of
interest, as muons arrive from above. This positioning requirement
makes tunnels and natural cavities suitable locations for deploying
muon detectors, providing opportunities to apply the method. This
imaging technique has been successfully applied to various geophysical
applications \citep{bonomi2020applications}, including imaging of the
density distribution of volcanic edifices
\citep{bajou2023high,rosas2017three}, detecting hidden cavities
\citep{morishima2017discovery, cimmino20193d}, reconstructing the
bedrock morphology underneath an \mbox{active} glacier
\citep{nishiyama2019bedrock}, mapping underground structures
\citep{borselli2022three}, mining exploration of ore bodies
\citep{beni2023transmission}, and detecting the presence of
discontinuities produced by karstic networks \citep{Lesparre2016},
among other applications \citep{cohu2023furnace}.

The continuous acquisition of muon data has shown promising potential
to monitor hydrothermal activity in volcanoes, such as La Soufri\`ere
de Guadeloupe in the Lesser Antilles, France
\citep{jourde2016monitoring,le2019abrupt}. Muography has also been
utilized for a time-lapse study of mass density changes over time
\citep{olah2023muon}, and for monitoring gravity-driven debris flows in
the Sakurajima volcano, Japan \citep{olah2021muographic}. The mentioned
time-lapse and monitoring applications suggest that Muography could be
key to bridge the current gap in groundwater storage estimations
between the regional satellite and the punctual borehole measurements
previously mentioned. Variations in groundwater content directly
translate into fluctuations in the amount of matter encountered by
muons as they reach underground detectors. This direct interplay
between subsurface mass changes and muon rate measured could result in
significant variations due to evolving water storage changes. 

Previous work has established a correlation between middle-atmosphere
temperature variations and the muon rate registered underground
\citep{osprey09}. This finding highlighted the sensitivity of muography
to atmospheric conditions and emphasized the importance of considering
such external factors in the interpretation of underground muon
measurements \citep{tramontini2019middle}, since this effect becomes
increasingly significant as the overburden opacity of the site
increases \citep{grashorn2010atmospheric}.

In this work, we use Muography to study groundwater variations in the
surroundings of the Mont Terri Underground Rock Laboratory (URL,
Switzerland). After correcting for the atmospheric temperature effect,
we analyze the muon rate variations from a series of continuous muon
measurements, and we compare the corrected muon rates with different
datasets associated to groundwater content. We propose a model to
explain the differences in muon rate sensitivity to groundwater
variations observed in the different regions scanned. We end with
discussions about the implications for future experiments and potential
uses of Muography in hydrogeophysics in general.

\section{Experimental setup} \label{sec2}

The Mont Terri Underground Rock Laboratory (URL) is a research facility
located in the Jura Mountains in Switzerland
(Figure~\ref{fig:intro:geology_scenario}). The facility has been
operational since 1996 and has been used for various research projects
related to hydrogeological, geochemical, and geotechnical
characterisation of a clay formation focused on radioactive waste and
CO$_2$ storage, including more recently heat storage 
\citep{bossart2018mont}. The geological profile and hydrogeological
classification at the Mont Terri URL is presented in
Figure~\ref{fig:intro:telescope}a. The study region contains two
distinct karstified aquifers, separated by a shaly aquitard. The local
aquifer, located to the north of the valley, presents higher inflow
rates compared to the regional aquifer located in the southern part. 
Based on the hydrogeological classification by
\citet{marschall2004mont}, the local aquifer is primarily composed of
sandy limestones that exhibit fissures and karstic features, with
inflow rates ranging from $12~\mathrm{l}{\cdot}\mathrm{s}^{-1}$ to
$200~\mathrm{l}{\cdot}\mathrm{s}^{-1}$ due to localised
water-conducting features. Conversely, the regional aquifer comprises
fissured and karstic limestones, with inflow rates exceeding
$1~\mathrm{l}{\cdot}\mathrm{s}^{-1}$. Acting as a barrier between these
aquifers, the aquitard is composed of shaly marls, limestones and marly
clay. Due to this distinct setup, the Mont Terri URL stands as an
exceptional site for conducting muography experiments within the
context of \mbox{hydrogeological} studies.

\begin{figure*}
\includegraphics{fig01}
\caption{\label{fig:intro:geology_scenario}(a)~Mont Terri's placemarks.
(b)~Aerial view of Mont Terri, with motorway tunnel and rock
laboratory. Image taken from the Mont-Terri URL website
\url{https://www.mont-terri.ch}.}
{\vspace*{1pc}}
\end{figure*}

\begin{figure*}
\includegraphics{fig02}
\caption{\label{fig:intro:telescope}(a)~Geological profile adapted from
\citet{Nussbaum2017}. Hydrogeological classification from
\citet{marschall2004mont}. Geological map based on the Geological Atlas
of Switzerland 1:25,000, map No.~40 ``St-Ursanne'', available through
\url{https://map.geo.admin.ch}. (b)~The muon detector deployed in the
Mont Terri Underground Rock Laboratory (URL).}
\end{figure*}

The muon detector used in this study was developed by the DIAPHANE
project \citep[e.g.,][]{marteau2012muons, marteau2017diaphane}. 
It consists of three detection matrices composed by $N_x=16$ and
$N_y=16$ plastic scintillators bars, in the horizontal and vertical
directions, which define $16 \times 16$ pixels of $5 \times
5~\text{cm}^2$ each (Figure~\ref{fig:intro:telescope}b). The distance
between the front and rear matrices is set to 100~cm. When a muon
passes through the three matrices of the detector in time coincidence,
an event is registered and the corresponding trajectory is
reconstructed from the sets of pixels fired in each matrix. Muon
detection is accomplished through the use of wavelength-shifting (WLS)
optical fibers that transport the photons generated by the
scintillators to the photomultipliers (PMTs) for detection. More
details on the optoelectronic chain can be found in
\citet{marteau2014implementation}. Due to the large rock volume under
study in relation to the size of the detector, we admit a point-like
approximation of the detector where events with pixels in the same
relative direction are considered to correspond to the same trajectory
\citep{lesparre2010geophysical}. This approximation yields a total of
$(2N_x-1)\times(2N_y-1)=961$ axes of observation $\hat{r}_{i,j}$, where
$(i,j)=(x,y)_{\mathrm{front}}-(x,y)_{\mathrm{rear}}$ and
$(x,y)_{\mathrm{front}}$ and $(x,y)_{\mathrm{rear}}$ represent the
$(x,y)$ pixels in the front and the rear matrices, respectively. Thus,
each $\hat{r}_{i,j}$ correspond to an unique angular direction.

\looseness=-1
In this work, we use a muon dataset acquired by the muon detector
deployed in the Mont Terri URL from October 2016 to February 2018. The
dataset encompasses 365 days of observations, accounting for data
acquisition interruptions for work at the site. The detector was
installed in the Niche IS, 220~m below the surface, with an inclination
of 65{\textdegree} and oriented 152{\textdegree}N. Due to the employed
setting, the range of length of rock traversed by the muons is of
200--500~m. Preceding the underground measurements, a calibration
experiment was conducted by measuring the open-sky muon flux at the
zenith \citep{lesparre2010geophysical}.
\vspace*{-.5pc}

\section{Methodology} \label{sec3}
\vspace*{-.2pc}
\subsection{Computing the muon image} \label{sec3.1}

The physical quantity estimated by muography is the opacity, $\varrho$
(mwe), which represents the amount of matter traversed by the muons
along their trajectories and is given by:
{\begin{equation} \label{eq:opacity}
\varrho(L) = \int_{L} \rho(l) d(l) = \langle\rho\rangle \times L,
\end{equation}}\unskip
where $\rho$ ($\text{g}{\cdot}\text{cm}^{-3}$) is the density of the
medium, $l$ is the coordinate measured along the trajectory of
longitude $L$~(m), and $\langle\rho\rangle$
($\text{g}{\cdot}\text{cm}^{-3}$) is the average density along the
trajectory. $\varrho$ is often expressed in mwe (meter water
equivalent), where 1~mwe corresponds to
$100~\mathrm{g}{\cdot}\mathrm{cm}^{-2}$.

A muon image consists of presenting the estimated values of average
density for each axis of observation of the muon detector $\langle\rho
(\hat{r}_{i,j})\rangle$. To obtain it, we first need to estimate the
integrated muon flux, $I~(\mathrm{cm}^{-2}{\cdot}\mathrm{sr}^{-1}{\cdot}
\mathrm{s}^{-1})$, which quantifies the flux of muons that have
energies above a certain threshold. The dependency of this integrated
flux with the opacity is expressed by
{\begin{equation}
I(\varrho,\theta) = \int_{E_{\mathrm{min}}(\varrho)}^{\infty}
\phi_{\mu}(E',\theta)\,\text{d}E',
\end{equation}}\unskip
where $\theta$ is the zenith angle and $E_{\mathrm{min}}(\varrho)$ is
the minimum energy required by a muon to traverse a given opacity
$\varrho$, and $\phi_{\mu}$ ($\text{cm}^{-2}{\cdot}\text{sr}^{-1}{\cdot}
\text{s}^{-1}{\cdot}\text{GeV}^{-1}$) is the differential muon flux,
which represents the flux of muons with a given energy value.
$E_{\mathrm{min}}$ is computed by solving numerically $E_{\mathrm{min}}
- \int_{0}^{\varrho} ({\text{d}E}/{\text{d}\varrho})\,\text{d}\varrho =
E_{\mu}$, where $E_{\mu}=0.10566$~GeV is the muon rest energy and
${\text{d}E}/ {\text{d}\varrho}$ accounts for the muon energy loss
through matter and is described by the attenuation constants provided
by the Particle Data Group tables \citep{PDG22}.\looseness=-1

\looseness=-1
We estimate the integrated muon flux along each axis of observation as
$I(\hat{r}_{i,j}) = ({N_{\mu}(\hat{r}_{i,j})})/({\rmDelta T \times
\mathcal{T}(\hat{r}_{i,j})})$, where $N_{\mu}(\hat{r}_{i,j})$ is the
number of muons detected along each axis of observation, $\rmDelta T (s)$
is the acquisition duration, and
$\mathcal{T}(\hat{r}_{i,j})~(\mathrm{cm}^2{\cdot}\mathrm{sr})$ is the
direction-dependent acceptance function, which quantifies the
detector's efficiency for detecting muons. $\mathcal{T}(\hat{r}_{i,j})$
is obtained from a calibration experiment at open-sky
\mbox{conditions}, where the measured muon rate is compared to the
theoretically expected integrated muon flux
$I^{\mathrm{open}\hyphen\mathrm{sky}} (\hat{r}_{i,j})$:
{\begin{equation}\label{eq:acceptance}
\mathcal{T}(\hat{r}_{i,j})=\frac{N_{\mu}^{\mathrm{calib}}(\hat{r}_{i,j})}
{\rmDelta T^{\mathrm{calib}} \times I^{\text{open-sky}}
(\hat{r}_{i,j})},
\end{equation}}\unskip
where $N_{\mu}^{\mathrm{calib}}(\hat{r}_{i,j})$ and $\rmDelta
T^{\mathrm{calib}} (s)$ are the number of muons registered along each
axis of observation during the calibration experiment, and the
calibration acquisition duration, respectively. The associated
confidence interval is determined through error propagation
\citep{bevington2003data}. To calculate
$I^{\mathrm{open}\hyphen\mathrm{sky}}(\hat{r}_{i,j})$ we follow
\citet{cohu2022atmospheric}: we use the CORSIKA code
\citep{heck1998corsika} to perform air shower simulations, and
incorporate the primary flux parameterization at the top of the
atmosphere described by \citet{horandel2003knee}.

Once the integrated muon flux is estimated, we use tabulated values of
the $(I,\varrho,\theta)$ dependence to determine the opacity
$\varrho(\hat{r}_{i,j})$. For each value of integrated flux and its
corresponding zenith angle ($I(\hat{r}_{i,j}),\theta(\hat{r}_{i,j})$)
we do a 3-D interpolation in the ($I,\varrho,\theta$) space. Finally,
we use a 2-meters horizontal resolution Digital Elevation Model to
compute the rock thickness traversed by the muons across the axes of
observation, $L(\hat{r}_{i,j})$. Then, the average density muon image
is given by $\langle \rho (\hat{r}_{i,j})\rangle=\varrho
(\hat{r}_{i,j})/L(\hat{r}_{i,j})$.

\subsection{Time series of muon rate} \label{sec3.2}

To study the influence of groundwater storage changes between the Mont
Terri URL and the surface on our muon measurements, we create time
series of the measured muon rate, that is, the number of muons measured
per unit of time. Our muon detector can discriminate 961 different axes
of observation. In order to improve the temporal resolution of the
time series, we merge the signals from adjacent directions together
when studying time variations \citep{jourde2016muon}. We compute the
average cosmic muon rate time series, $R$, using a 30-day width Hamming
moving average window \citep{hamming1998digital}. The associated
confidence interval is determined as the standard deviation of the mean
\citep{bevington2003data}. We assume that the time intervals between
consecutive muon detections follow an exponential distribution
\citep{leo1994techniques}, which is consistent with the assumption that
$R$ follows a Poisson process \citep{gaisser2016cosmic}.

Due to the significant amount of rock between the surface and the
tunnel, only relatively high-energy muons arrive to our detector.
\citet{tramontini2019middle} showed that in this case the muon rate
might be sensitive to temperature variations in the middle atmosphere
and we therefore need to account for this effect before interpreting
the time series in terms of groundwater variations. For this, we use
the effective temperature $T_{\mathrm{eff}}$~(K)
\citep{grashorn2010atmospheric}, which is a weighted average of the
atmospheric temperature profile, and the effective temperature
coefficient $\alpha_T$ \citep{barrett52}, which gives the relative
change in the muon rate $R$ given a relative change in
$T_{\mathrm{eff}}$. $\alpha_T$ depends on the zenith angle and the
opacity and we follow \citet{tramontini2019middle} to calculate its
value. To account for the muon rate variations caused by temperature
changes in the atmosphere we use the following equation:
{\vspace*{.25pc}}
{\begin{equation} \label{eq:filtered_muon_rate}
\left(\frac{\rmDelta R}{\langle R \rangle }\right)_{\mathrm{corrected}}
= \frac{\rmDelta R}{\langle R \rangle} - \alpha_T \times \frac{\rmDelta
T_{\mathrm{eff}}}{\langle T_{\mathrm{eff}}\rangle},
\end{equation}}\unskip
\par\vspace*{.35pc}\noindent
where $\rmDelta R = R - \langle R \rangle$, $\rmDelta T_{\mathrm{eff}} =
T_{\mathrm{eff}} - \langle T_{\mathrm{eff}} \rangle$, $\langle R
\rangle$ is the mean muon rate and $\langle T_{\mathrm{eff}} \rangle$
is the mean effective temperature. We compute $T_{\mathrm{eff}}$ using
the same 30-day width Hamming moving average window as $R$, and we
compute $\langle T_{\mathrm{eff}} \rangle$ and $\langle R \rangle$ for
each uninterrupted acquisition period separately.

\subsection{Time series of opacity} \label{sec3.3}

From the measured integrated muon flux, and the calibration experiment
that provides the experimental acceptance, we estimate the
corresponding opacity $\varrho$ (mwe), using the same tabulated values
of the $(I,\varrho,\theta)$ dependence that we use for computing the
muon image. As with $R$, we merge the signals from adjacent directions
together when studying time variations to improve the temporal
resolution of the time series. We compute the opacity variations time
series as $\rmDelta\varrho=\varrho-\langle\varrho \rangle$, where
$\langle\varrho\rangle$ is the mean opacity. We compute $\rmDelta\varrho$
using the same 30-day width Hamming moving average window as $R$, and
we compute $\langle \varrho\rangle$ for each uninterrupted acquisition
period \mbox{separately}.

\subsection{Time series of groundwater content} \label{sec3.4}

In this work, we use two parameters from the ERA5 and ERA5-Land climate
reanalysis datasets, which are produced by the European Centre for
Medium-Range Weather Forecasts (ECMWF). These datasets assimilate
global observations into atmospheric models to create a comprehensive
and consistent dataset using the laws of physics \citep{era5}. ERA5
provides hourly data on atmospheric, land-surface, and sea-state
parameters. The data are available on regular latitude-longitude grids
with a resolution of 0.25{\textdegree} ${\times}$ 0.25{\textdegree}
(${\sim}$31~km). The data covers atmospheric parameters at 37 pressure
levels ranging from 1000~hPa to 1~hPa. For various land applications,
ERA5-Land offers an improved accuracy by replaying the land component
of ERA5 climate reanalysis \citep{era5land}. ERA5-Land also provides an
enhanced resolution of 0.1{\textdegree} ${\times}$ 0.1{\textdegree}
(${\sim}$11~km) and offers parameters on the four levels of the ECMWF
surface model. We use these data to compute the time series of
effective temperature in the atmosphere and correct for this effect in
our muon time series, as presented previously.

From ERA5, we use the atmospheric temperature, spanning all 37
atmospheric pressure levels. From ERA5-LAND, we take the volumetric
soil water content estimated across the four layers of the ECMWF
surface model, $v_i$ ($\text{m}^3{\cdot}\text{m}^{-3}$) ($i=1$, 2, 3,
4), which together encompasses the first 289~cm of soil. The
corresponding thickness of each layer are $\rmDelta h_1=7$~cm, $\rmDelta
h_2=21$~cm, $\rmDelta h_3= 72$~cm and $\rmDelta h_4=189$~cm for the layer
1, 2, 3 and 4, respectively. This parameter is associated with soil
texture or classification, soil depth, and the underlying groundwater
level. Both parameters are listed four times a day (00:00, 06:00,
12:00, and 18:00~hr).

In addition, we analyse the stream flow of the Doubs River, measured at
the Ocourt station, and the rainfall registered at the Montenol station
for estimating the amount of groundwater within the study region
(Figure~\ref{fig:intro:geology_scenario}a). The distance between these
two stations is approximately 5~km. These datasets play an important
role in hydrological processes and can serve as indicators of
groundwater storage changes. Stream flow and rainfall data provide
information about groundwater discharge and recharge, respectively. In
order to compare the observed muon rate fluctuations to variations in
groundwater content, we define the weighted volumetric soil water in
the first 289~cm of soil, $v~(\text{m}^3{\cdot}\text{m}^{-3})$, as $v =
({\sum_{i=1}^{4}v_i \times \rmDelta h_i})/({\sum_{i=1}^{4} \rmDelta h_i})$.
We compute $v$ using a \mbox{30-day} width Hamming moving average
window. We define $\rmDelta v = v - \langle v \rangle$, where $\langle v
\rangle$ is the mean weighted volumetric soil and is computed for each
uninterrupted acquisition period\break separately.

\section{Results} \label{sec4}

From the calibration experiment, we measured a total acceptance of
$(1385 \pm 67)~\mathrm{cm}^2{\cdot}\mathrm{sr}$, which we then used to
calculate the muon image of average density for the whole measured
period. In Figure~\ref{fig:results:radiograph_geology} we present it,
projected over the surface given by the intersection of the topography
and each axis of observation, and superposed to the geological map
(Figure~\ref{fig:intro:telescope}a). We observe two distinct
regions at each side of the valley: a high-density region to the north
of the valley (HD) and a low-density region to its south (LD).
Additionally, we observe a smaller region characterized by high
densities located behind the southern ridge. Table~\ref{table:summary}
summarizes the volume of rock, average rock thickness, average muon
rate, average opacity, average density and the estimated acceptance for
the HD, LD, and combined (HD and LD) regions. Even if the LD region is
characterized by a lower density, it has a higher opacity because of
the larger amount of rock traversed by\break the muons.

\begin{figure*}
\includegraphics{fig03}
\caption{\label{fig:results:radiograph_geology}Muon image of average
density between the muon detector, located inside the Mont-Terri
tunnel, and the surface. The detector's axes of observation are
represented with gray lines. The average density values are plotted
over the topography, together with the geological map. The black, green
and blue dashed lines indicate the motorway tunnel, the HD region, and
the LD region, respectively. A nearest-neighbour interpolation is used
to account for the angular resolution of the axes of observation. The
geological map is based on the Geological Atlas of Switzerland
1:25,000, map No.~40 ``St-Ursanne'', and is available through
\url{https://map.geo.admin.ch}.}
{\vspace*{2pt}}
\end{figure*}

\begin{table*}%tab1
\caption{\label{table:summary}The scanned volume of rock $V$, the average
opacity $\langle\varrho\rangle$, the average muon rate $\langle
R\rangle$, the average rock thickness $\langle L\rangle$, the average
density $\langle\rho\rangle$ and the estimated acceptance $\mathcal{T}$
in each region}
\begin{tabular}{ccccccc}
\thead
Region & $V$ (m$^3$) & $\langle\varrho\rangle$ (mwe) & $\langle R\rangle$
($\text{d}^{-1}$) & $\langle L\rangle$ (m) & $\langle\rho\rangle$
($\text{g}{\cdot}\text{cm}^{-3}$) & $\mathcal{T}$
($\text{cm}^2{\cdot}\text{sr}$)\\
\endthead
High-density (HD) & \0$4.0 \times 10^6$ & $600 \pm 50$\0 & $605 \pm 7$\0 & $230 \pm 14$ & $2.82 \pm 0.2$ & \0$984 \pm 45$\\
Low-density (LD) & \0$7.7 \times 10^6$ & $800 \pm 200$ & $194 \pm 5$\0 & $339 \pm 81$ & \0$2.5 \pm 0.1$ & \0$401 \pm 22$\\
HD and LD & $11.0 \times 10^6$ & $670 \pm 160$ & $800 \pm 10$ & $270 \pm 72$ & \0$2.7 \pm 0.2$ & $1385 \pm 67$
\botline
\end{tabular}
\tabnote{For $\langle\varrho\rangle$, $\langle R\rangle$,
$\langle L\rangle$ and $\langle\rho\rangle$, the dispersion across the
axes of observation within each region is also provided.}
\end{table*}

Similar patterns can be observed between the average density values and
the geological formations (cf.\ 
Figure~\ref{fig:results:radiograph_geology}). The HD region corresponds
to sandy limestone formations that constitute the local aquifer (the
Ifenthal, Hauptrogenstein and Passwang formations). The LD region
corresponds to the local aquifer as well as the fissured and karstic
limestones that compose the regional aquifer (mainly the Saint-Ursanne
formation in the study region), along with the aquitard (the Barschwill
formation). Mean density is spatially variable. Mid-density values are
found on the crest of the mountain, which might also be linked to the
larger thickness of the weathered rock. The high-density values located
behind the southern ridge are partly correlated to the presence of the
Saint-Ursanne formation. Note that as the density values obtained are
averages through the muon paths, these values not only reflect the
density of the outcrop formations, but also of the geological
formations traversed in the subsurface. For instance, the LD region
also intersects the Ifenthal, Hauptrogenstein and Passwang formations
(Figure~\ref{fig:intro:telescope}a).

We compute the muon rate time series for the three selected regions
(cf.\ Methodology). As an example, the $\langle R \rangle$ and $\langle
T_{\mathrm{eff}} \rangle$ time series of the LD region are shown in
Figure~\ref{fig:results:monitoring_absolute_rate}. We can observe a
clear seasonal variation in the two time series, with lower muon rates
and effective temperatures during winter. To correct the time series
for the atmospheric temperature effect, we computed a specific
effective temperature coefficient $\alpha_{T}$ (cf.\ 
Table~\ref{table:alpha}). $\alpha_T$ remains constant within the
defined uncertainty range. This arises from the dependence of
$\alpha_T$ on the mean value of $E_{\mathrm{min}} \times \cos(\theta)$,
which is associated with the production and propagation of muon parent
particles through the atmosphere \citep{grashorn2010atmospheric,
gaisser2016cosmic}. While $E_{\mathrm{min}}$ is higher in the LD
region, the corresponding $\cos(\theta)$ is lower. Conversely, in the
HD region, $E_{\mathrm{min}}$ is lower, but $\cos(\theta)$ is higher.
As a result, these opposing effects on $E_{\mathrm{min}}\times
\cos(\theta)$ balance each other. The corrected muon rate time series
for each region are shown in Figure~\ref{fig:results:monitoring}, now
in terms of their relative changes $({\rmDelta R}/{\langle R
\rangle})_{\mathrm{corrected}}$ computed using
Equation~(\ref{eq:filtered_muon_rate}). The time series of the HD and
(HD and LD) regions seem stable after the temperature correction. But
the LD region presents significant muon rate variations throughout the
acquisition period, with amplitude variations from ${-}$4\% to 4\%
relative to the time series average value of 194 muons per day, and
well above the 95\% confidence interval.

\begin{figure*}
\includegraphics{fig04}
\caption{\label{fig:results:monitoring_absolute_rate}Muon rate $R$ in
the Low-Density (LD) region and effective temperature
$T_{\mathrm{eff}}$ time series represented as blue and orange lines,
respectively. The colored surface delimits the 95\% confidence
interval. To illustrate the effect of the correction derived from
Equation~(\ref{eq:filtered_muon_rate}), the corrected muon rate in the
LD region is shown as a gray line, with the corresponding 95\%
confidence interval outlined by dashed lines. The gray surface indicate
a period where the muon data acquisition was interrupted for work at
the Mont Terri URL.}
{\vspace*{.35pc}}
\end{figure*}

\begin{figure*}
\includegraphics{fig05}
\caption{\label{fig:results:monitoring}(a)~Rock volume and axes of
observation involved in each region studied. (b)~Deviation from the
mean of the filtered muon rate $({\rmDelta R}/{\langle R
\rangle})_{\mathrm{corrected}}$ time series. The colored surface
delimits the 95\% confidence interval. The gray surface indicate a
period where the muon data acquisition was interrupted for work at the
Mont Terri URL.}
\end{figure*}

\begin{table}%tab2
\caption{\label{table:alpha}Computed values of $\alpha_{T}$ for each
\mbox{region}\vspace*{-1pt}}
\begin{tabular}{cc}
\thead
Region & $\alpha_{T}$\\
\endthead
HD & $0.65 \pm 0.02_{\mathrm{stat}} \pm 0.03_{\mathrm{syst}}$\\
LD & $0.63 \pm 0.02_{\mathrm{stat}} \pm 0.03_{\mathrm{syst}}$\\
HD and LD & $0.65 \pm 0.02_{\mathrm{stat}} \pm 0.03_{\mathrm{syst}}$
\botline
\end{tabular}
\end{table}

To analyze the muon rate changes in the LD region, we present estimated
time series of the relative soil water content, the Doubs river stream
flow, and the rainfall precipitation
(Figure~\ref{fig:monitoring_hydro}, cf.\ Methodology). Note that the
axis of the muon rate variations has been inverted in order to make the
comparison with the relative soil water content straightforward. We
also present the absolute opacity variations estimated from the
corrected muon rate variations, together with the absolute volumetric
soil water content (see Methodology). While the river streamflow and
rainfall variations are not simple to analyze, the soil water content
variations show similarities between groundwater content and the muon
rate, and absolute opacity changes. Higher values of opacity (or,
conversely, lower muon rates) seem to follow the increase in soil water
content. However, we notice that in terms of groundwater storage, the
scales are very different. If we consider for example the first
increase in volumetric soil water content
(Figure~\ref{fig:monitoring_hydro}b), a variation of approximately
$0.06~\mathrm{m}^3{\cdot}\mathrm{m}^{-3}$ can be observed. Since the
total length of the considered layer is of 289~cm, this soil water
content variation corresponds to a change of 0.17~mwe of opacity.
However, the corresponding absolute opacity change is of 12~mwe, that
is, almost an order of magnitude larger. We discuss the possible
explanations for this misfit in the following Section.

\begin{figure*}
\includegraphics{fig06}
\caption{\label{fig:monitoring_hydro}(a)~Deviation from the mean of the
filtered muon rate $({\rmDelta R}/{\langle R
\rangle})_{\mathrm{corrected}}$ in the LD region and deviation from the
mean of the weighted volumetric soil water ${\rmDelta v}/{\langle v
\rangle}$ represented as blue and purple lines, respectively.
(b)~Opacity deviations from the mean value in the LD region, and
weighted volumetric soil water represented as green and purple lines,
respectively. (c)~Doubs river stream flow and rainfall datasets are
represented as purple and gray bars, respectively. The colored surfaces
delimits the 95\% confidence interval. The gray surface indicate a
period where the muon data acquisition was interrupted for work at the
Mont Terri URL.}
\end{figure*}

\section{Discussion} \label{sec5}

Our analyses suggest that the muon method is sensitive to the different
rock formations above the Mont Terri URL
(Figure~\ref{fig:results:radiograph_geology}), and that the corrected
muon rate variations are linked to groundwater variations in these
formations. But these variations are not observed in all the regions
scanned by the detector, and it is not clear whether the difference in
the variations observed are due to the presence of aquifers and
aquitards, or to a systematic sensitivity effect due to the experiment
design. To analyse the relation between the muon rate and
\mbox{groundwater} variations, we first bring into consideration
possible \mbox{external} sources of muon rate variations other than the
effective atmospheric temperature effect. We then propose a method to
quantify the influence of the water table variation in our particular
experiment, and then propose an interpretation for the{\break} results
obtained.

\subsection{External sources of muon rate variations} \label{sec5.1}

The natural flux of muons arriving to the Earth's surface can vary
according to other effects than the atmospheric temperature change
\citep{spurio2014particles}. We need to consider these effects before
concluding that the registered variations are due to density changes in
the scanned rock volume.

Due to the rather high opacity conditions under which we conducted this
experiment (Table~\ref{table:summary}), the impact of barometric
variations on the muon rate is negligible. Only low-energy muons of
around 10~GeV are influenced by atmospheric pressure fluctuations
\citep{lesparre2010geophysical}, and these are filtered out within the
first meters of rock in our experimental setup. In our configuration,
the minimum energy of detected muons ranges from ${\sim}$140~GeV to
${\sim}$410~GeV, when considering a medium of standard rock with a
density of $2.65~\mathrm{g}/\mathrm{cm}^3$. This correspond to a 
\mbox{barometric} coefficient $\beta_{p}$ below $10^{-2}$\%~hPa
\citep{sagisaka1986atmospheric}. As a result, the effect on ${\rmDelta
R}/{\langle R \rangle}$ is of the order of 0.15\%, calculated as the
product of $\beta_{p}$ and the difference between atmospheric pressure
and its mean value.

The Earth's magnetic field deflects low-energy cosmic rays, stopping
them from entering the atmosphere through the magnetosphere and
generating muons. This geomagnetic influence is significant for muons
arriving at sea level with energies up to ${\sim}$5~GeV
\citep{spurio2014particles}. Considering that muons with energies
between 1~GeV and 100~GeV loose energy as they traverse matter at a
rate of ${\sim}$0.25~GeV/mwe \citep{groom2001muon}, the muons affected by
the geomagnetic field variations would only be able to traverse
${\sim}$20~mwe. Similarly to the case of pressure variations, because of
the overburden opacity of the Mont Terri URL the geomagnetic effect is
negligible. In addition, temporal variations originating from
fluctuations in primary cosmic rays or in the geomagnetic field due to
solar wind typically manifest on significantly shorter (e.g., seconds
to hours) or longer (e.g., an ${\sim}$11-year solar cycle) time scales.

\subsection{Muon rate changes due to groundwater{\hfill\break}
variations} \label{sec5.2}

The temporal variations shown in Figure~\ref{fig:results:monitoring}
reveal a fundamental difference between the LD and HD regions. This
suggests a distinct cause affecting the muon rate in each region, as
opposed to a global effect such as the formerly mentioned, or a
malfunctioning of the muon detector. The differences in the muon rate
behavior could be related to the particular geological formations
present in each region, and the difference in inflow rates and
groundwater storage associated to these formations. According to
\citet{marschall2004mont}, the LD region has inflow rates higher than
1~l/s, and in the HD region the inflow rates are in the range between
12~l/s and 200~l/s, which represents a difference up to two orders of
magnitude. A lower inflow rate can result in groundwater accumulating
and discharging over longer time periods, within the LD region,
rendering the density changes due to groundwater variations visible
through muon observations. This stands in contrast to the HD region,
where the higher inflow rates would allow the infiltrated water to flow
rapidly without generating a density contrast for a significant amount
of time. However, it should be noted that the inflow rate measurements
were done locally within the Mont Terri URL. In karstic formations,
inflow rates can exhibit large spatial variability due to their
inherent heterogeneity \citep{goldscheider2007methods}. Therefore, the
difference in muon rate variations may not be solely attributed to the
difference in groundwater variations in the two regions, but could also
be influenced by local-scale karstic heterogeneity.

Unfortunately, there are no instrumented boreholes in the region to
compare the muon rate changes in the LD zone with water table
variations. This is why we resorted to 3 parameters of hydrological
interest depicted in Figure~\ref{fig:monitoring_hydro}: precipitation,
stream flow of a nearby river (cf.\ 
Figure~\ref{fig:intro:geology_scenario} for locations), and a
model-based estimation of soil water content in the near subsurface
(cf.\ Methodology). The latter is the parameter of most interest because
it provides a volumetric mass change, which is what directly affects
the muons traversing the rock. Its drawback is that it only
contemplates the first 289~cm of soil, and it is a regional estimate
from satellite observations. It is important to note that soil
hydraulic properties can exhibit significant variability at local
scales, ranging from meters to hundreds of meters
\citep{Vereecken2016}. Therefore, satellite-derived models should be
interpreted with caution, particularly when comparing them to
measurements conducted at the spatial resolution of muography analysis.
Nonetheless, we found that this estimate correlates well with the local
Doubs stream flow variations, and therefore we use it as a proxy for
the groundwater variations. We present it as a time series of relative
changes (${\rmDelta v}/{\langle v\rangle}$) in the same way as the muon
rate variations, and of absolute values to compare to the absolute
opacity estimates (Figure~\ref{fig:monitoring_hydro}b).

In Figures~\ref{fig:monitoring_hydro}a and~\ref{fig:monitoring_hydro}b,
we observe similar variation patterns between the muon rate and
inferred opacity, and the groundwater variations. Despite not being
coincident throughout the whole acquisition period, the anticorrelation
between these curves is clearer in the first uninterrupted period of
muon acquisition, and presents the same trend during the second one.
Three distinct peaks in volumetric soil water content (November 2016,
March 2017 and January 2018) correspond to muon rate minima, which
appear with a slight time shift. This behaviour is consistent with our
expectations from muon observations, where an increase in groundwater
content augments the opacity along muon trajectories, conversely
resulting in a reduction in the associated muon rate. The relative soil
water decrease around November 2017, however, does not have a
corresponding muon increase. This discrepancy may be due to local-scale
karstic heterogeneity, where variations in groundwater storage or
preferential flow paths at smaller scales could affect the muon rate
response, regardless of the overall soil water content trends. Besides
the correlation trends that follow the expected behaviour, the absolute
opacity estimates shown in Figure~\ref{fig:monitoring_hydro}b allow us
to estimate the groundwater storage changes directly from the muon
measurements. As stated in the Results section, although the trends are
similar, the amplitude of these variations are very different from the
groundwater storage changes assessed by the satellite observations. 

To understand this discrepancy, we consider a simple theoretical model
to compute the impact of a groundwater content change in the muon
measurements, along the different axes of observation
(\mbox{Figure}~\ref{fig:discussion:waterLayer}). We consider a layer of
water-saturated rock of height $\rmDelta z$~(m), located at a depth
$d$~(m) following the topography. The saturated layer height is given
by $\rmDelta z = w/p$, where $w$~(m) and $p$ correspond to the amount of
infiltrated water, and the layer porosity, respectively. We then
compute $\rmDelta L$~(m); the length of the muon trajectory $L$ that lies
inside the considered layer. $\rmDelta L$ thus corresponds to the
distance traveled by muons, in a certain axis of observation, inside
the water-saturated layer. Figure~\ref{fig:discussion:waterLayer}
illustrates the model with an example, where an axis of observation
intersects a water-saturated layer of height $\rmDelta z = 2$~m in two
segments, resulting in muons traversing a total distance of $\rmDelta L
=\rmDelta L_1 + \rmDelta L_2 \simeq 30$~m within the water-saturated layer,
where $\rmDelta L_1$ and $\rmDelta L_2$ represent the distance travelled by
muons in each segment. We define the Amplification Factor $\mathcal{A}$
($\text{mwe}{\cdot}\text{m}^{-1}$) as the ratio between the opacity
change, due to groundwater change along the axis of observation, to the
actual amount of infiltrated water, which is the same for every
viewpoint: $\mathcal{A} = \rmDelta \varrho/w$. $\rmDelta \varrho$ is given
by $\rmDelta \varrho = p \times \rmDelta L \times \rho_{\mathrm{water}}$,
where $\rho_{\mathrm{water}} =1~\mathrm{g}{\cdot}\mathrm{cm}^{-3}$ is
the density of water.

\begin{figure*}
\includegraphics{fig07}
\caption{\label{fig:discussion:waterLayer}Illustration of the
water-saturated layer model used to calculate the impact of groundwater
variations in muon measurements. A single axis of observation of the
muon detector is represented as an orange straight line, the topography
is represented as a green curve, and the water-saturated layer of
height $\rmDelta z$ at a depth $d$ is represented as a blue surface.
$\rmDelta L=\rmDelta L_1+\rmDelta L_2$ corresponds to the distance traveled
by muons in the water-saturated layer.}
{\vspace*{.85pc}}
\end{figure*}

\looseness=-1
In Figure~\ref{fig:discussion:ampFactor} we present the computed
Amplification Factor $\mathcal{A}=\mathcal{A}(w,d,p)$ in each axis of
observation for $w=0.3$~m, $d=10$~m and $p=0.3$. These values were
arbitrarily chosen in order to obtain an initial approximation of the
amplification effect. To better account for the topography of the site,
we synthetically increased the number of axes of observation from $31
\times 31$ to $61 \times 61$ to have a more precise mapping, but the
surface of investigation remains the same as in our real muon
experiment. We find that while $\mathcal{A} \sim
1~\mathrm{mwe}{\cdot}\mathrm{m}^{-1}$ in the HD region, indicating a
negligible amplification effect, $\mathcal{A} >
1~\mathrm{mwe}{\cdot}\mathrm{m}^{-1}$ in the LD region, that is, in the
region where significant muon rate variations are observed
(Figure~\ref{fig:discussion:ampFactor}). We assess the influence of
each parameter on the computation of $\mathcal{A}$ by modifying the
values of porosity ($N_{\mathrm{p}}=4$ linearly spaced values from 0.2 to 0.5),
infiltrated water ($N_{\mathrm{w}}=10$ linearly spaced values from 0.1~m to 3~m),
and the layer depth ($N_{\mathrm{d}}=33$ linearly spaced values from 6~m to
102~m). From these $N_{\mathrm{p}} \times N_{\mathrm{w}}\times N_{\mathrm{d}} = 1320$ computations, we
determine that $\mathcal{A}$ primarily depends on the depth at which
water starts to accumulate, $d$. This is shown in
Figure~\ref{fig:discussion:ampFactor_3D}a, where we present the mean
value of $\mathcal{A}$, $\langle\mathcal{A}\rangle$, of the study
region by considering all the axes of observation for each parameter
combination. Figure~\ref{fig:discussion:ampFactor_3D}b shows the
corresponding standard deviation of $\mathcal{A}$,
$\sigma_{\mathcal{A}}$.

\begin{figure*}
\includegraphics{fig08}
{\vspace*{.2pc}}
\caption{\label{fig:discussion:ampFactor}Amplification factor
expressing the ratio of opacity change in each axis of observation to
an increase in groundwater storage of 0.3~m (cf.\ 
Figure~\ref{fig:discussion:waterLayer}). The amplification factor
values are plotted over the topography. The green and blue dashed lines
indicate the HD and LD regions, respectively. A nearest-neighbour
interpolation is used to account for the angular resolution of the axes
of observation.}
\end{figure*}

\begin{figure*}
{\vspace*{-.3pc}}
\includegraphics{fig09}
{\vspace*{-.3pc}}
\caption{\label{fig:discussion:ampFactor_3D}(a)~Mean value of
$\mathcal{A}$, $\langle\mathcal{A}\rangle$, for the combined HD and LD
regions in function of depth of the water-saturated layer $d$ and the
infiltrated water $w$, and for different porosity values $p$.
(b)~Standard deviation of $\mathcal{A}$, $\sigma_{\mathcal{A}}$, for
the combined HD and LD regions in function of $d$ and $w$, and for
different porosity values~$p$.}
\end{figure*}

Based on this result, we compute the mean Amplification Factor as a
function of $d$, $\overline{\mathcal{A}} = \overline{\mathcal{A}}(d)$,
for the LD region, the HD region, and the combined HD and LD regions.
$\overline{\mathcal{A}}$ is computed as the weighted mean value of
$\langle\mathcal{A}\rangle$ for each region separately, when modifying
the amount of infiltrated water and the porosity of the rock as
described. The weighted mean is computed using inverse-variance weights
and is defined by:
{\vspace*{.25pc}}
{\begin{equation} \label{eq:ampFactor_vs_d_value}
\overline{\mathcal{A}}(d)=\frac{\sum_{i=1}^{N_{\mathrm{w}}}\sum_{j=1}^{N_{\mathrm{p}}}
\frac{\langle\mathcal{A}\rangle(w_i,p_j,d)}{\sigma_{\mathcal{A}}
(w_i,p_j,d)}}{\sum_{i=1}^{N_{\mathrm{w}}}\sum_{j=1}^{N_{\mathrm{p}}}\frac{1}
{\sigma_{\mathcal{A}}(w_i,p_j,d)}}.
\end{equation}}\unskip
The corresponding standard error is given by:
{\vspace*{.25pc}}
{\begin{equation} \label{eq:ampFactor_vs_d_error}
\sigma_{\overline{\mathcal{A}}}(d)=\sqrt{\frac{1}{\sum_{i=1}^{N_{\mathrm{w}}}
\sum_{j=1}^{N_{\mathrm{p}}}\frac{1}{\sigma_{\mathcal{A}}(w_i,p_j,d)}}}.
\end{equation}}\unskip
The results are shown in
Figure~\ref{fig:discussion:ampFactor_regions_mean_std}. We confirm that
the amplification effect is large in the LD region, with average values
larger than twice the amount of water infiltrated when the layer is
above 40~m depth. In the HD region, there is a negligible effect of the
geometry, since the axis of observation are almost perpendicular to the
saturated layer. Last, when considering the LD and HD regions together,
there is a mild average amplification, as expected from the relative
contribution of the LD region. 

\begin{figure*}
\includegraphics{fig10}
{\vspace*{-.2pc}}
\caption{\label{fig:discussion:ampFactor_regions_mean_std}Mean
Amplification Factor $\overline{\mathcal{A}}$ for the LD region (blue),
the HD region (green) and the combined HD and LD regions (black) as a
function of the depth of the water table. Each data point is calculated
considering all the axes of observation belonging to the considered
region. The error bars represent the standard error of the weighted
mean using variance-defined weights, and account for the dispersion due
to different porosities and water infiltrations.}
{\vspace*{-.4pc}}
\end{figure*}

Our geometrical considerations thus suggest that the sensitivity or not
to groundwater variations seem to depend largely on the geometry of
acquisition of the muon experiment. This is an exciting result for the
future applications of muography in hydrology, since it means that the
sensitivity could be increased by considering a good setup, where
ideally the detectors axes of observation should be oriented almost
parallel to the varying groundwater table. In our case, this model
allows to explain why we only observe muon rate changes in the LD
region. The model also suggests that the amplitude discrepancies
between the satellite and muon estimates of groundwater variation may
arise from the geometrical amplification factor. However, we note that
the misfit between the soil water content variations and the absolute
observed opacity values is not fully explained by the amplification
factor alone. Several other factors likely contribute to this
discrepancy. One is the uncertainty in the absolute observed opacity
(Figure~\ref{fig:monitoring_hydro}b). Another is that
the amplification factor only accounts for the
effect of a water-saturated layer parallel to the surface, while water
content variations may also occur in more complex geometries. Moreover,
local-scale{\break} features, such as preferential flow paths and karst
heterogeneity, may influence water distribution in ways not captured by
regional-scale satellite data.

\section{Conclusions} \label{sec6}

We used, for the first time, a continuous acquisition of Muography data
to estimate groundwater storage variations in the subsurface. After
correcting for the only significant external effect on the variation of
the incoming muon flux, the temperature changes in the middle
atmosphere, we found two distinct regions (LD and HD) where the muon
rate evolution over time is different: in the LD region we observe muon
rate changes that are correlated to estimates of volumetric soil water
content from satellite{\break} measurements, whereas the HD region does
not present significant variations over the period of measurements. We
attribute this difference to (1)~the difference in the geometry in
which the muon trajectory intersects the saturation-varying layer,
illustrated by the Amplification Factor in the Discussion section, and
(2)~to a less extent, the LD region's lower inflow rates potentially
facilitating more observable groundwater \mbox{phenomena}. 

The sensitivity analyses on the Amplification Factor reveal the large
impact of the geometry given by the trajectories of detected muons and
the \mbox{topography}. We thus point to the possibility of enhancing
the sensitivity of the method in groundwater storage estimations by
adequately setting the experiment geometry. From the study of the Mont
Terri URL, we find that the amplification factor primarily depends on
the depth of water accumulation, rather than the rock porosity or the
amount of infiltrated water. This underscores the added value of
muography in hydrogeophysics, potentially shedding light not only on
the groundwater storage variations, but also on assessing the depth of
the water table. \looseness=-1

{\vspace*{-2pt}}

\section*{Declaration of interests}

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

{\vspace*{-2pt}}

\section*{Funding}

This study of the MD experiment of the Mont Terri project
(\url{https://www.mont-terri.ch}) is funded by Swisstopo.

{\vspace*{-2pt}}

\section*{Acknowledgments}

We extend our
gratitude to Thierry Theurillat and Senecio Schefer for their technical
and logistical assistance at Mont Terri URL, and to Pauline Baland for
providing the geological map, as well as David Jaeggi and Vincent
Regard for providing the rainfall dataset.\ The ECMWF data is available
from \url{https://www.ecmwf.int/}. The Doubs river stream flow is
available from \url{https://www.hydrodaten.admin.ch/}.\ Muon data used
for all calculations are displayed in figures and are available in the
supplementary material. This study uses perceptually uniform color maps
to avoid visual distortion of the data and exclusion of readers with
colour-vision deficiencies \citep{crameri2020colormaps}.

\back{}

{\vspace*{-2pt}}

\section*{Supplementary data}
Supporting information for this article is available on the journal's
website under \printDOI\ or from the author.  

%Supplementary material for this article is supplied as a separate
%archive Table\_S1.xlsx, the \mbox{related} data is described in the document
%supplementary\_material.pdf.

\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{crgeos-277-suppl.pdf}}
\CDRsupplementaryTwotypes{supplementary-material}{\cdrattach{Table-S1.xlsx}}

\bibliographystyle{crgeos}
\bibliography{crgeos20240525}
\refinput{crgeos20240525-reference.tex}

\end{document}
