The flow in an unconfined double-porosity aquifer with a sloping base is investigated and equations are developed for its description. In order to obtain a relatively simple description of the problem, similar assumptions as in Moutsopoulos (2021) have been adopted; the pressure in the unsaturated zone, especially in the fractures' network, is considered to be atmospheric, and the Dupuit-Forchheimer approximation is invoked, reducing the dimensionality by eliminating the vertical direction. The derived equations have been linearized and solved analytically for a problem involving interaction between an inclined aquifer and an adjacent surface water body similar to the one examined by Akylas and Koussis (2007). The analytical solution has been checked against results obtained with state-of-the-art numerical codes. The agreement between the two approaches is excellent. The solution tools were used to gain insight in the influence of the aquifer's base inclination on the flow quantities.

Both active and passive full waveform acoustic loggings (FWAL), complemented by a flow log, were conducted in a borehole of an experimental site located in the Cher region (France).. The acoustic tool used for the FWAL experiments is a flexible monopole tool holding a pair of piezoelectric receivers and a magnetostrictive transducer. The tool was modified to perform both active and passive FWAL. For passive acoustic logging, several runs were recorded to obtain a set of acoustic noise sections. As the noise is simultaneously recorded by two receivers of the tool, an interference noise section was elaborated by correlating or deconvolving the pair of signals and then summing these pairs of acoustic traces at each depth. This procedure, which can be interpreted as an interferometry analysis, points out the presence of low-frequency waves identified as Stoneley waves. The velocity and RMS amplitude of the Stoneley wave were computed at each depth. It is shown that: 1- the Stoneley wave velocity obtained in passive mode can be used to estimate the shear velocity of the formation, 2- the RMS amplitude and velocity variations of the Stoneley waves are strongly correlated with the variations of the flowmeter.

Mixed Finite Element (MFE) method is a robust numerical technique for solving elliptic and parabolic partial differential equations (PDEs). However, MFE can generate solutions with strong unphysical oscillations and/or large numerical diffusion for hyperbolic type PDEs. For its part, Discontinuous Galerkin (DG) finite element method is well adapted to solve hyperbolic systems and can accurately reproduce solutions involving sharp fronts. Therefore, the combination of DG and MFE is a good strategy for solving hyperbolic/parabolic problems such as advection – diffusion/dispersion equations. The classical formulation of the two methods is based on operator and time splitting allowing for separate solutions to advection with an explicit scheme and to dispersion with an implicit scheme. However, this kind of approach has the following drawbacks: (i) it lacks efficiency, as two systems with different unknowns are solved at each time step, (ii) it induces errors generated by the splitting, (iii) it can be CPU wise-expensive because of the CFL constraint, and (iv) it cannot be employed for steady-state transport simulations.To overcome these difficulties, we develop in this work a fully implicit edge/face centered DG-MFE formulation where the two methods share the same unknowns. In this formulation, the DG method is developed on lumping regions associated with the mesh edges/faces instead of mesh elements. Thus, the traces of concentration at mesh edges/faces, which are the Degrees Of Freedom (DOF) of the hybrid-MFE, are also part of the DOFs of the DG. The temporal discretization is based on the Crank–Nicolson method for both advection and dispersion. Numerical tests are performed to validate the new scheme by comparison against an analytical solution and to show its ability to handle steady-state transport simulations.The procedure is developed for 2D triangular meshes but can easily be extended to other 2D and 3D shape elements.

Two boreholes of an experimental site located in the Cher region (France) were investigated via Full Waveform Acoustic Logging (FWAL). The acoustic tool used for the FWAL experiments is a flexible monopole tool holding two pairs of piezoelectric receivers and a magnetostrictive transducer. The tool was modified to perform both active and passive FWAL. To our knowledge, this change is a novelty. For passive acoustic logging, several runs were recorded to obtain a set of acoustic noise sections from which noise Root Mean Squared (RMS) amplitude logs and spectral amplitude logs in different frequency bandwidths were computed. The acoustic logs resulting from passive acoustic monitoring were compared with P-wave acoustic velocity, core data, and a flowmeter log. It is shown that: (1) the distribution of noise frequencies in the 0–5 kHz is strongly correlated with the variations of the flowmeter, (2) the distribution of noise frequencies and noise RMS amplitude is correlated with the lithology (core description), and the P-wave velocity log. As the noise is simultaneously recorded by two receivers of the tool, an interference noise section was elaborated by correlating and summing pairs of acoustic traces at each depth. This procedure, which can be interpreted as an interferometry analysis, points out the presence of low-frequency waves identified as Stoneley waves. It is shown that the Stoneley wave velocity obtained in passive mode can be used to estimate the shear velocity of the formation.

