# Defining the sensitivity of cosmic ray muons to groundwater storage changes

Matías Tramontini, Marina Rosas-Carbajal, Fabio Iván Zyserman, Laurent

Longuevergne, Christophe Nussbaum, Jacques Marteau

► To cite this version:

Matías Tramontini, Marina Rosas-Carbajal, Fabio Iván Zyserman, Laurent Longuevergne, Christophe Nussbaum, et al. Defining the sensitivity of cosmic ray muons to groundwater storage changes. Comptes Rendus. Géoscience, 2024, 356 (G1), pp.177-194. 10.5802/crgeos.277. insu-04797563

HAL Id: insu-04797563

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

Submitted on 22 Nov 2024

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

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Matías Tramontini, Marina Rosas-Carbajal, Fabio Iván Zyserman, Laurent Longuevergne, Christophe Nussbaum and Jacques Marteau

Defining the sensitivity of cosmic ray muons to groundwater storage changes

Volume 356 (2024), p. 177-194

Online since: 20 November 2024

https://doi.org/10.5802/crgeos.277

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

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Research article

# Defining the sensitivity of cosmic ray muons to groundwater storage changes

Matías Tramontini<sup>a, *, a</sup>, Marina Rosas-Carbajal<sup>b</sup>, Fabio Iván Zyserman<sup>a</sup>, Laurent Longuevergne<sup>b, c</sup>, Christophe Nussbaum<sup>b, d</sup> and Jacques Marteau<sup>c, e</sup>

<sup>a</sup> CONICET - Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata, Argentina

<sup>b</sup> Université de Paris, Institut de Physique du Globe de Paris, CNRS, UMR 7154, F-75238 Paris, France

<sup>c</sup> Geosciences Rennes - UMR 6118, Univ. Rennes, CNRS, 35000 Rennes, France

<sup>d</sup> Swiss Geological Survey at Swisstopo, Seftigenstrasse 264, CH-3084 Wabern, Switzerland

<sup>e</sup> Institut de Physique des 2 Infinis de Lyon, UMR 5822, CNRS-IN2P3, Université de Lyon, Université Claude Bernard Lyon 1, France

E-mail: mtramontini@fcaglp.unlp.edu.ar (M. Tramontini)

**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.

**Keywords.** Muography, Groundwater variations, Hydrogeophysics.

**Funding.** Swisstopo.

Manuscript received 28 June 2024, revised 21 October 2024, accepted 29 October 2024.

## 1. 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

<footnote>* Corresponding author</footnote>
"availability of life-sustaining resources" [Gaillardet et al., 2018]. 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 [Gleeson et al., 2014]. 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 [Alley et al., 2002].

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 [Hermans et al., 2023]. These are indirect methods that infer different physical properties of the subsurface, such as electrical resistivity [McGarr et al., 2021], seismic velocity [Blazevic et al., 2020], density [Crossley et al., 2013], spontaneous potential [Grobbe et al., 2021], and hydrogen content of pore fluid [Schmidt and Rempe, 2020], 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² to 5000 km² [Castellazzi et al., 2016], 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 [Schmidt et al., 2008], 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 [Longuevergne et al., 2010]. These techniques are effective for monitoring groundwater fluctuations at regional and basin scales [Leblanc et al., 2009, Syed et al., 2008], though their spatial resolution is limited and not suitable for detailed local studies, typically on the order of hundreds of meters.

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 [Evans et al., 2016, Köhli et al., 2015]. Similarly, the attenuation of underground cosmic-ray electron flux has been proven to be sensitive to variations in soil water content following rainfall events [Taketa et al., 2022]. 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 [Bonechi et al., 2020]. 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 [Bonomi et al., 2020], including imaging of the density distribution of volcanic edifices [Bajou et al., 2023, Rosas-Carbajal et al., 2017], detecting hidden cavities [Mishima et al., 2017, Cimmino et al., 2019], reconstructing the bedrock morphology underneath an active glacier [Nishiyama et al., 2019], mapping underground structures [Borselli et al., 2022], mining exploration of ore bodies [Beni et al., 2023], and detecting the presence of discontinuities produced by karstic networks [Lesparre et al., 2016], among other applications [Cohu et al., 2023].
The continuous acquisition of muon data has shown promising potential to monitor hydrothermal activity in volcanoes, such as La Soufrière de Guadeloupe in the Lesser Antilles, France [Jourde et al., 2016a, Le Gonidec et al., 2019]. Muography has also been utilized for a time-lapse study of mass density changes over time [Oláh et al., 2023], and for monitoring gravity-driven debris flows in the Sakurajima volcano, Japan [Oláh et al., 2021]. 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 [Osprey et al., 2009]. 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 [Tramontini et al., 2019], since this effect becomes increasingly significant as the overburden opacity of the site increases [Grashorn et al., 2010].

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.

