# Monitoring unsaturated water flow using magnetic resonance soundings

Anatoly Legchenko, Jean-Michel Baltassat, Mohamad Abbas, Arnaud Isch, Nadia Amraoui, Jacques Deparis, Pauline Kessouri, Clémence Ryckebusch, Mohamed Azaroual

► To cite this version:

Anatoly Legchenko, Jean-Michel Baltassat, Mohamad Abbas, Arnaud Isch, Nadia Amraoui, et al. Monitoring unsaturated water flow using magnetic resonance soundings. Journal of Hydrology, 2022, 612, pp.128318. 10.1016/j.jhydrol.2022.128318. insu-03762895

HAL Id: insu-03762895

https://insu.hal.science/insu-03762895

Submitted on 29 Aug 2022

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.
Journal Pre-proofs

Research papers

Monitoring unsaturated water flow using magnetic resonance soundings

Legchenko Anatoly, Baltassat Jean-Michel, Abbas Mohamad, Isch Arnaud, Amraoui Nadia, Azaroual Mohamed, Deparis Jacques, Ryckebusch Clemence, Kessouri Pauline

PII: S0022-1694(22)00890-3

DOI: https://doi.org/10.1016/j.jhydrol.2022.128318

Reference: HYDROL 128318

To appear in: *Journal of Hydrology*

Received Date: 27 March 2022

Revised Date: 9 July 2022

Accepted Date: 3 August 2022

Please cite this article as: Anatoly, L., Jean-Michel, B., Mohamad, A., Arnaud, I., Nadia, A., Mohamed, A., Jacques, D., Clemence, R., Pauline, K., Monitoring unsaturated water flow using magnetic resonance soundings, *Journal of Hydrology* (2022), doi: https://doi.org/10.1016/j.jhydrol.2022.128318

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

© 2022 Published by Elsevier B.V.

<figure><img src="image_1.png" /></figure>
Monitoring unsaturated water flow using magnetic resonance soundings

¹Legchenko Anatoly, ²Baltassat Jean-Michel, ³Abbas Mohamad, ⁴Isch Arnaud, ⁵Amraoui Nadia, ⁶Azaroual Mohamed, ⁷Deparis Jacques, ⁸Ryckebusch Clemence, ⁹Kessouri Pauline

¹ Univ. Grenoble Alps, IRD, IGE, Grenoble, France, E-mail: [anatoli.legtchenko@ird](mailto:anatoli.legtchenko@ird)

² BRGM, Orléans, France, E-mail: jm.baltassat@brgm.fr

³ Univ. Orléans, CNRS, ISTO, UMR 7327, F-45071, Orléans, France, E-mail: [mohamad.abbas@cnrs-orleans.fr](mailto:mohamad.abbas@cnrs-orleans.fr)

⁴ Univ. Orléans, CNRS, ISTO, UMR 7327, F-45071, Orléans, France, E-mail: [arnaud.isch@gmail.com](mailto:arnaud.isch@gmail.com)

5 BRGM, Orléans, France, E-mail: n.amraoui@brgm.fr

⁶ Univ. Orléans, CNRS, BRGM, ISTO, UMR 7327, F-45071, Orléans, France, E-mail : [mohamed.azaroual@cnrs-orleans.fr](mailto:mohamed.azaroual@cnrs-orleans.fr)

⁷ BRGM, Orléans, France, E-mail: j.deparis@brgm.fr

8 BRGM, ISTO, Orléans, France, E-mail: c.ryckebusch@brgm.fr

⁹ BRGM, Orléans, France, E-mail: p.kessouri@brgm.fr

**Key words:** MRS, time-lapse, unsaturated water flow, limestone, Villamblain, O-ZNS observatory

Abstract
We present an innovative hydrogeophysical approach for non-invasive quantification of the unsaturated water flow. For the water content measurements, we apply the Magnetic Resonance Soundings (MRS) method in the time-lapse mode. For inversion of MRS measurements, we approximate the subsurface by a horizontally stratified media. Laterally, the MRS estimated water content is averaged over the area below the measuring setup. During this study, we used a square loop of 75 m side-length. The time-lapse inversion allows visualizing variations in the water content in the vadose zone down to twenty meters. We tested our approach at the Villamblain test site (France). For that, we carried out the MRS monitoring in 1999-2000 and 2020-2021. During these two observation windows, MRS shows seasonal variations of the water content but also a strong dependence of the water content on the rainfall that was different between two data sets. Comparison of MRS results with the rainfall records shows a good correspondence in-between. For measuring and processing MRS data, we used commercially available MRS instrument and interpretation software.

# 1. Introduction

Numerical modeling of hydrological processes requires measuring various input variables, such as hydro-meteorological data and mass and heat exchanges between the topsoil and the atmosphere. In aquifers and shallow soils, these measurements are well established and frequently used. The vadose zone (VZ) plays an important role in the transmission of water and pollutants from the surface to groundwater. Therefore, an understanding of the unsaturated flow is crucial for the sustainable management of water resources and agricultural activities. However, quantification of the unsaturated flow in the VZ is a more difficult task than that in aquifers and the corresponding methods are less developed.

In the VZ, hydrogeologists use centimetric scale measurements of the hydraulic properties and water content by single-point techniques, such as neutron probes, tensiometers, and time domain reflectometry (TDR), (e.g. Isch et al., 2019; Skierucha et al., 2012; Verhoef et al., 2006;
Evett, 2003; Robinson et al., 2003). However, these techniques are intrusive, applied in shallow
depths or a small area around borehole body, and result in local information at only several
points that may not represent the spatial water content distribution at the field scale. Non-
invasive geophysical techniques represent alternative solutions for monitoring water content
and water flow within the VZ (e.g. Binley et al., 2015; Vereecken et al., 2006; Rubin and
Hubbard, 2005).

Geophysical techniques have shown a promising potential in providing large-scale and relatively high-resolution images of hydrogeological parameters and processes (e.g. Binley et al., 2015; 2010). Different authors report surface, borehole, and cross-hole geophysical techniques applied to hydrogeological applications. The most frequently used methods are the ground penetrating radar (GPR) (e.g. Dafflon et al., 2011; Lunt et al., 2005), seismic techniques (e.g. Blazevic et al., 2020), electrical resistivity (ER) techniques (e.g. Johnson et al., 2015; Robinson et al., 2008), gravity (e.g. Eppelbaum, 2019), self-potential (e.g. Abbas et al., 2017; Eppelbaum, 2019; Jardani et al., 2012), and induced polarization (e.g. Johnson et al., 2010). Time-lapse geophysical measurements allow monitoring temporal and spatial variations of the water content following natural or artificial infiltration processes (e.g., Uhlemann et al., 2017). Comparison of the geophysical datasets with the local hydrogeological parameters at specific locations allows getting petrophysical relationships and improves imaging preferential water flow pathways in both continuous and fractured geological formations (e.g., De Jong et al., 2020; Gance et al., 2016; Wehrer and Slater, 2015; Steelman and Endres, 2011; Brunet et al., 2010; Cassiani et al., 2009b; Looms et al., 2008; Linde et al., 2006; Kowalsky et al., 2005; Huisman et al., 2003; Binley et al., 2002a; Hubbard et al., 2001a).