Understanding subsurface flow, especially in partly karstified rock formations mainly housing water through a few preferential pathways, is still challenging. This point is the consequence of the poor accessibility of the subsurface and lack of accurate depictions of water bearing bodies and distributions. This notwithstanding, highly-resolved geophysical investigations bring new images of the subsurface. A 3-Dseismic surveywith shots and wavemonitoring at the surface is carried out over a subsurface karstified reservoir located at theHydrogeological Experimental Site (HES) of theUniversity of Poitiers (France). Processing the 3-D data, in association with wave velocity calibration from vertical seismic profiles ( VSP) recorded via geophones in wells, renders a 3-D velocity block. The velocity block is then converted into pseudo-porosity values revealing three high-porosity, presumably water-productive, layers, at depths of 35-40, 85-87, and 110-115 m. In addition, full wave acoustic logging (FWAL) can detect, close to wells, porous or open bodies that are too small for being captured by the spatial resolution of 3-D seismic images. A FWAL can also confirmor invalidate data fromVSP recorded via hydrophones. The block of pseudo-porosities is compared to a different representation of the subsurface in the form of hydraulic conductivity distributions (or hydraulic diffusion) obtained by slug tests or by inversion of transient interference testing between wells. The inverted hydraulic conductivity maps do not match up the distribution of porous bodies identified by seismic data. This poses the question of guiding conventional inversions on the basis of a prior guess as the subsurface structure obtained via geophysical investigations.

This study is geared towards evaluating the hydrological information that can be extracted from spring water temperature variations in shallow and thin aquifers of headwater catchments. A series of temperature variations (from 2013 to 2017) in four spring water sources at the Strengbach Critical Zone Observatory (CZO) was analyzed and interpreted by relying upon a coupled two-dimensional unsaturated-saturated flow and heat transfer model specifically designed for the study. Temperature variations of spring waters at the Strengbach catchment obey a seasonal non-phased and attenuated cyclicity compared with temperature inputs from the air or very shallow waters in soils. Simulation results show that, within the shallow subsurface horizons, heat transfer is mainly controlled by thermal conduction and not by fluid flow. The results also emphasize that the depth and temperature value of the thermal invariance zone are key parameters to surface temperature atten-uation patterns in the first meters beneath the surface. The average maximum porosity (or saturated water content) is another important parameter impacting heat transfers in the shallow subsurface. Stronger attenua-tions and delays of the thermal signal occur as the porosity increases. These different results indicate that temperature variations in the spring waters, especially the attenuation and phase shift of the signals compared with incoming air-surface water temperatures, do not bring information on hydrological transfers. In very shallow systems, characteristics, such as flow velocities and water residence times, cannot be inferred from temperature variations. Nevertheless, this study suggests that spring water temperature variations could become a way of identifying the usually poorly known saturated water content and its spatial variability at the catchment scale.

The paper underlines the contributions of Ghislain de Marsily (GdM) to the identification of aquifers heterogeneity using inverse methods mainly for modeling subsurface flow. Inverse methods require an objective function to express the goodness of fit of the chosen model, a parameterization to describe the spatial distribution of model parameters, and a minimization algorithm. The resulting inverse problem, which consists in seeking model parameters’ values that render model outputs close to the observations, is usually unstable. GdM developed seminal ideas for the two key inversion issues that are: to stabilize the inverse problem through regularization, and to parameterize it to reproduce the natural heterogeneity of the subsurface with a limited number of parameters. GdM conducted pioneering works that are the basis of current parameterization methods relying upon adaptive zonation and/or interpolation based on pilot points. We take here the opportunity to highlight the GdM’s contributions inspiring currently used techniques.