## 2. Experimental setup

The Mont Terri Underground Rock Laboratory (URL) is a research facility located in the Jura Mountains in Switzerland (Figure 1). 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₂ storage, including more recently heat storage [Bossart et al., 2018]. The geological profile and hydrogeological classification at the Mont Terri URL is presented in Figure 2a. 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 Marschall et al. [2004], the local aquifer is primarily composed of sandy limestones that exhibit fissures and karstic features, with inflow rates ranging from 12 l·s⁻¹ to 200 l·s⁻¹ due to localised water-conducting features. Conversely, the regional aquifer comprises fissured and karstic limestones, with inflow rates exceeding 1 l·s⁻¹. 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 hydrogeological studies.

The muon detector used in this study was developed by the DIAPHANE project [e.g., Marteau et al., 2012, 2017]. 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 × 16 pixels of 5 × 5 cm² each (Figure 2b). 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 Marteau et al. [2014]. Due to the large rock volume under study in relation to the size of the detector, we admit a point-like approximation of the de-
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (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 https://www.mont-terri.ch.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** (a) Geological profile adapted from Nussbaum et al. [2017]. Hydrogeological classification from Marschall et al. [2004]. Geological map based on the Geological Atlas of Switzerland 1:25,000, map No. 40 "St-Ursanne", available through https://map.geo.admin.ch. (b) The muon detector deployed in the Mont Terri Underground Rock Laboratory (URL).</figcaption></figure>

detector where events with pixels in the same relative direction are considered to correspond to the same trajectory [Lesparre et al., 2010]. 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)_{\text{front}} - (x, y)_{\text{rear}}$ and $(x, y)_{\text{front}}$ and $(x, y)_{\text{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.

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° and oriented 152°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 [Lesparre et al., 2010].
# 3. Methodology

## 3.1. Computing the muon image

The physical quantity estimated by muography is the opacity, $\rho$ (mwe), which represents the amount of matter traversed by the muons along their trajectories and is given by:

$$ \rho(L) = \int_L \rho(l) d(l) = \langle \rho \rangle \times L, \quad (1) $$

where $\rho$ (g·cm⁻³) is the density of the medium, $l$ is the coordinate measured along the trajectory of longitude $L$ (m), and $\langle \rho \rangle$ (g·cm⁻³) is the average density along the trajectory. $\rho$ is often expressed in mwe (meter water equivalent), where 1 mwe corresponds to 100 g·cm⁻².

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$ (cm⁻²·sr⁻¹·s⁻¹), 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

$$ I(\rho, \theta) = \int_{E_{\min}(\rho)}^{\infty} \phi_{\mu}(E', \theta) dE', \quad (2) $$

where $\theta$ is the zenith angle and $E_{\min}(\rho)$ is the minimum energy required by a muon to traverse a given opacity $\rho$, and $\phi_{\mu}$ (cm⁻²·sr⁻¹·s⁻¹·GeV⁻¹) is the differential muon flux, which represents the flux of muons with a given energy value. $E_{\min}$ is computed by solving numerically $E_{\min} - \int_0^\theta (dE/d\rho) d\rho = E_\mu$, where $E_\mu = 0.10566$ GeV is the muon rest energy and $dE/d\rho$ accounts for the muon energy loss through matter and is described by the attenuation constants provided by the Particle Data Group tables [Workman et al., 2022].

We estimate the integrated muon flux along each axis of observation as $I(\hat{r}_{i,j}) = (N_\mu(\hat{r}_{i,j}))/(\Delta 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, $\Delta T(s)$ is the acquisition duration, and $\mathcal{T}(\hat{r}_{i,j})$ (cm²·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 conditions, where the measured muon rate is compared to

the theoretically expected integrated muon flux $I^{\text{open-sky}}(\hat{r}_{i,j})$:

$$ \mathcal{T}(\hat{r}_{i,j}) = \frac{N_{\mu}^{\text{calib}}(\hat{r}_{i,j})}{\Delta T^{\text{calib}} \times I^{\text{open-sky}}(\hat{r}_{i,j})}, \quad (3) $$

where $N_\mu^{\text{calib}}(\hat{r}_{i,j})$ and $\Delta T^{\text{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 [Bevington and Robinson, 2003]. To calculate $I^{\text{open-sky}}(\hat{r}_{i,j})$ we follow Cohu et al. [2022]: we use the CORSIKA code [Heck et al., 1998] to perform air shower simulations, and incorporate the primary flux parameterization at the top of the atmosphere described by Hörandel [2003].

Once the integrated muon flux is estimated, we use tabulated values of the $(I, \rho, \theta)$ dependence to determine the opacity $\rho(\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, \rho, \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 = \rho(\hat{r}_{i,j})/L(\hat{r}_{i,j})$.

## 3.2. Time series of muon rate

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 [Jourde et al., 2016b]. We compute the average cosmic muon rate time series, $R$, using a 30-day width Hamming moving average window [Hamming, 1998]. The associated confidence interval is determined as the standard deviation of the mean [Bevington and Robinson, 2003]. We assume that the time intervals between consecutive muon detections follow an exponential distribution [Leo, 1994], which is consistent with the assumption that $R$ follows a Poisson process [Gaisser et al., 2016].
Due to the significant amount of rock between the surface and the tunnel, only relatively high-energy muons arrive to our detector. Tramontini et al. [2019] 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_{eff}$ (K) [Grashorn et al., 2010], which is a weighted average of the atmospheric temperature profile, and the effective temperature coefficient $\alpha_T$ [Barrett et al., 1952], which gives the relative change in the muon rate $R$ given a relative change in $T_{eff}$. $\alpha_T$ depends on the zenith angle and the opacity and we follow Tramontini et al. [2019] to calculate its value. To account for the muon rate variations caused by temperature changes in the atmosphere we use the following equation:

$$ \left( \frac{\Delta R}{\langle R \rangle} \right)_{\text{corrected}} = \frac{\Delta R}{\langle R \rangle} - \alpha_T \times \frac{\Delta T_{\text{eff}}}{\langle T_{\text{eff}} \rangle}, \quad (4) $$

where $\Delta R = R - \langle R \rangle$, $\Delta T_{eff} = T_{eff} - \langle T_{eff} \rangle$, $\langle R \rangle$ is the mean muon rate and $\langle T_{eff} \rangle$ is the mean effective temperature. We compute $T_{eff}$ using the same 30-day width Hamming moving average window as $R$, and we compute $\langle T_{eff} \rangle$ and $\langle R \rangle$ for each uninterrupted acquisition period separately.

### 3.3. Time series of opacity

From the measured integrated muon flux, and the calibration experiment that provides the experimental acceptance, we estimate the corresponding opacity $\rho$ (mwe), using the same tabulated values of the $(I, \rho, \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 $\Delta\rho = \rho - \langle\rho\rangle$, where $\langle\rho\rangle$ is the mean opacity. We compute $\Delta\rho$ using the same 30-day width Hamming moving average window as $R$, and we compute $\langle\rho\rangle$ for each uninterrupted acquisition period separately.

### 3.4. Time series of groundwater content

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 [Hersbach et al., 2020]. 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° × 0.25° (~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 [Muñoz-Sabater et al., 2021]. ERA5-Land also provides an enhanced resolution of 0.1° × 0.1° (~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$ (m$^3$·m$^{-3}$) ($i$ = 1, 2, 3, 4), which together encompasses the first 289 cm of soil. The corresponding thickness of each layer are $\Delta h_1 = 7$ cm, $\Delta h_2 = 21$ cm, $\Delta h_3 = 72$ cm and $\Delta 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 1a). 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$ (m$^3$·m$^{-3}$), as $v = (\sum_{i=1}^4 v_i \times \Delta h_i)/(\sum_{i=1}^4 \Delta h_i)$. We compute $v$ using a 30-day width Hamming moving
average window. We define $\Delta\nu = \nu - \langle\nu\rangle$, where $\langle\nu\rangle$ is the mean weighted volumetric soil and is computed for each uninterrupted acquisition period separately.

## 4. Results

From the calibration experiment, we measured a total acceptance of (1385 ± 67) cm²·sr, which we then used to calculate the muon image of average density for the whole measured period. In Figure 3 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 2a). 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 1 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 the muons.

Similar patterns can be observed between the average density values and the geological formations (cf. Figure 3). 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 2a).

We compute the muon rate time series for the three selected regions (cf. Methodology). As an example, the $\langle R \rangle$ and $\langle T_{\text{eff}} \rangle$ time series of the LD region are shown in Figure 4. 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 2). $\alpha_T$ remains constant within the defined uncertainty range. This arises from the dependence of $\alpha_T$ on the mean value of $E_{\text{min}} \times \cos(\theta)$, which is associated with the production and propagation of muon parent particles through the atmosphere [Grashorn et al., 2010, Gaisser et al., 2016]. While $E_{\text{min}}$ is higher in the LD region, the corresponding $\cos(\theta)$ is lower. Conversely, in the HD region, $E_{\text{min}}$ is lower, but $\cos(\theta)$ is higher. As a result, these opposing effects on $E_{\text{min}} \times \cos(\theta)$ balance each other. The corrected muon rate time series for each region are shown in Figure 5, now in terms of their relative changes $(\Delta R/\langle R \rangle)_{\text{corrected}}$ computed using Equation (4). 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.

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 6, 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 6b),
<figure><img src="image_4.png" /></figure>

**Figure 3.** 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 https://map.geo.admin.ch.

<table><caption>**Table 1.** The scanned volume of rock *V*, the average opacity $\langle\rho\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</caption><thead><tr><th>Region</th><th><i>V</i> (m<sup>3</sup>)</th><th>⟨<i>ρ</i>⟩ (mwe)</th><th>⟨<i>R</i>⟩ (d<sup>-1</sup>)</th><th>⟨<i>L</i>⟩ (m)</th><th>⟨<i>ρ</i>⟩ (g·cm<sup>-3</sup>)</th><th><i>T</i> (cm<sup>2</sup>·sr)</th></tr></thead><tbody><tr><td>High-density (HD)</td><td>4.0 × 10<sup>6</sup></td><td>600 ± 50</td><td>605 ± 7</td><td>230 ± 14</td><td>2.82 ± 0.2</td><td>984 ± 45</td></tr><tr><td>Low-density (LD)</td><td>7.7 × 10<sup>6</sup></td><td>800 ± 200</td><td>194 ± 5</td><td>339 ± 81</td><td>2.5 ± 0.1</td><td>401 ± 22</td></tr><tr><td>HD and LD</td><td>11.0 × 10<sup>6</sup></td><td>670 ± 160</td><td>800 ± 10</td><td>270 ± 72</td><td>2.7 ± 0.2</td><td>1385 ± 67</td></tr></tbody></table>

For $\langle\rho\rangle$, $\langle R\rangle$, $\langle L\rangle$ and $\langle\rho\rangle$, the dispersion across the axes of observation within each region is also provided.

<table><caption>**Table 2.** Computed values of $\alpha_T$ for each region</caption><thead><tr><th>Region</th><th>α<sub>T</sub></th></tr></thead><tbody><tr><td>HD</td><td>0.65 ± 0.02<sub>stat</sub> ± 0.03<sub>syst</sub></td></tr><tr><td>LD</td><td>0.63 ± 0.02<sub>stat</sub> ± 0.03<sub>syst</sub></td></tr><tr><td>HD and LD</td><td>0.65 ± 0.02<sub>stat</sub> ± 0.03<sub>syst</sub></td></tr></tbody></table>

a variation of approximately 0.06 m<sup>3</sup>·m<sup>-3</sup> 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.

## 5. Discussion

Our analyses suggest that the muon method is sensitive to the different rock formations above the Mont Terri URL (Figure 3), 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
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Muon rate *R* in the Low-Density (LD) region and effective temperature *T*<sub>eff</sub> 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 (4), 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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Rock volume and axes of observation involved in each region studied. (b) Deviation from the mean of the filtered muon rate ($\Delta R/\langle R \rangle_{\text{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.</figcaption></figure>

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 groundwater

variations, we first bring into consideration possible external sources of muon rate variations other than the effective atmospheric temperature effect. We then propose a method to quantify the influence
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Deviation from the mean of the filtered muon rate ($\Delta R/\langle R\rangle$)corrected in the LD region and deviation from the mean of the weighted volumetric soil water $\Delta 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.</figcaption></figure>

of the water table variation in our particular experiment, and then propose an interpretation for the results obtained.

## 5.1. *External sources of muon rate variations*

The natural flux of muons arriving to the Earth's surface can vary according to other effects than the atmospheric temperature change [Spurio, 2014]. 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 1), 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 [Lesparre et al., 2010], 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 ~140 GeV to ~410 GeV, when considering a medium of standard rock with a density of 2.65 g/cm³. This correspond

to a barometric coefficient $\beta_p$ below $10^{-2}\%$ hPa [Sagisaka, 1986]. As a result, the effect on $\Delta 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 ~5 GeV [Spurio, 2014]. Considering that muons with energies between 1 GeV and 100 GeV lose energy as they traverse matter at a rate of ~0.25 GeV/mwe [Groom et al., 2001], the muons affected by the geomagnetic field variations would only be able to traverse ~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 ~11-year solar cycle) time scales.
## 5.2. Muon rate changes due to groundwater variations

The temporal variations shown in Figure 5 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 Marschall et al. [2004], 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 [Goldscheider and Drew, 2007]. 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 6: precipitation, stream flow of a nearby river (cf. Figure 1 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 [Vereecken et al., 2016]. 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 $(\Delta 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 6b).

In Figures 6a and 6b, 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 6b 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 (Figure 7). We consider a layer of water-saturated
rock of height $\Delta z$ (m), located at a depth $d$ (m) following the topography. The saturated layer height is given by $\Delta z = w/p$, where $w$ (m) and $p$ correspond to the amount of infiltrated water, and the layer porosity, respectively. We then compute $\Delta L$ (m); the length of the muon trajectory $L$ that lies inside the considered layer. $\Delta L$ thus corresponds to the distance traveled by muons, in a certain axis of observation, inside the water-saturated layer. Figure 7 illustrates the model with an example, where an axis of observation intersects a water-saturated layer of height $\Delta z = 2$ m in two segments, resulting in muons traversing a total distance of $\Delta L = \Delta L_1 + \Delta L_2 \approx 30$ m within the water-saturated layer, where $\Delta L_1$ and $\Delta L_2$ represent the distance travelled by muons in each segment. We define the Amplification Factor $\mathcal{A}$ (mwe·m⁻¹) 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} = \Delta\rho/w$. $\Delta\rho$ is given by $\Delta\rho = p \times \Delta L \times \rho_{\text{water}}$, where $\rho_{\text{water}} = 1 \text{ g·cm}^{-3}$ is the density of water.

In Figure 8 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 × 31 to 61 × 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$ mwe·m⁻¹ in the HD region, indicating a negligible amplification effect, $\mathcal{A} > 1$ mwe·m⁻¹ in the LD region, that is, in the region where significant muon rate variations are observed (Figure 8). We assess the influence of each parameter on the computation of $\mathcal{A}$ by modifying the values of porosity ($N_p = 4$ linearly spaced values from 0.2 to 0.5), infiltrated water ($N_w = 10$ linearly spaced values from 0.1 m to 3 m), and the layer depth ($N_d = 33$ linearly spaced values from 6 m to 102 m). From these $N_p \times N_w \times N_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 9a, 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 9b shows the corresponding standard deviation of $\mathcal{A}$, $\sigma_{\mathcal{A}}$.

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:

$$ \overline{\mathcal{A}}(d) = \frac{\sum_{i=1}^{N_w} \sum_{j=1}^{N_p} \frac{\langle \mathcal{A} \rangle(w_i, p_j, d)}{\sigma_{\mathcal{A}}(w_i, p_j, d)}}{\sum_{i=1}^{N_w} \sum_{j=1}^{N_p} \frac{1}{\sigma_{\mathcal{A}}(w_i, p_j, d)}}. \quad (5) $$

The corresponding standard error is given by:

$$ \sigma_{\overline{\mathcal{A}}}(d) = \sqrt{\frac{1}{\frac{\sum_{i=1}^{N_w} \sum_{j=1}^{N_p} \frac{1}{\sigma_{\mathcal{A}}(w_i, p_j, d)}}{\sum_{i=1}^{N_w} \sum_{j=1}^{N_p} \frac{1}{\sigma_{\mathcal{A}}(w_i, p_j, d)}}}}. \quad (6) $$

The results are shown in Figure 10. 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.

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 6b). Another is that the amplification factor only accounts for the
<figure><img src="image_8.png" /></figure>

**Figure 7.** 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 $\Delta z$ at a depth $d$ is represented as a blue surface. $\Delta L = \Delta L_1 + \Delta L_2$ corresponds to the distance traveled by muons in the water-saturated layer.

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

**Figure 8.** 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 7). 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.

effect of a water-saturated layer parallel to the surface, while water content variations may also occur in more complex geometries. Moreover, local-scale

features, such as preferential flow paths and karst heterogeneity, may influence water distribution in ways not captured by regional-scale satellite data.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Mean value of *A*, ⟨*A*⟩, 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 *A*, σ_A, for the combined HD and LD regions in function of *d* and *w*, and for different porosity values *p*.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Mean Amplification Factor $\bar{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.</figcaption></figure>

**6. Conclusions**

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 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 phenomena.

The sensitivity analyses on the Amplification Factor reveal the large impact of the geometry given by the trajectories of detected muons and the topography. We thus point to the possibility of enhancing the sensitivity of the method in groundwater
storage estimations by adequately setting the exper-
iment 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 muogra-
phy in hydrogeophysics, potentially shedding light
not only on the groundwater storage variations, but
also on assessing the depth of the water table.

Declaration of interests

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

Funding

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

Acknowledgments

We extend our gratitude to Thierry Theurillat and
Senecio Schefer for their technical and logistical as-
sistance 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 https:
//www.ecmwf.int/. The Doubs river stream flow is
available from https://www.hydrodaten.admin.ch/.
Muon data used for all calculations are displayed in
figures and are available in the supplementary mate-
rial. This study uses perceptually uniform color maps
to avoid visual distortion of the data and exclusion
of readers with colour-vision deficiencies [Crameri
et al., 2020].

Supplementary data

Supporting information for this article is available on
the journal's website under https://doi.org/10.5802/
crgeos.277 or from the author.

**References**

Alley, W. M., Healy, R. W., LaBaugh, J. W., and Reilly, T. E. (2002). Flow and storage in groundwater systems. *Science*, 296(5575), 1985–1990.

Bajou, R., Rosas-Carbajal, M., Tonazzo, A., and Marteau, J. (2023). High-resolution structural imaging of volcanoes using improved muon tracking. *Geophys. J. Int.*, 235(2), 1138–1149.

Barrett, P. H., Bollinger, L. M., Cocconi, G., Eisenberg, Y., and Greisen, K. (1952). Interpretation of cosmic-ray measurements far underground. *Rev. Mod. Phys.*, 24(3), 133–178.

Beni, T., Borselli, D., Bonechi, L., et al. (2023). Transmission-based muography for ore bodies prospecting: A case study from a Skarn Complex in Italy. *Nat. Resour. Res.*, 32(4), 1529–1547.

Bevington, P. R. and Robinson, D. K. (2003). *Data Reduction and Error Analysis*. McGraw Hill, New York.

Blazevic, L. A., Bodet, L., Pasquet, S., Linde, N., Jougnot, D., and Longuevergne, L. (2020). Time-lapse seismic and electrical monitoring of the vadose zone during a controlled infiltration experiment at the ploemeur hydrological observatory, France. *Water*, 12(5), article no. 1230.

Bonechi, L., D'Alessandro, R., and Giammanco, A. (2020). Atmospheric muons as an imaging tool. *Rev. Phys.*, 5, article no. 100038.

Bonomi, G., Checchia, P., D'Errico, M., Pagano, D., and Saracino, G. (2020). Applications of cosmic-ray muons. *Prog. Part. Nucl. Phys.*, 112, article no. 103768.

Borselli, D., Beni, T., Bonechi, L., et al. (2022). Three-dimensional muon imaging of cavities inside the temperino mine (Italy). *Sci. Rep.*, 12(1), article no. 22329.

Bossart, P., Bernier, F., Birkholzer, J., et al. (2018). Mont Terri rock laboratory, 20 years of research: introduction, site characteristics and overview of experiments. In *Mont Terri Rock Laboratory, 20 Years*, pages 3–22. Springer, Switzerland.

Castellazzi, P., Martel, R., Galloway, D. L., Longuevergne, L., and Rivera, A. (2016). Assessing groundwater depletion and dynamics using grace and insar: Potential and limitations. *Groundwater*, 54(6), 768–780.

Cimmino, L., Baccani, G., Noli, P., et al. (2019). 3D muography for the search of hidden cavities. *Sci.*
*Rep.*, 9(1), article no. 2974.

Cohu, A., Chevalier, A., Nechyporuk, O., Franzen, A., Sauerwald, J., Ianigro, J.-C., and Marteau, J. (2023). First 3D reconstruction of a blast furnace using muography. *J. Instrument.*, 18(07), article no. P07004.

Cohu, A., Tramontini, M., Chevalier, A., Ianigro, J.-C., and Marteau, J. (2022). Atmospheric and geodesic controls of muon rates: A numerical study for muography applications. *Instruments*, 6(3), article no. 24.

Crameri, F., Shephard, G. E., and Heron, P. J. (2020). The misuse of colour in science communication. *Nat. Commun.*, 11(1), article no. 5444.

Crossley, D., Hinderer, J., and Riccardi, U. (2013). The measurement of surface gravity. *Rep. Prog. Phys.*, 76(4), article no. 046101.

Evans, J. G., Ward, H. C., Blake, J. R., et al. (2016). Soil water content in southern england derived from a cosmic-ray soil moisture observing system - cosmos-uk. *Hydrol. Process.*, 30(26), 4987-4999.

Gaillardet, J., Braud, I., Hankard, F., et al. (2018). Ozcar: The French network of critical zone observatories. *Vadose Zone J.*, 17(1), 1-24.

Gaisser, T. K., Engel, R., and Resconi, E. (2016). *Cosmic Rays and Particle Physics*. Cambridge University Press, UK.

Gleeson, T., Moosdorf, N., Hartmann, J., and van Beek, L. V. (2014). A glimpse beneath Earth's surface: Global hydrogeology maps (glhymps) of permeability and porosity. *Geophys. Res. Lett.*, 41(11), 3891-3898.

Goldscheider, N. and Drew, D. (2007). *Methods in Karst Hydrogeology: IAH: International Contributions to Hydrogeology*, volume 26. CRC Press, Boca Raton.

Grashorn, E., De Jong, J., Goodman, M., Habig, A., Marshak, M., Mufson, S., Osprey, S., and Schreiner, P. (2010). The atmospheric charged kaon/pion ratio using seasonal variation methods. *Astropart. Phys.*, 33(3), 140-145.

Grobbe, N., Mordret, A., Barde-Cabusson, S., et al. (2021). A multi-hydrogeophysical study of a watershed at kaiwi coast (OAhu, Hawaii), using seismic ambient noise surface wave tomography and self-potential data. *Water Resour. Res.*, 57(4), article no. e2020WR029057.

Groom, D. E., Mokhov, N. V., and Striganov, S. I. (2001). Muon stopping power and range tables

10 mev-100 tev. *At. Data Nucl. Data Tables*, 78(2), 183-356.

Hamming, R. (1998). *Digital Filters*. Dover Publications, New York.

Heck, D., Knapp, J., Capdevielle, J., et al. (1998). CORSIKA: A Monte Carlo code to simulate extensive air showers. *Report FZKA*, 6019(11), article no. FZKA-6019.

Hermans, T., Goderniaux, P., Jougnot, D., et al. (2023). Advancing measurements and representations of subsurface heterogeneity and dynamic processes: towards 4D hydrogeology. *Hydrol. Earth Syst. Sci.*, 27(1), 255-287.

Hersbach, H., Bell, B., Berrisford, P., et al. (2020). The ERA5 global reanalysis. *Q. J. R. Meteorol. Soc.*, 146(730), 1999-2049.

Hörandel, J. R. (2003). On the knee in the energy spectrum of cosmic rays. *Astroparticle Phys.*, 19(2), 193-220.

Jourde, K., Gibert, D., Marteau, J., de Bremond d'Ars, J., Gardien, S., Girerd, C., and Ianigro, J.-C. (2016a). Monitoring temporal opacity fluctuations of large structures with muon radiography: a calibration experiment using a water tower. *Sci. Rep.*, 6, article no. 23054.

Jourde, K., Gibert, D., Marteau, J., de Bremond d'Ars, J., and Komorowski, J.-C. (2016b). Muon dynamic radiography of density changes induced by hydrothermal activity at the La Soufrière of Guadeloupe volcano. *Sci. Rep.*, 6, article no. 33406.

Köhli, M., Schrön, M., Zreda, M., Schmidt, U., Dietrich, P., and Zacharias, S. (2015). Footprint characteristics revised for field-scale soil moisture monitoring with cosmic-ray neutrons. *Water Resour. Res.*, 51(7), 5772-5790.

Le Gonidec, Y., Rosas-Carbajal, M., Bremond d'Ars, J. D., Carlus, B., Ianigro, J.-C., Kergosien, B., Marteau, J., and Gibert, D. (2019). Abrupt changes of hydrothermal activity in a lava dome detected by combined seismic and muon monitoring. *Sci. Rep.*, 9(1), article no. 3079.

Leblanc, M. J., Tregoning, P., Ramillien, G., Tweed, S. O., and Fakes, A. (2009). Basin-scale, integrated observations of the early 21st century multiyear drought in southeast Australia. *Water Resour. Res.*, 45(4), article no. W04408.

Leo, W. R. (1994). *Techniques for Nuclear and Particle Physics Experiments: A How-to Approach*. Springer Science & Business Media, Germany.
Lesparre, N., Cabrera, J., and Marteau, J. (2016). 3-D density imaging with muon flux measurements from underground galleries. *Geophys. J. Int.*, 208(3), 1579-1591.

Lesparre, N., Gibert, D., Marteau, J., Déclais, Y., Carbonne, D., and Galichet, E. (2010). Geophysical muon imaging: feasibility and limits. *Geophys. J. Int.*, 183(3), 1348-1361.

Longuevergne, L., Scanlon, B. R., and Wilson, C. R. (2010). Grace hydrological estimates for small basins: Evaluating processing approaches on the high plains aquifer, USA. *Water Resour. Res.*, 46(11), article no. W11517.

Marschall, P., Croisé, J., Schlichenrieder, L., Boisson, J. Y., Vogel, P., and Yamamoto, S. (2004). Synthesis of hydrogeological investigations at the Mont Terri site (phases 1 to 5). In Heitzmann, P., editor, *Mont Terri Project – Hydrogeological Synthesis, Omotonic Flow*. Office fédéral des eaux et de la géologie OFEG, Switzerland.

Marteau, J., de Bremond d'Ars, J., Gibert, D., Jourde, K., Gardien, S., Girerd, C., and Ianigro, J.-C. (2014). Implementation of sub-nanosecond time-to-digital convertor in field-programmable gate array: applications to time-of-flight analysis in muon radiography. *Meas. Sci. Technol.*, 25(3), article no. 035101.

Marteau, J., de Bremond d'Ars, J., Gibert, D., Jourde, K., Ianigro, J.-C., and Carlus, B. (2017). DIAPHANE: Muon tomography applied to volcanoes, civil engineering, archaelogy. *J. Instrument.*, 12(02), article no. C02008.

Marteau, J., Gibert, D., Lesparre, N., Nicollin, F., Noli, P., and Giacoppo, F. (2012). Muons tomography applied to geosciences and volcanology. *Nucl. Instr. Meth. Phys. Res. Sec. A*, 695, 23-28.

McGarr, J. T., Wallace, C. D., Ntarlagiannis, D., Sturmer, D. M., and Soltanian, M. R. (2021). Geophysical mapping of hyporheic processes controlled by sedimentary architecture within compound bar deposits. *Hydrol. Process.*, 35(9), article no. e14358.

Morishima, K., Kuno, M., Nishio, A., et al. (2017). Discovery of a big void in khufu's pyramid by observation of cosmic-ray muons. *Nature*, 552(7685), 386-390.

Muñoz-Sabater, J., Dutra, E., Agustí-Panareda, A., et al. (2021). ERA5-Land: a state-of-the-art global reanalysis dataset for land applications. *Earth Syst.

Sci. Data, 13(9), 4349-4383.

Nishiyama, R., Ariga, A., Ariga, T., et al. (2019). Bedrock sculpting under an active alpine glacier revealed from cosmic-ray muon radiography. *Sci. Rep.*, 9(1), article no. 6970.

Nussbaum, C., Kloppenburg, A., Caër, T., and Bossart, P. (2017). Tectonic evolution around the Mont Terri rock laboratory, northwestern Swiss Jura: constraints from kinematic forward modelling. *Swiss J. Geosci.*, 110(1), 39-66.

Oláh, L., Gallo, G., Hamar, G., et al. (2023). Muon imaging of volcanic conduit explains link between eruption frequency and ground deformation. *Geophys. Res. Lett.*, 50(2), article no. e2022GL101170.

Oláh, L., Tanaka, H. K., and Hamar, G. (2021). Muographic monitoring of hydrogeomorphic changes induced by post-eruptive lahars and erosion of Sakurajima volcano. *Sci. Rep.*, 11(1), article no. 17729.

Osprey, S., Barnett, J., Smith, J., et al. (2009). Suden stratospheric warmings seen in MINOS deep underground muon data. *Geophys. Res. Lett.*, 36(5), article no. L05809.

Rosas-Carbajal, M., Jourde, K., Marteau, J., Deroussi, S., Komorowski, J.-C., and Gibert, D. (2017). Three-dimensional density structure of La Soufrière de Guadeloupe lava dome from simultaneous muon radiographies and gravity data. *Geophys. Res. Lett.*, 44(13), 6743-6751.

Sagisaka, S. (1986). Atmospheric effects on cosmic-ray muon intensities at deep underground depths. *Il Nuovo Cimento C*, 9, 809-828.

Schmidt, L. and Rempe, D. (2020). Quantifying dynamic water storage in unsaturated bedrock with borehole nuclear magnetic resonance. *Geophys. Res. Lett.*, 47(22), article no. e2020GL089600.

Schmidt, R., Flechtner, F., Meyer, U., Neumayer, K.-H., Dahle, C., König, R., and Kusche, J. (2008). Hydrological signals observed by the GRACE satellites. *Surv. Geophys.*, 29(4), 319-334.

Spurio, M. (2014). *Particles and Astrophysics*. Springer, Switzerland.

Syed, T. H., Famiglietti, J. S., Rodell, M., Chen, J., and Wilson, C. R. (2008). Analysis of terrestrial water storage changes from grace and gldas. *Water Resour. Res.*, 44(2), article no. W02433.

Taketa, A., Nishiyama, R., Yamamoto, K., and Iguchi, M. (2022). Radiography using cosmic-ray electro-
magnetic showers and its application in hydrology.
*Sci. Rep.*, 12(1), article no. 20395.

Tramontini, M., Rosas-Carbajal, M., Nussbaum, C., Gibert, D., and Marteau, J. (2019). Middle-atmosphere dynamics observed with a portable muon detector. *Earth Space Sci.*, 6(10), 1865–1876.

Vereecken, H., Schnepf, A., Hopmans, J., et al. (2016).
Modeling soil processes: Review, key challenges,
and new perspectives. *Vadose Zone J.*, 15(5), 1–57.

Workman, R. L., Burkert, V. D., Crede, V., et al. (2022).
Review of particle physics (particle data group).
*Prog. Theor. Exp. Phys.*, 2022(8), article no. 083C01.