However, estimating hydraulic properties by the aforementioned geophysical methods is not
a trivial process, and the accuracy of this approach is dependent on the petrophysical
relationships that link the measured physical parameters to the subsurface properties of interest
(e.g., Topp et al., 1980; Archie, 1942). This introduced difficulties reflected in the need of
expensive and time-consuming laboratory calibration of the empirical relationships that are
generally only suitable for homogeneous structures, and affected by not only the water content,
but to a higher extent by geological characteristics such as the porosity and tortuosity (Glover,
2016).

The Magnetic Resonance Sounding (MRS) is a geophysical method directly and selectively sensitive to groundwater. It is a large scale and time-saving tool for hydrogeological applications, which is used for characterizing aquifers in terms of the specific yield and hydraulic conductivity at the field scale (e.g. Legchenko et al., 2021; Valois et al., 2018; Vilhelmsen et al., 2014; Legchenko et al., 2013; Vouillamoz et al., 2012; Nielsen et al., 2011; Boucher et al., 2009; Chalikakis et al., 2008; Lubczynski and Roy, 2003; Legchenko et al., 2002). The MRS method relies on the capability of hydrogen protons contained in water molecules to generate a decaying oscillating magnetic field in response to electromagnetic pulses generated by the surface loop. The strength of the response is related to the number of excited protons, and therefore, directly proportional to the water content. Moreover, the decay time is proportional to the pore size, and therefore, can infer information about the hydraulic conductivity of the medium (e.g. Hertrich et al., 2007; Legchenko et al., 2002; Legchenko and Shushakov, 1998; Schirov et al., 1991). The MRS is an emerging technique, which was used for several hydrogeological applications including the VZ (e.g. Legchenko et al., 2020; Falzone and Keating, 2016; Legchenko et al., 2014; Vouillamoz et al., 2014; Herckenrath et al., 2012; Walsh et al., 2012; Costabel and Yaramanci, 2011; Descloitres et al., 2008; Legchenko et al., 2008).

The goal of our study was to investigate efficiency of the MRS method applied to
quantification of the water content in the unsaturated geological formation. In this paper, we
report results of the water content monitoring using time-lapse MRS measurements in a highly
heterogeneous VZ characterized by altered, fractured and karstified limestone formations. We show that, in non-magnetic rocks like chalk and limestone, the MRS can measure seasonal variations of water content in the VZ at a depth between 0 and 20 m.

## 2. Geology and hydrogeology of the experimental site

The experimental site is located in France, near Villamblain, 30 km northwest of Orleans (Fig. 1). This area is a sub-catchment of the Beauce aquifer system (9700 km²). The Beauce aquifer is the largest drinking water reserve in France and cover an area which is characterized by intensive agricultural activities and irrigation mainly in summer (Aldana et al., 2021; Lejars et al., 2012; Flipo et al., 2012). The regional hydraulic gradient is about 0.1% with the mean annual recharge estimated at 110 mm (Schnebelen et al., 1999). Calcareous soils and cryoturbated materials cover about 48% of this area. The shallow subsurface is composed of highly permeable soils with a vertical soil water flow regime (Bruand et al., 1997). Previous studies report an absence of run-off in Villamblain (Schnebelen et al., 1999).
<figure><img src="image_2.png" /></figure>

Figure 1. Location of the MRS test site in Villamblain with respect to the Beauce aquifer and the Seine and Loire basins. Blue circles report main agriculture boreholes locations. Poiseaux and Bricy stations provide respectively the rainfall data and GWL. The map of the Beauce aquifer is modified from Flipo et al., (2012).

Fig. (2) shows the annual rainfall (blue bars) and the groundwater level (GWL) (red line) recorded between 1965 and 2021. The Orléans - Bricy (Loiret - France) meteorological station (20 km from the study site) provides the rainfall data and the Poiseaux station (4 km from the
study site) the GWL (Fig. 1). The GWL annual dynamics overlaid by multiannual trends correlates with the rainfall data, but we also observe a discrepancy in-between. The scale difference may explain this discrepancy: while the rainfall is a local event, the GWL depends on many factors at the watershed scale.

<figure><img src="image_3.png" /><figcaption>Figure 2. The annual rainfall (blue bars) and the groundwater level (GWL) (red line) collected between 1965 and 2021.</figcaption></figure>

The MRS test site is located near the O-ZNS observatory (1.5 km) with a few boreholes drilled around (https://plateformes-pivots.eu/o-zns/). Boreholes show a highly heterogeneous geological structure characterized by altered, fractured and karstified limestone facies, along with macro and micro porosity (Aldana et al., 2021). In the VZ down to 20 m, the subsurface comprises three main lithological units with significant vertical and lateral heterogeneities (Fig. 3a):
1) Heterogeneous top soil layer with a thickness of approximately one meter and characterized by a silt loam texture (IUSS Working Group WRB, 2015).

2) Highly heterogeneous, incoherent and marly limestone layer (2-8 m deep), cryoturbated in its upper part during the Quaternary (Ould Mohamed and Bruand, 1994), (Fig. 3a). Fig. (3b) shows strong variations of the saturated hydraulic conductivity within this formation composed of cemented aggregates with a high content of silt and clay alternate with unconsolidated coarse sand and gravel.

3) Hard limestone rock formation (8-20 m) composed of a microcrystalline texture having high carbonate content (>90%) and high bulk density values (Aldana et al., 2021). The formation comprises weathering patterns, fissures, matrix alterations, vugs, open fractures, and karst networks intensifying below fifteen meters. The altered limestone units observed between eight and twenty meters deep are interrupted by well-bedded massive limestone banks.