Summary Two boreholes of an experimental site located in the the Cher region (France) were investigated via full waveform acoustic logging (FWAL). The acoustic tool used for the FWAL experiments is a flexible monopole tool holding two pairs of piezoelectric receivers and a magnetostrictive transducer. The tool has been modified to perform both active and passive FWAL. To our knowledge, such a change is a novelty. For passive acoustic logging, several runs have been recorded to obtain a set of acoustic noise sections. Noise root mean squared (RMS) amplitude logs and spectral amplitude logs in different frequency bandwidths have been computed. The acoustic logs resulting from passive acoustic monitoring have been compared with P-wave acoustic velocity, core data and flowmetry. It is shown that: 1- the distribution of noise frequencies in the 0 -5 kHz are strongly correlated with the variations of the flowmeter, 2-the distribution of noise frequencies and noise RMS amplitude are correlated with the lithology (core description) and the P-wave velocity log.

We present an integrated petrological, petrophysical, and hydrogeological study of the critical zone (CZ) developed in the Hercynian granitic basement of the Strengbach watershed ( Vosges Massif, France) to characterize its deep architecture and water circulation levels. For this purpose, six boreholes (50-120mdepth), from which three are cored, and three piezometers (10-15mdepth) were drilled to define the vertical extension and lateral variability of the main CZ horizons. The Strengbach watershed is composed of a topsoil horizon of limited vertical extension (0.81.2 m), a mobile saprolite level, and an in-place fractured bedrock. The latter is subdivided into a few meters thick saprock horizon, defined by open sub-horizontal fractures and a deeper fractured bedrock horizon with steeply dipping fractures (>50 degrees). In the north-facing slope, the vertical extension of the mobile saprolite horizon increases from approximate to 1-2 m at the top of the slope to approximate to 9 m downstream, close to the valley bottom. In contrast, the south-facing and more easterly slope shows a mobile saprolite horizon with limited vertical extension (approximate to 2-3 m thick). Such a difference is associated with the existence of a knickpoint in the river bed, separating a downstream zone marked by currently active erosion from an upstream one, less prone to erosion, with preserved reliefs formed around 20 ka ago. The water circulation scheme within the Strengbach watershed involves two different systems: a subsurface circulation within the shallow aquifer, corresponding to the mobile saprolite horizon and the saprock, and a deeper circulation in the fractured bedrock. The water circulation in the fractured bedrock is controlled by fractures of regional orientations, linked to the Vosges massif and the Rhine Graben Tertiary tectonics, and partly to reactivated Hercynian fracture zones. The unaltered bedrock was not reached by any of the three cores. These results from the Strengbach CZ demonstrate theimportance of integrating geological history of the watershed, either the long-termgeological bedrock evolution or the Quaternary erosion patterns, to better understand and model the CZ hydrological functioning at the watershed scale.

A rare dataset of in-situ Be-10 from high-resolution depth profiles, soils, rock outcrops, and stream sediments is combined with geochemical analysis and modelling of regolith evolution to understand the variability of denudation rates in a mountain watershed (Strengbach critical zone observatory). High-resolution depth profiles are key to detect the presence of mobile regolith and to highlight how it affects the critical zone evolution. The modelling of regolith evolution and Be-10 concentrations along depth profiles allow us to estimate both the cosmic ray exposure age (19 kyr) and the mean denudation rate (22 mm kyr(-1)) of the regolith without any steady-state assumption on Be-10 concentrations. Comparison with maximum denudation rates inferred from topsoil samples collected from the surface of the depth profiles and calculated using the temporal steady-state assumption of Be-10 concentrations highlights an overestimation of denudation by a factor of two. Maximum spatially averaged denudation rates determined from stream sediment samples also likely overestimate denudation rates by a factor of two. These biases are significant for investigating the geomorphological evolution and we propose a method to correct denudation rates using the inherited Be-10 concentrations and the cosmic ray exposure age deduced from the high-resolution depth profiles. A key result is also that a steady state of Be-10 concentrations and a steady state of regolith thickness are two different equilibrium states that do not necessarily coincide. The comparison between locally corrected and spatially averaged denudation rates indicates that the watershed geomorphology is not in a topographic steady state but is modulated by regressive fluvial erosion. Nonetheless, our study demonstrates that even in a watershed where the steady-state assumption of Be-10 concentrations is not verified, the spatial variations of in-situ Be-10 concentrations in sediments still carry qualitatively relevant information on the geomorphological evolution of landscapes.