<figure><img src="image_4.png" /></figure>
Figure 3. (a) The geological cross-section observed in the O-ZNS site; (b) physical properties measured in rock cores (granulometry, calcium carbonate percentage (red curve with dots), saturated volumetric water content (red curve) and density (black curve).

An open pit constructed for this laboratory project allows visual observation of the shallow geological structure in the first 17 meters.

<figure><img src="image_5.png" /><figcaption>Figure 4. Photos of the limestone formation observed during construction of a 4 m diameter and 17 m deep pit at the O-ZNS Observatory.</figcaption></figure>

# 3. Background

## 3.1. Forward modeling and inversion
The MRS method is a geophysical application of the Nuclear Magnetic Resonance (NMR) phenomenon in the Earth's magnetic field. For measuring, a wire loop on the soil surface is energized by a pulse of oscillating current. The frequency of the current is equal to the Larmor frequency ($\omega_0 = \gamma B_0$), where $B_0$ is the Earth's magnetic field and $\gamma$ is the gyromagnetic ratio.

The Larmor frequency is the signature of the hydrogen atoms in groundwater, which renders the MRS method selectively sensitive to groundwater. After the pulse cutoff, groundwater generates an exponentially decaying electromagnetic field that induces the corresponding voltage (NMR signal) in the measuring loop called. The pulse is characterized by the pulse moment $q = I_0\tau$, where moment $I_0$ and $\tau$ are the current amplitude and the pulse duration. The initial amplitude of the MRS signal versus pulse moment $e(q)$ is proportional to the investigated volume $V$ and the volumetric water content $w$. It can be computed using the following integral equation (Legchenko et Valla, 2002):

$$e(q) = \frac{\omega_0}{I_0} \int_V B_1 M_0 \sin\left(\frac{\gamma B_1 q}{2I_0}\right) w d\nu, \quad (1)$$

where $B_1$ is the active component of the loop magnetic field, $M_0$ is the macroscopic nuclear magnetization and $\gamma$ is the gyromagnetic ratio for protons.

We consider the 1-D case and approximate Eq. (1) by a matrix equation:

$$A\mathbf{w} = \mathbf{e}. \quad (2)$$

The matrix $A = [a_{i,j}]$ represents the kernel of the integral equation, the vectors $\mathbf{w} = [w_j]$ and $\mathbf{e} = [e_i]$ correspond to the water content and the amplitude of the MRS signal respectively. $J$ is the number of layers in the inverse model and $I$ is the number of pulse moments ($i=1,2,...,I$; and $j=1,2,...,J$). When processing one individual sounding, the solution of equation Eq. (2) is a vertical distribution of the volumetric water content $w(z)$ (Legchenko and Shushakov, 1998).

For the inversion, we use the Tikhonov regularization method and minimize the Tikhonov functional:
$$ \left\{ \begin{array}{l} \| \mathbf{A} \mathbf{w} - \mathbf{e} \| + \alpha \left\| \frac{\mathrm{d} \mathbf{w}}{\mathrm{d} z} \right\| = \min, \\ \| \mathbf{A} \mathbf{w} - \mathbf{e} \| \le \varepsilon \end{array} \right. \qquad (3) $$

where $\varepsilon$ is an estimate of the data error and $\alpha$ is the smoothing factor.

For investigating the uncertainty in the MRS inversion using the singular value decomposition (SVD) (Aster et al., 2011), we present the matrix $\mathbf{A}$ as a product of three orthogonal matrices: $\mathbf{U}$, $\mathbf{V}$ and $\mathbf{S}$ (Aster et al., 2005):

$$ \mathbf{A} = \mathbf{U} \mathbf{S} \mathbf{V}^T. \qquad (4) $$

The uncertainty caused by experimental noise is:

$$ \Delta \mathbf{w} = \sigma \sqrt{\mathrm{diag}(\mathbf{V} \mathbf{F} \mathbf{S}^{-2} \mathbf{V}^T)}, \qquad (5) $$

were $\sigma^2$ are independent and identically distributed normal data errors and the filter factor $\mathbf{F}$ is a diagonal matrix representing regularization.

Inversion provides the volumetric water content $w_j$ and the thickness $\Delta z_j$ of each layer in the MRS inverse model and allows estimating the equivalent water column:

$$ H_{\mathrm{MRS}} = \sum_j (w_j \Delta z_j). \qquad (6) $$

The water column is a more stable parameter than $w_j$ and $\Delta z_j$ computed separately (Legchenko et al., 2004).

The time-lapse inversion comprises solution of the matrix equation:

$$ \tilde{\mathbf{A}} \tilde{\mathbf{w}} = \tilde{\mathbf{e}}. \qquad (7) $$

In Eq. (7), the data set $\tilde{\mathbf{e}}$ is composed of $K$ individual soundings separated by time intervals $\Delta t_k$ ($k=1,2,...,K$):

$$ \tilde{\mathbf{e}} = \begin{bmatrix} \mathbf{e}_1 \\ \vdots \\ \mathbf{e}_K \end{bmatrix}. \qquad (8) $$

The matrix $\tilde{\mathbf{A}}$ for the time-lapse inversion includes $K$ matrixes $\mathbf{A}_k$ corresponding to each individual sounding and the zero-matrixes:
$$ \tilde{A} = \begin{bmatrix} A_1 & 0 & 0 \\ 0 & \ddots & 0 \\ 0 & 0 & A_K \end{bmatrix}. \qquad (9) $$

The volumetric water content is a two-dimensional function $w(z,t)$ approximated by the solution vector:

$$ \tilde{w} = \begin{bmatrix} w_1 \\ \vdots \\ w_K \end{bmatrix}. \qquad (10) $$

We carry out the time-lapse inversion using the Tikhonov regularization method (Legchenko et al., 2020) and assume the smoothness of the solution ($w(z,t)$) versus both the depth ($z$) and the time ($t$).

$$ \left\{ \begin{array}{l} \| \tilde{A} \tilde{w} - \tilde{e} \| + \alpha_z \| \frac{d \tilde{w}}{dz} \| + \alpha_t \| \frac{d \tilde{w}}{dt} \| = \min. \\ \| \tilde{A} \tilde{w} - \tilde{e} \| \le \tilde{\epsilon} \end{array} \right. \qquad (11) $$

Within the aquifer, rock formations are always saturated, and the water content depends on the rock porosity. It does not vary in time. In this case, measurements of the GWL provide necessary hydrological data about groundwater dynamics.

In the VZ, the GWL monitoring does not provide complete hydrological information and the solution $w(z,t)$ provided by the monitoring of the VZ shows seasonal variations in the water content. It allows estimating the water storage variations by computing the equivalent water column versus time.

## 3.2. Data acquisition

Twice, in 1999 - 2000 and in 2020 - 2021, the BRGM (Bureau de Recherches Géologiques et Minières) carried out a one-year MRS monitoring of the water content at the Villamblain test site (Fig. 1). The MRS data set comprises results of MRS measurements performed with a coincident 75×75 m² square loop always set at the same position (33 soundings in 1999 - 2000 and 10 soundings in 2020 - 2021). The Larmor frequency varied around 2011 Hz in 1999 and around 2039 Hz in 2020. Between the 26th of April 1999 and the 15th of March 2000, the
NUMIS instrument (IRIS Instruments, France) was used; and between the 23rd of July 2020 and 30th of June 2021, the NUMISpoly system (IRIS Instruments, France) was used. MRS measurements provided the volumetric water content in the subsurface with a time step from two weeks to one-month.

For monitoring, the accuracy of MRS measurements is an important issue. The way of verifying the accuracy between different instruments could be a comparison of measurements performed under similar conditions. Fig. (5) shows the annual rainfall recorded before and during monitoring. The records show close values in the years preceding our monitoring (1998 and 2019 respectively). Between the 1st of January and 31st of July 1999, the cumulative rainfall was 416 mm, and for the same period of 2020, it was 307 mm. This showed only 100 mm difference in the rainfall for the 1.5 years before the MRS measurements (considering both data sets, 1999 and 2020). We consider the hydrologic conditions before the monitoring comparable and we expect similar MRS signals

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

Figure 5. The annual rainfall recorded before the monitoring time (1998 and 2019: 600 and 607 mm, respectively) and that corresponding to the monitoring time (1999 and 2020: 837 and 535 mm, respectively).
Fig. (6a) shows amplitude of the MRS signal plotted versus pulse moment. We compare two consecutive soundings performed with the NUMIS instrument on the 19th of July 1999 with one sounding performed using the NUMIS^poly^ instrument on the 23rd of July 2020. This example confirms that the MRS signals measured with the two systems are consistent as expected in the above section considering the limited discrepancies of rainfall in 1998 versus 2019 shown in Fig. (5). In 1999 and in 2020, the water content in the subsurface may be slightly different, and this may explain the minor difference in the amplitudes measured in July 1999 and in July 2020.

Fig. (6b) presents MRS signals with the minimal and the maximal amplitudes observed during our work. The difference is significant and the experimental errors cannot be the reason. Thus, a large rainfall in 1999 (837 mm) in comparison with 2020 (535 mm) is consistent with the observed large difference in MRS amplitudes.

<figure><img src="image_7.png" /><figcaption>Figure 6. An example of the MRS signals versus pulse moment observed during our study: a) measurements under similar hydro-meteorological conditions; b) two soundings showing the maximum and the minimum MRS signals. Solid lines show the theoretical signal computed</figcaption></figure>
after the corresponding inverse models. For all these soundings, the noise level was similar
and it is shown only, for one sounding per graph, by the corresponding color in each graph.

**4. Results**

Figs. (7a and 7b) show the water content and the relaxation time $T_2^*$ versus depth derived from the MRS inversion. The water content profiles correspond to the minimal and the maximal MRS signals shown in Fig. (6b). Inversion fits experimental data with the misfit being smaller than the noise estimate (Fig. 6b, solid lines). The diagonal elements of the model resolution matrix (Fig. 7c) confirm that the inverse models are well resolved with the error bars estimated with the SVD. Note that the SVD-estimate corresponds to the uncertainty caused by the experimental noise and not to the accuracy of the water content measured with MRS (Aster et al., 2011; Legchenko, 2021). Fig.(7d) shows the resistivity profile and the lithological log of two nearby boreholes (X0119 and X0107 – see Fig. 1). Fig. (7a) shows that MRS does not point to the GWL. Water in the capillary fringe and limited resolution of the MRS inversion explain the smooth increase of the water content toward the GWL. The electrical resistivity of 10 ohm.m observed in the first 10 m is consistent with the powdery limestone with high clay content observed in the boreholes and with the geological description made by the O-ZNS observatory (Fig. 3). MRS cannot measure short signals corresponding to bound water and consequently does not see water in clay and the waterless interval shown by MRS corroborates well with other observations. Below 25 m, MRS shows a lower water content corresponding to a smaller specific yield of the Beauce aquifer related with the Gatinais molasses and other marly layers observed in the lower limestone of X0107 borehole below 25 m. Larger MRS signal corresponds to more water in the VZ and a long relaxation time (Fig. 7b) points to additional water stored preferably in large pores and in karst voids.
<figure><img src="image_8.png" /><figcaption>Figure 7. a) The water content profiles corresponding to the minimum (black) and the maximum (red) MRS signals; b) corresponding relaxation times $T_2^*$; c) corresponding diagonal elements of the model resolution matrix; d) the resistivity profile provided by electrical measurements and the lithological log of the boreholes X0119 (03622X0119/FP3-25) and X0107 (03622X017/F) located near MRS station (see Fig. 1). GWL at both dates are within 30 cm of each other.</figcaption></figure>

Time-lapse inversion provides the water content versus time recorded with MRS in 1999-2000 (Fig. 8a) and in 2020-2021 (Fig. 8b). Figs. (8c) and (8d) show the daily rainfall recorded during the monitoring period.
<figure><img src="image_9.png" /></figure>

Figure 8. MRS inverse models of the water content versus time (color scale) computed down to 22 m and a daily rainfall corresponding to the monitoring time: a,c) 1999-2000 data set; b,d) 2020-2021 data set.

For validating the inversion, we show the time-lapse MRS data sets (Fig. 9). We normalize the pulse moment values considering the maximum pulse moment over all soundings and, for each sounding, plot the measured amplitude versus normalized pulse moment (red dots). Fig. (9a) shows the 1999-2000 data set and Fig. (9b) the 2020-2021 data set. The theoretical amplitudes computed after the inverse model (black lines) fit well with the measured ones.
<figure><img src="image_10.png" /></figure>

Figure 9. For each sounding, measured amplitude of the MRS signal is plotted versus normalized pulse moment of the corresponding data set and the inversion fit: a) 1999-2000 data set; b) 2020-2021 data set.

Fig. (10) shows the equivalent water column ($H_{MRS}$) versus time estimated with the MRS in 1999-2000 (red circles) and in 2020-2021 (blue crosses). We plot these data starting from the 1st of January. Therefore, we set the zero-time to the 1st of January 1999 and the 1st of January 2020. In this coordinate system, the MRS measurements started on the 26th of April 1999 and on the 23rd of July 2020, the water columns ($H_{MRS}$) are shifted relative to the 1st of January (the zero-time) of the corresponding year.
<figure><img src="image_11.png" /><figcaption>Figure 10. The water column estimated with MRS in 1999-2000 and in 2020-2021 plotted versus time. The zero-time corresponds to the 1st of January 1999 and the 1st of January 2020.</figcaption></figure>

MRS measurements carried out the 19th of July 1999 and the 23rd of July 2020 (Fig. 6a, day 200) reveal about equal $H_{MRS}$ (Fig. 10). Between the 1st of August and the 31st of December, the cumulative rainfall was 421 mm in 1999 and 228 mm in 2020. The meteorological data show that the cumulative reference Penman-Monteith potential evapotranspiration (ETP) was significantly higher in 2020 than in 1999 (1044 mm in 2020 against 804 mm in 1999). This means that the water infiltration was larger in 1999-2000 in comparison with 2020-2021. Thus, the larger $H_{MRS}$ estimated with MRS in 1999 compared to 2020 is consistent with more rainfall water and smaller ETP recorded in 1999. We plot $H_{MRS}$, the monthly rainfall and the GWL versus time considering both data sets (Fig. 11). The increase of the $H_{MRS}$ observed in 1999 follows the rainfall, and the GWL rises with about three months delay after the intensive rainfall. In 2020, the correlation between the GWL and the rainfall, as well as the $H_{MRS}$, is less pronounced. We can partly explain this observation by an insufficient sampling rate of MRS measurements in 2020.
<figure><img src="image_12.png" /></figure>

Figure 11. The equivalent water column ($H_{MRS}$), the monthly rainfall and the GWL versus time: a) 1999 – 2000 data sets; b) 2020 – 2021 data set. In these graphs, the vertical scale for plotting the MRS water column is set different for showing more details in figure 11b.

We use the 1999 data set for verifying correlations between the MRS and other observations. For that, we plot the cumulative rainfall computed between the 1st of September and the 31st of December 1999 versus monitoring time (Fig. 12). The water accumulation ($H_{MRS}$) increases following the cumulative rainfall. The GWL rises with a delay relative to the cumulative rainfall, thus confirming that MRS measures water in the VZ and not the GWL variations. Note that the spike of rainfall in September 1999 (Fig. 11a) does not produce the corresponding spike of the $H_{MRS}$. We explain it by an insufficient sampling rate during our MRS monitoring.
<figure><img src="image_13.png" /><figcaption>Figure 12. The water column estimated with MRS (red), the cumulative rainfall (black) and the depth of the GWL (blue) observed in Villamblain in 1999-2000.</figcaption></figure>

**5. Discussion**