In this study, a novel experimental setup is proposed for which a column filled with glass beads and parallelepiped-shaped limestone beams is used to reconstruct a multiple fracture limestone media. The proposed setup produces asymmetric breakthrough curves (BTCs) that are consistent with the shape expected from the past field and lab-scale studies. Three transport experiments have been conducted under fast, medium, and slow flow velocity conditions. The research focuses on parameter and state estimation using Bayesian inference via Markov Chain Monte Carlo (MCMC) sampler, investigating the degree to which three models of transport through fractured media can reproduce the experimental results under the three flow conditions. The first transport model, named ADE, is based on the equivalent porous medium (EPM) approach and corresponds to the linear advection dispersion equation (ADE). The second model, named FOMIM (first-order mobile immobile), is based on the mobile/immobile approach and uses the dual porosity model with a linear first-order transfer between mobile and immobile regions. The third model, named NLMIM (non-linear mobile-immobile), uses a nonlinear transfer function between these two regions. The results of the three models show that almost all the unknown model input parameters can be well-estimated with narrow confidence intervals using the MCMC method. With respect to state estimation, the ADE model fails to reproduce correctly the tail of the BTCs observed under slow and medium flow conditions. The FOMIM model improves the tailing of the BTCs, but significant discrepancies remain between simulated and measured concentrations. The NLMIM model with velocity-dependent parameters is the only model that captures BTCs under all three conditions of slow, medium, and fast flow velocities.

Solute transport models based on the resolution of the 3-D Advection-Dispersion (AD) equation are frequently plagued by several numerical problems, which add to the high computational cost. The hydrological model NIHM (Normally Integrated Hydrological Model) was recently proposed as a tool simulating the hydrological responses of watersheds with shallow saturated aquifers by coupling surface flow and a low-dimensional subsurface system, including the vadose zone. In this paper, we couple the low-dimensional flow model NIHM with a transport module solving the AD to propose an approach that enables to reduce the dimensionality of both the flow and transport problems. In NIHM, the low-dimensionality in the subsurface compartment results from an integration along the local direction normal to the bedrock of the aquifer. NIHM was previously evaluated and applied to actual hydrosystems-without addressing mass transfers-and it showed its ability to capture various hydrological responses even from complex systems. However, the relevance of a low-dimensional approach to transport is not proven yet as the model reduction could also render approximated velocity fields inappropriate to mass transfer problems. The accuracy and computational efficiency of the proposed model have been thoroughly examined through various synthetic test cases under different hydrodynamic conditions to assess the influence of the reduction of dimensionality on solute transport simulations. The findings of this study demonstrate that the reduction of dimension remains suited to predicting solute transport behaviors in shallow subsurface systems while providing an important gain in computation time. This might be promising for various applications dealing with groundwater quality.

Qatar’s water resource has been largely overexploited, leading to the severe depletion of its aquifers and degradation of water quality due to saline intrusions. Qatar envisions employing regional aquifers to store water via forced injection of desalinated water and thus increase available from a few days to two months. A strategy for the implementation of forced injections is proposed based on a spatially distributed model of groundwater flow at the scale of the whole country. The model is based on calibration under steady-state flow conditions and for a two-dimensional single regional aquifer due to the lack of data. Injection scenarios include various mean injection rates at the scale of the whole system and are interpreted under the assumption that the additional storage should feed 2.7 M inhabitants for two months at a rate of 100 L/person/day. When this water supply stock is reached, the model is run to define the infiltration rate, which allows the stock to remain constant over time as a result of an even balance between infiltrations, withdrawals and also leaks or inlets through the boundary conditions of the system.