The volumetric water content provided by the time-lapse inversion of the 1999-2000 data (Fig. 8a) allows representing the VZ in four principal layers at the depths of 0-2 m, 2-10 m, 10-16 m and 16-20 m. In the winter of 1999-2000, the MRS data show water accumulation between 9 and 12 m deep. The MRS does not show this accumulation in 2020-2021 (Fig. 8b). This observation is consistent with meteorological data showing smaller rainfall and higher ETP recorded in 2020 (535 and 1044 mm respectively) relative to 1999 (837 and 804 mm respectively).

Both, in 1999-2000 and 2020-2021, the MRS data show a small amount of water in the layer between 2 and 10 m deep. In boreholes, the description of the lithologies encountered between 2 and 8 m deep shows the presence of powdery limestone with high clay and silt content and calcareous sand interbeds (Aldana et al. 2021). The 8-10 m layer represents the transition between these facies and the hard limestone rock with karstic features, containing no or very
little clay. The MRS cannot detect water in clay and we may interpret the MRS data as a low
permeable clay material between 2 and 10 m (Legchenko, 2021). However, this formation may
also correspond to a permeable formation (coarse sand, gravel, karst) that allows a rapid water
transition. When not saturated, this formation has a low storage capacity, which could explain
the MRS results. Indeed, Mazzilli et al. (2020) report the MRS monitoring of a heavy rain event
in a karst environment. The authors observed a large MRS signal only during a few days after
the rain event. This large signal corresponds to an intensive water flow that refills the karst
channels. During our monitoring, we measured with a time step of over two weeks between the
soundings and an insufficient sampling may explain why we do not see water in the Villamblain
karst after a heavy rain event. The MRS results can thus be explained by both hypotheses:

(1) a low permeable powdery limestone with high clay content; or

(2) a very permeable formation, both corresponding well to the measured signals.

However, if the formation (2-8 m) has a very low permeability, then we can expect run-off
events. Knowing that there is no run-off in Villamblain, we interpret the layer between 2 and
10 m as a highly permeable formation with a high clay content in the limestone matrix. This
clay was observed in boreholes and explain a low electrical resistivity revealed by geophysical
measurements. Boreholes do not show developed karst channels in this permeable zone. The
powdery limestone observed in boreholes is intercalated with clay, coarse sand and gravel
deposits. This material is permeable and contains water in large pores only after heavy rainfall
episodes that are rare and the corresponding water stock in this upper layer may be missed with
our monitoring. Measurements in core samples shown in Fig. (3) confirm the heterogeneous
geological formation, but they may be not fully representative for a much larger scale of the
MRS setup. Note that in this heterogeneous geology, one has to be careful when extending
observations made in one site to another located 1.5 km apart.
Another interpretation assumes horizontal flow channels. It is based on the observations made in the O-ZNS boreholes (Aldana et al. 2021). Borehole logs show a low permeable powdery limestone between 2 and 8 m deep, which does not allow rapid vertical transfer of rainwater. However, the MRS results show increasing water content below 10 m following the rain season (Figs. 8 and 12). We can explain this observation by horizontal flows above and below the powdery limestone formation that allow water circulation without runoff and without transfer through the powdery limestone. In this case, we assume hydraulic connections between the surface and channels above 2 m and channels below 10 m deep. These connections may exist somewhere off the monitoring site and thus are not detected with the MRS.

Figs. (11 and 12) demonstrate that the increase of the MRS estimated water column correlates with the rainfall and not with the rise of the groundwater level. Hence, a water accumulation in the VZ explains the rise of the MRS signal.

Results presented in Figs. (8 and 11) allow optimizing the MRS monitoring network and provide additional components for an unsaturated water flow modeling. Water accumulation in the VZ observed with MRS in 1999 (Fig. 8a) suggests at least four formations with different hydraulic properties. Fig. (8b) shows only three of them. The hydraulic system is not linear and its behavior depends on the quantity of the rainfall. For example, Fig. (11a) shows a large water accumulation in 1999 and a strong correlation between the $H_{MRS}$, the rainfall and the GWL. However, we do not observe these correlations in 2020 (Fig. 11b). The observation strategy requires denser sampling, especially during heavy rain events. An accurate water flow modeling requires at least a few years of MRS monitoring. This will allow getting data corresponding to large and small rainfalls, which is necessary when modeling non-linear systems.

## 6. Conclusions

Presenting this study, we show that time-lapse MRS measurements is an efficient observation tool for providing seasonal variations of the water quantity stored in a thick VZ (0–
20 m). MRS monitoring allows recording the volumetric water content averaged over the loop size (75×75 m²). We observe a good correlation of MRS results with the hydrological and meteorological observations.

For demonstration, we present results acquired in our experimental site near Villamblain and the O-ZNS observatory in France. We performed a one-year MRS monitoring in 1999 and another one in 2020. In these years, the hydro-meteorological conditions were different: the rainfall records show 837 mm in 1999 and 535 mm in 2020 and the ETP was 804 and 1044 mm respectively. When the rainfall was high (1999), the MRS showed a VZ composed of four layers with different water content. In 2020, the rainfall was smaller and the MRS showed only three layers. These observations suggest that the unsaturated water flow depends on the rainfall intensity.

Our data do not show fast water flow during and immediately after heavy rain events reported by Mazzilli et al., (2020). In order to better analyze large and small rainfall effects, and provide the necessary constraints to what appears as a non-linear flow modelling, we recommend a dense sampling rate during intensive rainfall events.

## Acknowledgements

This work was carried out in the framework of the collaborative research program of BRGM (https://www.brgm.fr/) and IRD-IGE (https://www.ige-grenoble.fr/). The authors acknowledge the support of the Institut des Sciences de la Terre d'Orléans, UMR 7327 Univ. d'Orléans (ISTO, https://www.isto-orleans.fr/) in the framework of the O-ZNS project, which is also a part of the PIVOTS project. We are thankful to the Région Centre-Val de Loire (ARD 2020 program and CPER 2015-2020) and the French Ministry of Higher Education and Research (CPER 2015-2020 and the BRGM public service) for their financial support. The fieldwork was co-funded by the European Union with the European Regional Development Fund and by the national research program Labex VOLTAIRE (ANR-10-LABX-100-01).
References

Abbas, M., Jardani, A., Soueid Ahmed, A., Revil, A., Brigaud, L., Bégassat, Ph., & Dupon, J. P., 2017. Redox potential distribution of an organic-rich contaminated site obtained by the inversion of self-potential data. Journal of Hydrology, 554 (2017) 111–127. https://doi.org/10.1016/j.jhydrol.2017.08.053.

Aldana, C., Isch, A., Bruand, A., Azaroual, M., Coquet, Y., 2021. Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer. Vadose Zone Journal, 20, e20127. https://doi.org/10.1002/vzj2.20127

Archie, G. E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of American Institute of Mining Metallurgical Engineers, 146. 54-62. https://doi.org/10.2118/942054-G.

Aster, R C, Borchers, B, Thurber, C H, 2011. Parameter Estimation and Inverse Problems.
Elsevier Academic Press. (376 pp., eBook ISBN: 9780123850492).

Binley, A., Hubbard, S. S., Huisman, J. A., Revil, A., Robinson, D. A., Singha, K., & Slater, L. D., 2015. The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales. Water Resour. Res., 51:3837–3866. https://doi.org/10.1002/2015WR017016.

Binley, A., Cassiani, G., & Deiana, R., 2010. Hydrogeophysics - Opportunities and Challenges.
Bollettino di Geofisica Teorica ed Applicata, 51, 4, 267-284.

Binley, A., Cassiani, G., Middleton, R., & Winship, P., 2002a. Vadose zone flow model parameterisation using cross-borehole radar and resistivity imaging. J. Hydrol, 267:147–159. https://doi.org/10.1016/S0022-1694(02)00146-4.

Blazevic, L., Bodet, L., Pasquet, S., Linde, N., Jougnot, D., 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. https://doi.org/10.3390/w12051230.
Brunet, P., Clément, R., & Bouvier, C., 2010. Monitoring soil water content and deficit using Electrical Resistivity Tomography (ERT) – A case study in the Cevennes area, France. Journal of Hydrology, 380, 146-153. https://doi.org/10.1016/j.jhydrol.2009.10.032.

Boucher, M., G. Favreau, M. Descloitres, J-M. Vouillamoz, S. Massuel, Y. Nazoumou, B. Cappelaere, A. Legchenko, 2009. Contribution of geophysical surveys to groundwater modelling of a porous aquifer in semiarid Niger: An overview. C. R. Geoscience, 341, 10-503 11, 800-809. doi:10.1016/j.crte.2009.07.008.

Bruand, A., Creusot, G., Quétin, P., Darthout, R., Raison, L., Courtemanche, P., Gaillard, H., 1997. Variabilité de la recharge de la nappe de Beauce: Rôle de l'irrigation et des caractéristiques du sol. Étude et Gestion des Sols Association Française pour l'Etude des Sols, 229-245.

Cassiani, G., Ferraris, S., Giustiniani, M., Deiana, R., & Strobbia, C., 2009b. Time-lapse surface-to-surface GPR measurements to monitor a controlled infiltration experiment. Boll. Geof. Teor. App., 50, 209-226.

Chalikakis, K., M. R Nielsen, and A. Legchenko, 2008. MRS applicability for a study of glacial sedimentary aquifers in Central Jutland, Denmark. J. Appl. Geophys., 66, 3-4, 176-530 187. https://doi.org/10.1016/j.jappgeo.2007.11.005.

Costabel, S., and U. Yaramanci, 2011. Relative hydraulic conductivity in the vadose zone from magnetic resonance sounding - Brooks-Corey parameterization of the capillary fringe. Geophysics, 76, 3, B1-B11. https://doi.org/10.1190/1.3552688.

Dafflon, B., Irving, J., & Barrash, W., 2011. Inversion of multiple intersecting high-resolution crosshole GPR profiles for hydrological characterization at the Boise hydrogeophysical research site. J. Appl. Geophys., 73:305-314. https://doi.org/10.1016/j.jappgeo.2011.02.001.

De Jong, S., Heijenk, R., Nijland, W., & Meijde, M., 2020. Monitoring Soil Moisture Dynamics Using Electrical Resistivity Tomography under Homogeneous Field Conditions. Sensors (Basel, Switzerland). https://doi.org/10.3390/s20185313.
Descloitres, M., L. Ruiz, M. Sekhar, A. Legchenko, J-J. Braun, M. S. Mohan Kumar, S. Subramanian, 2008. Characterization of seasonal local recharge using Electrical Resistivity Tomography and Magnetic Resonance Sounding. Hydrol. Processes, 22, 3, 384-394, 557. https://doi.org/10.1002/hyp.6608.

Eppelbaum, L.V., 2019. Geophysical Potential Fields: Geological and Environmental Applications. Elsevier Inc, (476 pp., Hardcover ISBN: 9780128196465, eBook ISBN: 9780128116869). Evett, S., 2003. Soil Water Measurement by Time Domain Reflectometry. In Encyclopedia of Water Science; Marcel Dekker, Inc.: New York, NY, USA; ISBN 0-8247-4241-9.

Falzone, S., and K. Keating, 2016. Algorithms for removing surface water signals from surface nuclear magnetic resonance infiltration surveys. Geophysics, 81, 4, WB97-568 WB107. https://doi.org/10.1190/geo2015-0386.1.

Flipo, N., Monteil, C., Poulin, M., de Fouquet, C., Krimissa, M., 2012. Hybrid fitting of a hydrosystem model: Long-term insight into the Beauce aquifer functioning (France). Water Resour. Res. 48, W05509. https://doi.org/10.1029/2011WR011092.

Gance, J., Malet, J. P., Supper, R., Sailhac, P., Ottowitz, D., & Jochum, B., 2016. Permanent electrical resistivity measurements for monitoring water circulation in clayey landslides. J. Appl. Geophys., 126, 28-115. https://doi.org/10.1016/j.jappgeo.2016.01.011.

Glover, P., 2016. Archie's law – a reappraisal. Solid Earth, 7. 1157-1169. DOI: 10.5194/se-7-1157-2016.

Herckenrath, D., E. Auken, L. Christiansen, A. A. Behroozmand, P. Bauer-Gottwein, 2012. Coupled hydrogeophysical inversion using time-lapse magnetic resonance sounding and time-lapse gravity data for hydraulic aquifer testing: Will it work in practice? Water Resour. Res., 48, W01539, doi:10.1029/2011WR010411.

Hertrich, M., M. Braun, T. Günther, A. G. Green, U. Yaramanci, 2007. Surface nuclear magnetic resonance tomography, IEEE Trans. Geosci. Remote Sens, 45, 11, 3752-3759, 591 doi:10.1109/TGRS.2007.903829.
Hubbard, S. S., Chen, J. S., Peterson, J., Majer, E. L., Williams, K., Swift, D. J., Mailloux, B., & Rubin, Y., 2001a. Hydrogeological characterization of the South Oyster Bacterial Transport Site using geophysical data. Water Resour. Res., 37(10), 2431–2456. https://doi.org/10.1029/2001WR000279.

Huisman, J. A., Hubbard, S. S., Redman, J. D., & Annan, P. A., 2003. Measuring soil water content with ground penetrating radar: A review. Vadose Zone J., 2(4), 476–491. https://doi.org/10.2113/2.4.476.

Isch, A., Montenach, D., Hammel, F., Ackerer, P., Coquet, Y., 2019. A Comparative Study of Water and Bromide Transport in a Bare Loam Soil Using Lysimeters and Field Plots. Water, 11, 1199. https://doi.org/10.3390/w11061199.

IUSS Working Group WRB, 2015. World Reference Base for Soil Resources 2014, update 2015. International soil classification system for naming soils and creating legends for soil maps. World Soil Resources Reports No. 106. FAO, Rome.

Jardani, A., Revil, A., & Dupont, J. P., 2012. Stochastic joint inversion of hydrogeophysical data for salt tracer test monitoring and hydraulic conductivity imaging. Advances in Water Resources. https://doi.org/10.1016/j.advwatres.2012.08.005.

Johnson, T., Versteeg, R., Thomle, J., Hammond, G., Chen, X., & Zachara, J., 2015. Four-dimensional electrical conductivity monitoring of stage-driven river water intrusion: Accounting for water table effects using a transient mesh boundary and conditional inversion constraints. Water Resour Res 51:6177–6196. https://doi.org/10.1002/2014WR016129.

Johnson, T. C., Versteeg, R. J., Ward, A., Day-Lewis, F. D., & Revil, A., 2010. Improved hydrogeophysical characterization and monitoring through parallel modeling and inversion of time-domain resistivity and induced-polarization data. Geophysics, 75, 27–41. https://doi.org/10.1190/1.3475513.

Kowalsky, M., Finsterle, S., Peterson, J., Hubbard, S., Rubin, Y., Majer, E., Ward, A., & Gee, G., 2005. Estimation of field-scale soil hydraulic and dielectric parameters through joint inversion of GPR and hydrological data. Water Resour. Res., 41, W11425, https://doi.org/10.1029/2005WR004237.
Legchenko, A. V., Shushakov, O. A., 1998. Inversion of surface NMR data. Geophysics, 63 (1):75-84. https://doi.org/10.1190/1.1444329.

Legchenko, A, Valla, P, 2002. A review of the basic principles for proton magnetic resonance sounding measurements. J. Appl. Geophys., 50:3-19. http://dx.doi.org/10.1016/S0926-9851(02)00127-1.

Legchenko, A., Baltassat, J-M., Beauce, A. and Jean Bernard, J., 2002. Nuclear magnetic resonance as a geophysical tool for hydrogeologists. J. Appl. Geophys., 50, 21-46. https://doi.org/10.1016/S0926-9851(02)00128-3.

Legchenko, A, Baltassat, J-M, Bobachev, A, Martin, C, Robin, H, Vouillamoz, J-M, 2004. Magnetic resonance sounding applied to aquifer characterization. J. Ground Water 42 (3):363-373. http://dx.doi.org/10.1111/j.1745-6584.2004.tb02684.x.

Legchenko, A., M. Ezersky, C. Camerlink, A. Al-Zoubi, K. Chalikakis, and J-F. Girard, 2008. Locating water-filled karst caverns and estimating their volume using magnetic resonance soundings. Geophysics, 73, 5, 51-61, doi:10.1190/1.2958007.

Legchenko, A, 2013. Magnetic resonance imaging for groundwater. Wiley-ISTE (158 pp., Online ISBN:9781118649459). https://doi.org/10.1002/9781118649459.

Legchenko, A, Comte J-C, Ofterdinger, U, Vouillamoz, J-M, Lawson, F M A, Walsh, J, 2017. Joint use of singular value decomposition and Monte-Carlo simulation for estimating uncertainty in surface NMR inversion. J. Appl. Geophys., 144:28-36. https://doi.org/10.1016/j.jappgeo.2017.06.010.

Legchenko, A, Baltassat, J-M, Duwig, C, Boucher, M, Girard, J-F, Soruco, A, Beauce, A, Mathieu, F, Legout, C, Descloitres, M, Flores Avilés, G B, 2020. Time-lapse magnetic resonance sounding measurements for numerical modeling of water flow in variably saturated media. J. Appl. Geophys., 175, Article 103984. https://doi.org/10.1016/j.jappgeo.2020.103984.
Legchenko, A., 2021. Surface NMR for Hydrogeology. IOP Publishing Ltd (168 pp., Online
ISBN: 978-0-7503-3155-5, Print ISBN: 978-0-7503-3153-1°.
https://doi.org/10.1088/978-0-7503-3155-5.

Lejars, C., Fusillier, J.L., Bouarfa, S., Coutant, C., Brunel, L., Rucheton, G., 2012. Limitation of agricultural groundwater uses in Beauce (France): what are the impacts on farms and on the food-processing sector ? *Irrigation and Drainage* 61, 54-64.
https://doi.org/10.1002/ird.1659.

Linde, N., Binley, A., Tryggvason, A., Pedersen, L. B., & Revil, A., 2006. Improved hydro-geophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data. Water Resour. Res., 42: W12404.
https://doi.org/10.1029/2006WR005131.

Looms, M. C., Binley, A., Jensen, K. H., Nielsen, L., & Hansen, T. M., 2008. Identifying unsaturated hydraulic parameters using an integrated data fusion approach on cross-borehole geophysical data. Vadose Zone J., 7:238-248. https://doi.org/10.2136/vadose zonej2007.0087.

Lubczynski, M., and Roy, J., 2003. Hydrogeological interpretation and potential of the new magnetic resonance sounding (MRS) method. J. Hydrol., 283, 19-40, doi:10.1016/S0022-666 1694(03)00170-7.

Lunt, I., Hubbard, S. S., & Rubin, Y., 2005. Soil moisture estimation using ground-penetrating radar reflection data. J. Hydrol., 307, 254-269.
https://doi.org/10.1016/j.jhydrol.2004.10.014.

Mallet, C., Isch, A., Laurent, G., Jodry, C., and Azaroual, M., 2022. Integrated static and dynamic geophysical and geomechanical data for transport property characterization. International Journal of Rock Mechanics and Mining Sciences. 22, 153, pp., 105050.
HAL Id: hal-03587248. https://hal.archives-ouvertes.fr/hal-03587248.
Mazzilli, N., Chalikakis, K., Simon, D., Carrière, S. D., Legchenko, A., 2020. Surface Nuclear Magnetic Resonance Monitoring Reveals Karst Unsaturated Zone Recharge Dynamics during a Rain Event. Water 2020, 12, 3183. https://doi.org/10.3390/w12113183.

Nielsen, M. R., Hagensen, T. F., Chalikakis, K., and Legchenko, A., 2011. Comparison of transmissivities from MRS and pumping tests in Denmark. Near Surf. Geophys., 9, 211-223, doi:10.3997/1873-0604.2010071.

Robinson, D., Binley, A., Crook, N., Day-Lewis, F., Ferré, T., Grauch, V.J.S, Knight, R., Knoll, M., Lakshmi, V., Miller, R., Nyquist, J., Pellerin, L., Singha, K., & Slater, L., 2008. Advancing process-based watershed hydrological research using near-surface geophysics: A vision for, and review of, electrical and magnetic geophysical methods. Hydrol. Processes, 22, 3604-3635. https://doi.org/10.1002/hyp.6963.

Robinson, D., Jones, S., Wraith, J., Or D., & Friedman, S., 2003. A review of advances in dielectric and electrical conductivity measurement in soils using time domain reflectometry. Vadose Zone J., 2, 444-475. https://doi.org/10.2113/2.4.444.

Ould Mohamed, S., Bruand, A., 1994. Morphology and origin of secondary calcite in soils from Beauce, France, in: Developments in Soil Science. Elsevier, pp. 27-36. https://doi.org/10.1016/S0166-2481(08)70395-7

Rubin, Y., & Hubbard, S. S., 2005. Hydrogeophysics, 523 pp., Springer, N. Y.

Schirov, M., A. Legchenko, and G. Creer, 1991. New direct non-invasive ground water detection technology for Australia. Exploration Geophysics, 22, 333-338.

Schnebelen, N, Ledoux, E, Bruand, A, Creuzot, G, 1999. Stratification hydrogéochimique et écoulements verticaux dans l'aquifère des calcaires de Beauce (France) : un système anthropisé à forte variabilit spatiale et temporelle. Comptes rendus hebd. séances acad. sci. Paris, 329:421-428. https://doi.org/10.1016/S1251-8050(00)80066-X.

Skierucha, W., Wilczek, A., Szypłowska, A., Sławiński, C., Lamorski, K., Skierucha, W., Wilczek, A., Szypłowska, A., Sławiński, C., & Lamorski, K., 2012. A TDR-based soil moisture monitoring system with simultaneous measurement of soil temperature and electrical conductivity. Sensors, 12:13545-13566. https://doi.org/10.3390/s121013545.
Steelman, C., & Endres, A., 2011. Comparison of Petrophysical Relationships for Soil Moisture Estimation using GPR Ground Waves. Vadose Zone Journal, 10. 270. https://doi.org/10.2136/vadose zonej2010.0040.

Topp, G. C., Davis, J. L., & Annan, A. P., 1980. Electromagnetic determination of soil water content: Measurements in coaxial transmission lines. Water Resour. Res., 16:574–582. https://doi.org/10.1029/WR016i003p00574.

Uhlemann, S., Chambers, J., Wilkinson, P., Maurer, H., Merritt, A., Meldrum, P., Kuras, O., Gunn, D., Smith, A. & Dijkstra, T., 2017. Four-dimensional imaging of moisture dynamics during landslide reactivation. Journal of Geophysical Research: Earth Surface, 122, 398-418. https://doi.org/10.1002/2016JF003983.

Valois, R., Vouillamoz, J.M., Lun, S., and Arnout, L., 2018. Mapping groundwater reserves in northwestern Cambodia with the combined use of data from lithologs and time-domain-electromagnetic and magnetic-resonance soundings. Hydrogeology Journal. 10.1007/s10040-018-1726-1.

Vereecken, H., Binley, A., Cassiani, G., Revil, A., & Titov, K., 2006. Applied Hydrogeophysics, 383 pp., Springer, Dordrecht, Netherlands.

Verhoef, A., Fernández-Gálvez, J., Diaz-Espejo, A., Main, B.E. & El-Bishti, M., 2006. The diurnal course of soil moisture as measured by various dielectric sensors: Effects of soil temperature and the implications for evaporation estimates. Journal of Hydrology, 321, 147-162. https://doi.org/10.1016/j.jhydrol.2005.07.039.

Vilhelmsen, T. N., A. A. Behroozmand, S. Christensen, and T. H. Nielsen, 2014. Joint inversion of aquifer test, MRS, and TEM data. Water Resour. Res., 50, 3956–3975, 758 doi:10.1002/2013WR014679.

Vouillamoz, J-M., J. Hoareau, M. Grammare, D. Caron, L. Nandagiri, A. Legchenko, 2012. Quantifying aquifer properties and freshwater resource in coastal barriers: a hydrogeophysical approach applied at Sasihithlu (Karnataka state, India). Hydrol. Earth Syst. Sci., 16, 4387–4400, doi:10.5194/hess-16-4387-2012.
Vouillamoz, J.M., Lawson, F.M.A., Yalo, N., Descloitres, M., 2014. The use of magnetic resonance sounding for quantifying specific yield and transmissivity in hard rock aquifers: The example of Benin. J. Appl. Geophys. 107, 16–24. http://dx.doi.org/10.1016/j.jappgeo.2014.05.012.

Walsh, D., Grunewald, E., Zhang, H., Ferre, P. and Hinnell, A., 2012. Recent advancements in NMR for characterizing the vadose zone. In: Recent advancements in NMR for characterizing the vadose zone. 5thInternational Workshop on Magnetic Resonance in the Subsurface, September 25 - 27, 2012, Hannover, Germany.

Wehrer, M., & Slater, L., 2015. Characterization of water content dynamics and tracer breakthrough by 3-D electrical resistivity tomography (ERT) under transient unsaturated conditions. Water Resources Research, 51. https://doi.org/10.1002/2014WR016131.

## Figure captions

Figure 1. Location of the MRS test site in Villamblain with respect to the Beauce aquifer and the Seine and Loire basins. Blue circles report main agriculture boreholes locations. Poiseaux and Bricy stations provide respectively the rainfall data and GWL. The map of the Beauce aquifer is modified from Flipo et al., (2012).

Figure 2. The annual rainfall (blue bars) and the groundwater level (GWL) (red line) collected between 1965 and 2021.

Figure 3. (a) The geological cross-section observed in the O-ZNS site; (b) physical properties measured in rock cores (granulometry, calcium carbonate percentage (red curve with dots), saturated volumetric water content (red curve) and density (black curve).

Figure 4. Photos of the limestone formation observed during construction of a 4 m diameter and 17 m deep pit at the O-ZNS Observatory. Figure 5. The annual rainfall recorded before the
monitoring time (1998 and 2019: 600 and 607 mm, respectively) and that corresponding to the
monitoring time (1999 and 2020: 837 and 535 mm, respectively).

Figure 6. An example of the MRS signals versus pulse moment observed during our study:
a) measurements under similar hydro-meteorological conditions; b) two soundings showing
the maximum and the minimum MRS signals. Solid lines show the theoretical signal computed
after the corresponding inverse models. For all these soundings, the noise level was similar
and it is shown only, for one sounding per graph, by the corresponding color in each graph.

Figure 7. a) The water content profiles corresponding to the minimum (black) and the maximum (red) MRS signals; b) corresponding relaxation times $T_2^*$; c) corresponding diagonal elements of the model resolution matrix; d) the resistivity profile provided by electrical measurements and the lithological log of the boreholes X0119 (03622X0119/FP3-25) and X0107 (03622X017/F) located near MRS station (see Fig. 1). GWL at both dates are within 30 cm of each other.