Understanding subsurface flow, especially in fractured rocks only housing water through a few preferential pathways, is still challenging. The point is mainly associated with the poor accessibility of the subsurface and the lack of accurate representations for both heterogeneity and spatial distribution of water bearing bodies. This notwithstanding, highly-resolved geophysical investigations bring new images of the subsurface. This is exemplified over a fractured limestone aquifer at the site scale (for example, that of the radius of influence of an extraction well). On an experimental site, situated in the Cher region (France), two boreholes have been drilled for field experiments. Full Waveform Acoustic Logging (FWAL) and seismic experiments were conducted. Hybrid seismic imaging, which consists in combining refraction and reflection seismic results, has been carried out. Based on a four-step procedure, the processing of refracted and reflected waves provided two sections. After assemblage, these sections produced in a first step an extended time reflectivity section starting from the surface and, in a second step, a section over depth after calibration with Vertical Seismic Profile (VSP) and acoustic data. However, even the Very High Resolution (VHR) seismic methods do not have a sufficient vertical resolution to describe accurately the geological formation. The acoustic sections were processed to separate the different wave fields, to extract the criss-cross events and to build a criss-cross index log. A log of fracturation index, based on both criss-cross index and P-wave velocity measurements, was computed to detect the presence of fractures. After calibration, and under the assumption that the slower the P-wave velocity, the higher the permeability – porosity, a 3D seismic block of reflection can inform on preferential areas where flow should occur. At the scale of an open wellbore, acoustic loggings that measure wave velocities over a short distance within the well also inform on open features crosscut by the well. Finally, flow log measurements confirm the occurrence of flowing horizons that were previously marked by both seismic and acoustic data. Seismic and acoustic data are therefore suited to image contrasted hydraulic properties over fractured subsurface systems usually poorly documented.

In the context of element migration in clay-rich media, self-diffusion coefficients of interlayer cations in swelling clay minerals obtained from molecular simulations are rarely used by macroscopic models predicting cation-exchange processes. Based on experiments and simulations, this study aims at (i) making a connection between molecular and sample scale processes to predict the dynamics of cation-exchange reactions between the interlayer space of millimetre disks of vermiculite and aqueous reservoirs, and (ii) assessing the role played by both self-diffusion and selectivity coefficients on this process. Time-resolved cation exchange experiments were performed using Ca-saturated vermiculite disks immersed in aqueous reservoirs with different NaCl or SrCl2 salinities. The results were reproduced via a finite-volume model constrained by (i) cation self-diffusion coefficients calculated by molecular dynamics simulations and (ii) interlayer selectivity coefficients drawn from "batch" cation-exchange isotherms. Results showed that considering the averaged values for both the cation-exchange selectivity coefficients and self-diffusion coefficients of the slowest interlayer cation led to good agreement between the experiments and simulations, validating the modelling strategy for the connection between the molecular and laboratory time scales. A sensitivity test regarding the influence of the two input parameters on the overall results was then performed. This study underlined a constrained upscaling strategy to better assess the role played by different intrinsic parameters of the clay/water systems (molecular self-diffusion coefficients in the interlayer space vs. selectivity coefficient) on the diffusion of cations during cation-exchange reaction in clay-rich media.

A 3D seismic survey was done on a near surface karstic reservoir located at the hydrogeological experimental site (HES) of the University of Poitiers (France). The processing of the 3D data led to obtaining a 3D velocity block in depth. The velocity block was converted in pseudo porosity. The resulting 3D seismic pseudo-porosity block reveals three high-porosity, presumably-water-productive layers, at depths of 30–40, 85–87 and 110–115 m. This paper shows how full wave acoustic logging (FWAL) can be used to validate the results obtained from the 3D seismic survey if the karstic body has a lateral extension over several seismic. If karstic bodies have a small extension, FWAL in open hole can be fruitfully used to: detect highly permeable bodies, thanks to measurements of acoustic energy and attenuation; detect the presence of karstic bodies characterized by a very strong attenuation of the different wave trains and a loss of continuity of acoustic sections; confirm the results obtained by vertical seismic profile (VSP) data. The field example also shows that acoustic attenuation of the total wavefield as well as conversion of downward-going P-wave in Stoneley waves observed on VSP data are strongly correlated with the presence of flow.