Figure 8. MRS inverse models of the water content versus time (color scale) computed down to 22 m and a daily rainfall corresponding to the monitoring time: a,c) 1999-2000 data set; b,d) 2020-2021 data set.

Figure 9. For each sounding, measured amplitude of the MRS signal is plotted versus normalized pulse moment of the corresponding data set and the inversion fit: a) 1999-2000 data set; b) 2020-2021 data set.

Figure 10. The water column estimated with MRS in 1999-2000 and in 2020-2021 plotted versus time. The zero-time corresponds to the 1st of January 1999 and the 1st of January 2020.

Figure 11. The equivalent water column ($H_{MRS}$), the monthly rainfall and the GWL versus time: a) 1999 – 2000 data sets; b) 2020 – 2021 data set. In these graphs, the vertical scale for plotting the MRS water column is set different for showing more details in figure 11b.
Figure 12. The water column estimated with MRS (red), the cumulative rainfall (black) and the depth of the GWL (blue) observed in Villamblain in 1999-2000.

* We present an innovative approach for non-invasive quantification of the water content in the unsaturated limestone

* We confirm the possibility of time-lapse MRS measurements of the water content

* In our test site, MRS shows seasonal and annual variations in the water content in the unsaturated zone

* We use existing hardware and software that renders our method easily available for users