Abstract. Understanding the variability of the chemical composition of surface waters is a major issue for the scientific community. To date, the study of concentration–discharge relations has been intensively used to assess the spatiotemporal variability of the water chemistry at watershed scales. However, the lack of independent estimations of the water transit times within catchments limits the ability to model and predict the water chemistry with only geochemical approaches. In this study, a dimensionally reduced hydrological model coupling surface flow with subsurface flow (i.e., the Normally Integrated Hydrological Model, NIHM) has been used to constrain the distribution of the flow lines in a headwater catchment (Strengbach watershed, France). Then, hydrogeochemical simulations with the code KIRMAT (i.e., KInectic Reaction and MAss Transport) are performed to calculate the evolution of the water chemistry along the flow lines. Concentrations of dissolved silica (H4SiO4) and in basic cations (Na+, K+, Mg2+, and Ca2+) in the spring and piezometer waters are correctly reproduced with a simple integration along the flow lines. The seasonal variability of hydraulic conductivities along the slopes is a key process to understand the dynamics of flow lines and the changes of water transit times in the watershed. The covariation between flow velocities and active lengths of flow lines under changing hydrological conditions reduces the variability of water transit times and explains why transit times span much narrower variation ranges than the water discharges in the Strengbach catchment. These findings demonstrate that the general chemostatic behavior of the water chemistry is a direct consequence of the strong hydrological control of the water transit times within the catchment. Our results also show that a better knowledge of the relations between concentration and mean transit time (C–MTT relations) is an interesting new step to understand the diversity of C–Q shapes for chemical elements. The good match between the measured and modeled concentrations while respecting the water–rock interaction times provided by the hydrological simulations also shows that it is possible to capture the chemical composition of waters using simply determined reactive surfaces and experimental kinetic constants. The results of our simulations also strengthen the idea that the low surfaces calculated from the geometrical shapes of primary minerals are a good estimate of the reactive surfaces within the environment.

Magnetic Resonance Sounding (MRS) measurements are acquired at 16 stations in the Strengbach headwater catchment (Vosges Mountains – France). These data, rendering the vertical distribution of water contents in the subsurface, are used to show their potential in conditioning a hydrological model of the catchment, as described in the article “Magnetic resonance sounding measurements as posterior information to condition hydrological model parameters: Application to a hard-rock headwater catchment” – Journal of Hydrology (2020). Acquisition protocols follow a free induction decay scheme. Data are filtered by applying a band-pass filter at the Larmor frequency. A filter removing the 50 Hz noise is also applied with the exception of data at a Larmor frequency close to the 50 Hz harmonic. The signal envelopes are then fitted by a decaying exponential function over time to estimate the median characteristic relaxation time of each MRS sounding.

In headwater catchment, the calibration of hydrological models is complex due to the scarcity of data in mountainous areas. Here, an innovative methodology is developed to condition hydrological model parameters by using magnetic resonance sounding (MRS) measurements in combination with stream flow rate data. MRS has the specificity in the various geophysical imaging techniques of being mainly sensitive to the vertical distribution of water content among the subsurface. In a way very similar to hydraulic head observations, these local distributions of water content may serve as information in a hydrological model to pattern subsurface flow by seeking model parameters. Simulations are run with different sets of parameters of a hydrological model. Each simulation provides as an output a 4-D map (3-D spatial plus time) of the vertical water content distributions over the whole catchment and their fluctuations over time. This output is then used to simulate the MRS signal that would be produced by the estimated water content. The simulated MRS signal is compared to measured MRS data to determine which hydrological simulations (which model parameters) are close to observations. The approach is applied on a hard-rock headwater catchment housing a very shallow and thin aquifer where an MRS survey covers the whole studied site. Hydraulic parameters of an integrated hydrological model of the catchment are spatially distributed by zones with uniform values, the prior delineation of the zones being guided by pedological studies. As MRS measurements supply local but spatially distributed information, the method conditions the various zones on their parameter values in a much better way than the classical (in headwater catchments) measure of the stream flow rate at the outlet of the system. Finally, hydrological simulation and time-dependent MRS forward calculations can help identifying possible locations for MRS stations to monitor the transient behavior of the hydrological state of the catchment.