Arnaud Isch
email: arnaud.isch@gmail.com

Yves Coquet

Bouamama Abbar

Carlos Aldana

Mohamad Abbas

Ary Bruand

Mohamed Azaroual

A comprehensive experimental and numerical analysis of water flow and travel time in a highly heterogeneous vadose zone

Keywords: rainfall, TT, water travel time, WTL, water table level, VZ, vadose zone Limestone aquifer, Hydraulic Properties, Unimodal and bimodal approaches, HYDRUS-1D, Vertical and lateral heterogeneity, Meteorological conditions

 Simulation of water flow and TT in the VZ using unimodal and bimodal approaches  Parameterization of HYDRUS-1D using hydraulic properties laboratory measurements  Experimental hydraulic properties were more accurately described with the DBM model  Results obtained with unimodal and bimodal models presented strong differences  Large influence of VZ vertical heterogeneity and meteorological conditions

Introduction

An increasing number of Earth's Critical Zone (CZ) networks and observatories have emerged in recent years (Gaillardet et al., 2018;[START_REF] Guo | Critical Zone Research and Observatories: Current Status and Future Perspectives[END_REF][START_REF] Zacharias | A Network of Terrestrial Environmental Observatories in Germany[END_REF], stimulated by societal challenges like the prediction of natural hazards, carbon storage or water resources management. These observatories are dedicated to deciphering the complex hydrological and biogeochemical processes which occur at various temporal and spatial scales in the CZ [START_REF] Bogena | The TERENO-Rur Hydrological Observatory: A Multiscale Multi-Compartment Research Platform for the Advancement of Hydrological Science[END_REF][START_REF] Brantley | Susquehanna Shale Hills Critical Zone Observatory: Shale Hills in the Context of Shaver's Creek Watershed[END_REF][START_REF] Galle | AMMA-CATCH, a Critical Zone Observatory in West Africa Monitoring a Region in[END_REF][START_REF] Jensen | HOBE: The Danish Hydrological Observatory[END_REF][START_REF] Li | The Effect of Lithology and Agriculture at the Susquehanna Shale Hills Critical Zone Observatory[END_REF][START_REF] Molénat | OMERE: A Long-Term Observatory of Soil and Water Resources, in Interaction with Agricultural and Land Management in Mediterranean Hilly Catchments[END_REF][START_REF] O'geen | Southern Sierra Critical Zone Observatory and Kings River Experimental Watersheds: A Synthesis of Measurements, New Insights, and Future Directions[END_REF][START_REF] Pisinaras | Hydrologic and Geochemical Research at Pinios Hydrologic Observatory: Initial Results[END_REF][START_REF] Seyfried | Reynolds Creek Experimental Watershed and Critical Zone Observatory[END_REF]. The Vadose Zone (VZ), which extends from the soil surface down to the aquifer, is considered as a key component of the CZ and is of major importance in the hydrological and biogeochemical cycles. Consequently, an increasing body of work aims at improving the knowledge of the VZ functioning, which includes numerous complex and coupled physical, geochemical, and microbial processes [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Arora | Understanding and Predicting Vadose Zone Processes[END_REF][START_REF] Vereecken | Modeling Soil Processes: Review, Key Challenges, and New Perspectives[END_REF]. This constitutes a fundamental milestone for a proper assessment of groundwater recharge and in managing contamination issues for the preservation of aquifers. This is especially true since the increasing use of groundwater resources worldwide and the degradation of the quality of the water resources caused by intensive agriculture have led to the enforcement of several stringent regulatory measures in several countries geared towards achieving a good chemical status for all water bodies in years ahead [START_REF] Vero | A framework for determining unsaturated zone water quality time lags at catchment scale[END_REF].

In this context, the estimation of the water travel time through the VZ is crucial as it constitutes an important component of the total time lag defining the inherent hydrologic response delay of water bodies to a change in land management practices. This response time is often conceptualized as consisting in both a vertical (within the VZ) and a lateral (within the saturated zone) component [START_REF] Sousa | A simple method to assess unsaturated zone time lag in the travel time from ground surface to receptor[END_REF]. The knowledge of the water travel time through the VZ is considered as a prerequisite for quantifying the recharge rate and the time of migration of contaminants from the soil surface down to the aquifer [START_REF] Mattern | Estimating travel time of recharge water through a deep vadose zone using a transfer function model[END_REF]. However, the water travel time in the VZ is known to be highly variable in space and time. Its quantification is recognized as highly complex and influenced by several factors including the size and geometry of the VZ, the properties of the VZ materials, the physical, chemical and biological processes governing water flow and solute transport within the VZ and the pattern of the meteorological inputs [START_REF] Jeong | A method of estimating sequential average unsaturated zone travel times from precipitation and water table level time series data[END_REF][START_REF] Leij | Solute transport[END_REF][START_REF] Sousa | A simple method to assess unsaturated zone time lag in the travel time from ground surface to receptor[END_REF][START_REF] Sprenger | Travel times in the vadose zone: Variability in space and time[END_REF][START_REF] Vero | Consequences of varied soil hydraulic and meteorological complexity on unsaturated zone time lag estimates[END_REF]. Numerous experimental methods are currently used to estimate the water travel time through the VZ, such as lysimeters (mainly focused on the soil) or field experiments using radioactive or stable isotopes, chloride, nitrate, bromide or fluorobenzoate as tracers [START_REF] Asadollahi | Transport and Water Age Dynamics in Soils: A Comparative Study of Spatially Integrated and Spatially Explicit Models[END_REF][START_REF] Chen | Water, nitrate and atrazine transfer through the unsaturated zone of the Chalk aquifer in northern France[END_REF][START_REF] Ju | Characterization of flow dynamics around highly-utilized agricultural wells in a fractured-rock aquifer: Assessment of uncertainties lying on groundwater age-dating[END_REF][START_REF] Sprenger | Travel times in the vadose zone: Variability in space and time[END_REF][START_REF] Turkeltaub | Recharge and Nitrate Transport Through the Deep Vadose Zone of the Loess Plateau: A Regional-Scale Model Investigation[END_REF]. However, these experimental approaches are known to be fairly expensive, labor intensive, time consuming and may also have limited predictive capacity. As a consequence, numerical models are being increasingly used for the study of water flow, solute transport and travel time assessments within the VZ [START_REF] Bouraoui | Modelling mitigation options to reduce diffuse nitrogen water pollution from agriculture[END_REF][START_REF] Konikow | The Secret to Successful Solute-Transport Modeling[END_REF][START_REF] Schoups | Model complexity control for hydrologic prediction[END_REF][START_REF] Szymkiewicz | Numerical Analysis of Recharge Rates and Contaminant Travel Time in Layered Unsaturated Soils[END_REF][START_REF] Vero | Consequences of varied soil hydraulic and meteorological complexity on unsaturated zone time lag estimates[END_REF]. These models need an estimation of the VZ hydrodynamic parameters and may partly rely on pedotransfer functions [START_REF] Patil | Pedotransfer Functions for Estimating Soil Hydraulic Properties: A Review[END_REF][START_REF] Schaap | Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions[END_REF][START_REF] Szymkiewicz | Estimation of Conservative Contaminant Travel Time through Vadose Zone Based on Transient and Steady Flow Approaches[END_REF][START_REF] Vereecken | Using Pedotransfer Functions to Estimate the van Genuchten-Mualem Soil Hydraulic Properties: A Review[END_REF][START_REF] Vero | A framework for determining unsaturated zone water quality time lags at catchment scale[END_REF][START_REF] Zhang | Estimation of saturated hydraulic conductivity with pedotransfer functions: A review[END_REF]. Some authors used inverse modeling approaches based on transient measurements of water content, matric head and tracer content in the VZ [START_REF] Haws | Single-porosity and dual porosity modeling of water flow and solute transport in subsurface-drained fields using effective field-scale parameters[END_REF][START_REF] Jacques | Calibration of Richards' and convection-dispersion equations to field-scale water flow and solute transport under rainfall conditions[END_REF][START_REF] Sprenger | Travel times in the vadose zone: Variability in space and time[END_REF][START_REF] Sprenger | Estimating flow and transport parameters in the unsaturated zone with pore water stable isotopes[END_REF]Varvaris et al., 2021aVarvaris et al., , 2021b)). However, hydraulic data acquisition is mostly focused on the soil compartment and the numerical simulations run for the estimation of the travel time are often limited to the use of a unimodal approach, relying notably on the Mualem-van Genuchten model to describe the VZ hydraulic properties [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF].

The main objective of this work was to provide a new methodological approach for the study of the water flow and the estimation of the water travel time at the scale of a highly heterogeneous VZ. To this end, the hydraulic properties (water retention and hydraulic conductivity curves) of fifteen samples taken from three cored boreholes and representative of the lithologies (soil, powdery limestone, calcareous sand, limestone rock) encountered within the VZ of a vulnerable limestone aquifer (Beauce, France) were first determined in the laboratory. These experimental data were then used for the parameterization of analytical models describing the measured hydraulic properties using single (unimodal) and dual porosity (bimodal) approaches. The travel time within the VZ has been estimated by means of a virtual bromide tracing experiment performed over a period of 55 years using the HYDRUS-1D software (Šimůnek et al., 2016). The appearance of the initial breakthrough, the peak concentration and the total exit of the tracer at the maximum water table level of the aquifer were considered as indicators of the water travel time. The impact of the choice of the analytical model used to describe the experimental hydraulic properties, the vertical and lateral heterogeneity of the VZ and the meteorological conditions on the simulation of the water flow and travel time through the VZ has been explored.

Materials & Methods

The Beauce limestone aquifer

The multi-layered Beauce limestone aquifer is located in the center of France and extends over 9700 km 2 (de Frutos Cachorro et al., 2017). It is mainly composed of limestones from upper Oligocene to lower Miocene but also consists in clay, marl, powdery limestone, calcareous sand or chalk facies depending on the location. This mostly unconfined aquifer is spread between the Seine (northeast) and the Loire (southwest) rivers and constitutes one of the largest groundwater reservoirs in France with an average water stock of 20 billion m 3 and high interannual variations [START_REF] Coz | Gestion durable d'une ressource en eaux souterraines Cas de la nappe de Beauce[END_REF]. The Beauce represents the main cereal-producing region in Europe [START_REF] Graveline | Combining flexible regulatory and economic instruments for agriculture water demand control under climate change in Beauce[END_REF]. Its land use consists essentially in agriculture (74%), 50% of which is irrigated. Indeed, the region is one of the driest in France with an annual average precipitation of 600 mm [START_REF] Lejars | Limitation of agricultural groundwater uses in Beauce (France): what are the impacts on farms and on the food-processing sector[END_REF]. This results in an intense irrigated farming of between 120 000 and 240 000 ha that withdraws from the groundwater resource between 150 and 450 million m 3 of water per year [START_REF] Graveline | Combining flexible regulatory and economic instruments for agriculture water demand control under climate change in Beauce[END_REF]. The sensitivity of the water resource to natural recharge and the observation of an overall lowering of the water table level since the early 1990s have led to the implementation of several local water policy measures to control water withdrawals in connection with the French Water Law (Légifrance, 1992), which is now the translation of the European Water Framework Directive (European Commission, 2000). Over the past decades, the Beauce intensive agriculture has also impacted the groundwater quality with nitrates and pesticides (notably atrazine and bentazon) whose concentrations still exceed regulatory limits (DDT, 2016) despite government policies and regulatory measures that have been put into effect since the early 1990s (European Commission, 2009, 1991).

Location, sampling and description of the VZ materials

The experimental site is located at Villamblain, 30 km northwest of Orléans (France) (DMS coordinates: X = 48°1'5.131"; Y=1°34'55.333"). Three cored boreholes (B1, B2 and B3) separated by a maximum distance of 10.8 m from each other (Figure 1 As stated by [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF], compaction caused by hydraulic percussion led to unreliable hydraulic properties data for the soil samples. Consequently, additional soil samples were collected manually from the ground surface after digging a pit at a few meters from the three boreholes (Figure 1) and down to a depth of 1 m. The VZ lithologies were described based on visual examination of the undisturbed core samples. It revealed a high vertical (within a single borehole) and lateral (between the three boreholes) lithological heterogeneity [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Isch | Material Characteristics, Hydraulic Properties, and Water Travel Time through the Heterogeneous Vadose Zone of a Cenozoic Limestone Aquifer[END_REF]. A soil with a thickness between 0.9 and 1.8 m was observed at the top of the VZ. This soil is typical of the Beauce region and is referred to as a loamy clay Calcisol [START_REF] Duval | Notice explicative de la carte pédologique de Villamblain (Beauce) au 1/10 000e, Contrat de recherche site expérimental de Villamblain, rapport période 1991-1992[END_REF][START_REF] Michot | Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography[END_REF] according to the French soil reference system [START_REF] Baize | Référentiel pédologique 2008[END_REF], or a Hypereutric Cambisol with a silt loam texture according to the World Reference Base for Soils (IUSS Working Group WRB, 2015). Three main soil layers were observed during the visual description of the core samples and the soil pit (0.0-0.3 m; 0.3-0.6 m; > 0.6 m) [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Isch | Material Characteristics, Hydraulic Properties, and Water Travel Time through the Heterogeneous Vadose Zone of a Cenozoic Limestone Aquifer[END_REF]. The thickness of the soil layers was in good agreement with that described in other studies carried out near the experimental site [START_REF] Michot | Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography[END_REF][START_REF] Ould Mohamed | Estimating Long-Term Drainage at a Regional Scale Using a Deterministic Model[END_REF]. The soil is developed on a highly heterogeneous Miocene lacustrine fragmented powdery limestone facies (referred to as "powdery limestone" in this work) that has been cryoturbated in its upper part during the Quaternary [START_REF] Michot | Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography[END_REF][START_REF] Ould Mohamed | Morphology and origin of secondary calcite in soils from Beauce, France[END_REF] and which also contains calcareous sand interbeds [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF]. The thickness of this layer is from 5.2 to 6.8 m. Finally, the last main stratigraphic facies is a 12.2 to 13.4 m thick massive, fractured or weathered hard rock (Pithiviers limestone) [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Schnebelen | Stratification hydrogéochimique et écoulements verticaux dans l'aquifère des calcaires de Beauce (France): Un système anthropisé à forte variabilité spatiale et temporelle[END_REF].

Climate and water table level data

Meteorological data were collected from the Bricy weather station located about 20 km east of the study site. From 1966 to 2020, the mean annual rainfall (R) and reference Penman-Monteith potential evapotranspiration (ETP) [START_REF] Monteith | Evaporation and environment[END_REF] (1977-1984 and 2012-2014). The lowering of the WTL observed from 1989 to the 19/08/1992, which led to a WTL MIN of -22.21 m at that date, was linked to low annual R and R-ETP balance (below -300 mm) over the period. The same trend was observed during other periods corresponding to a lowering WTL (1970-1974, 2003-2011 and 2017-2020).

Measurements of the hydraulic properties of the vadose zone samples

Following the drilling phase, twelve undisturbed samples representative of the VZ materials were selected from the three cored boreholes, based on the observations made during the lithological description. Three types of lithology were sampled throughout the VZ profile (1-20 m deep) including soft sediments with four powdery limestone samples (P A , P B , P C and P D )

and two calcareous sand samples (I A and I B ), and hard materials with six limestone rock samples (R A , R B , R C , R D , R E and R F ). Additionally, three undisturbed soil samples (S A , S B and S C ) were collected from a soil pit. Their hydraulic properties were determined according to the same procedure as for the other VZ samples.

The undisturbed samples were re-cored in the laboratory for adaptation to the experimental device, except sample I B (calcareous sand deposit) which was repacked to the observed field bulk density. The unsaturated hydraulic properties of the samples, i.e., the water retention curve (h) and the hydraulic conductivity curve K(h), were measured within the pF (log |h|, h in cm) value range 0.5-3.0 by means of a triaxial system used by applying the multistep outflow method [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF]. Calculation of the hydraulic conductivity (K) from outflow data was based on the method given by [START_REF] Gardner | Calculation of Capillary Conductivity from Pressure Plate Outflow Data[END_REF] assuming that the matric head in the samples varies linearly with depth. The measurement of the water retention curve was extended to the pF range 4.0-6.0 using the WP4C Dewpoint Potentiometer (METER Group ® ). A comprehensive description of the measurement system, experimental protocol and calculation methods can be found in [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF].

The experimental water retention and hydraulic conductivity data were then fitted to the (h)

and K(h) curves with three analytical models using the RETC software [START_REF] Van Genuchten | The RETC code for quantifying hydraulic functions of unsaturated soils[END_REF]. The fitting procedure is described in § 2.5.4.

Simulation of water flow and solute transport in the vadose zone

Theoretical background and general assumptions

The simulation of water flow and solute transport in the VZ was made using the HYDRUS-1D software (Šimůnek et al., 2016). The one-dimensional vertical water flow in the VZ was described by the Richards equation [START_REF] Richards | Capillary conduction of liquids through porous mediums[END_REF] (1):

∂𝜃 ∂𝑡 = ∂ ∂𝑧 [ 𝐾 ( ∂ℎ ∂𝑧 + 1 )] -𝑆 (1)
with  the volumetric water content (cm 3 /cm 3 ), t the time (d), z the coordinate along the vertical axis pointing positively upwards (cm), h the matric head (cm), K the hydraulic conductivity (cm/d) and S a sink term (d -1 ) which accounts for root water uptake.

We used three analytical models for the description of the hydraulic properties, i.e., the water retention characteristic (h) and the hydraulic conductivity characteristic K(h), of the VZ samples.

The first one is the Mualem-van Genuchten (MVG) model which uses the van Genuchten's expression [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] to describe the water retention curve (2) and the statistical pore connection model established by [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] to predict the hydraulic conductivity from the water retention curve (3):

𝜃(ℎ) = { 𝜃 𝑟 + 𝜃 𝑠 -𝜃 𝑟 [1 + |𝛼ℎ| 𝑛 ] 𝑚 ℎ < 0 𝜃 𝑠 ℎ ≥ 0 with 𝑚 = 1 - 1 𝑛 𝑛 > 1 (2)
with et respectively the residual and saturated volumetric water content (cm 3 /cm 3 ),  an 𝜃 𝑟 𝜃 𝑠 empirical parameter related to the matric head at the inflection point of the retention curve (cm - 1 ) and n a pore size distribution parameter (-) which determines the slope of the curve at the inflection point.

𝐾(ℎ) = 𝐾 𝑠 𝑆 𝑙 𝑒 [ 1 -(1 -𝑆 𝑒 1 𝑚 ) 𝑚 ] 2 with 𝑆 𝑒 = 𝜃 -𝜃 𝑟 𝜃 𝑠 -𝜃 𝑟 (3)
with the saturated hydraulic conductivity (cm/d), the effective saturation (-) and l a pore 𝐾 𝑠 𝑆 𝑒 connectivity parameter (-). The latter was fixed at the value of 0.5 [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF].

The second model (BBC) used to describe the hydraulic properties implements the equations of Brooks andCorey (1966, 1964) to describe the water retention curve (4) and the pore-size distribution model of [START_REF] Burdine | Relative Permeability Calculations From Pore Size Distribution Data[END_REF] to predict the hydraulic conductivity (5):

𝑆 𝑒 = { |𝛼ℎ| -𝑛 ℎ < -1/𝛼 1 ℎ ≥ -1/𝛼 (4) 𝐾(ℎ) = 𝐾 𝑠 𝑆 2 𝑛 + 𝑙 + 2 𝑒 (5)
The third model is the Durner bimodal expression (DBM) [START_REF] Durner | Hydraulic conductivity estimation for soils with heterogeneous pore structure[END_REF] which extends the unimodal MVG model by considering the porous medium composed of two overlapping regions and fitting the bimodal behavior of the hydraulic properties by a MVG-type function for each of the two regions [START_REF] Köhne | Estimating the hydraulic functions of dual permeability models from bulk soil data[END_REF][START_REF] Dimitrov | Soil Hydraulic Parameters and Surface Soil Moisture of a Tilled Bare Soil Plot Inversely Derived from L-Band Brightness Temperatures[END_REF] (6) (7):

𝑆 𝑒 = 𝑤 1 [1 + (𝛼 1 ℎ) 𝑛 1 ] -𝑚 1 + 𝑤 2 [1 + (𝛼 2 ℎ) 𝑛 2 ] -𝑚 2 (6) K(𝑆 𝑒 ) = (𝑤 1 𝑆 𝑒1 + 𝑤 2 𝑆 𝑒2) 𝑙 (𝑤 1 𝛼 1[ 1 -(1 -𝑆 𝑒 𝑙/𝑚 1 1 ) 𝑚 1 ] + 𝑤 2 𝛼 2[ 1 -(1 -𝑆 𝑒 𝑙/𝑚 2 2 ) 𝑚 2 ]) 2 (𝑤 1 𝛼 1 + 𝑤 2 𝛼 2 ) 2 (7)
with w i the weighing factors of the two regions (bimodal pore system) which display separate hydraulic functions and empirical parameters  i , n i , m i and l (subscript i=1 and i=2).

The root water uptake model was that of [START_REF] Feddes | Simulation of Field Water Use and Crop Yield[END_REF] with no solute uptake. HYDRUS-1D implements the same linear interpolation scheme for the water stress response function as the one used in the SWATRE code [START_REF] Wesseling | Introduction of the occurrence of high groundwater levels and surface water storage in computer program SWATRE[END_REF]. The water stress response function parameters of grass were used [START_REF] Taylor | Physical edaphology: the physics of irrigated and nonirrigated soils[END_REF][START_REF] Wesseling | SWATRE: Instructions for input[END_REF].

The root depth of grass was fixed to 30 cm [START_REF] Schenk | The global biogeography of roots[END_REF].

The one-dimensional vertical solute transport in the VZ was described by the convectiondispersion equation (CDE) (8):

∂(𝜃𝐶) ∂𝑡 = ∂ ∂𝑧 ( 𝜃𝐷 ∂𝐶 ∂𝑧 ) -∂(𝑞𝐶) ∂𝑧 (8)
with C, the solute concentration in liquid phase (g/cm 3 ), q, the water flux density (cm/d) and D, the hydrodynamic dispersion coefficient (cm²/d) given by Bear (1972) (9):

𝐷 = 𝜆𝑉 + 𝐷 0 𝜏 𝜃 (9)
with  the material dispersivity (cm), V, the average water velocity in the pores of the material (cm/d), D 0 , the solute molecular diffusion coefficient in pure water (1.584 cm²/d for bromide [START_REF] Lide | CRC Handbook of Chemistry and Physics[END_REF]) and , the tortuosity factor in liquid phase (-) given by [START_REF] Millington | Permeability of porous solids[END_REF].

Representation of the vadose zone profiles

The VZ profiles were reconstituted in HYDRUS-1D for each of the three boreholes (B1, B2 and B3, Figure 1). The lithological heterogeneities along each VZ profile were reproduced based on the visual descriptions of the undisturbed cored samples [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF]. A 23 m deep profile composed of fourteen (B1) and fifteen (B2 and B3) different materials was created for each borehole (Table 1). Each material corresponds to a sample whose hydraulic properties have been determined in the laboratory. The last material (R F ) has been extended from the maximum drilling depth (20 m) to the maximum profile depth (23 m) in order to allow the WTL to move below 20 m deep.

The three VZ profiles were vertically discretized using (i) a 1-cm mesh size, from 0 to 1 m deep and from 10 to 0 cm above and from 0 to 10 cm below the limit between two successive materials; (ii) a 2-cm mesh size, from 50 to 10 cm above and from 10 to 50 cm below the limit between two successive materials; (iii) a 5-cm mesh size everywhere else. Ten observations nodes were implemented in the B2 profile. Each observation node was placed at the bottom of a given material and was associated with a specific color (Table 1). For the comparisons made between the three boreholes, ten observations nodes were implemented within the three VZ profiles with a specific color assigned to a given borehole (Table 1).

Initial and boundary conditions

As the water table level at the time of the start of the simulation was -19.57 m, the initial matric head profile was defined with h i = -100 cm for -18.57 < z < 0.00 m and h i varying linearly from -100 to +343 cm from z = -18.57 to z = -23.00 m. At the soil surface, a water flux was imposed as upper boundary condition using daily ETP and rainfall data (atmospheric boundary condition with surface layer (2 cm)). The daily variations of the water table level were used as lower boundary condition (variable pressure head).

For the solute transport modelling, a time-dependent concentration of bromide was applied at the soil surface boundary. More details about the concentration of the bromide in the incoming water are given in § 2.5.5. A Cauchy type condition was imposed at the soil surface with dispersive flux considered negligible and a Neumann type condition was imposed at the lower limit of the VZ profile.

Table 1: Representation of the VZ profile for each borehole (B1, B2 and B3): material number, sample identifier, depth interval and depth of the observation nodes.

For the study made on the B2 profile, ten observations nodes were implemented and each of them was placed at the bottom of a given material and associated with a specific color. For the comparisons made between the three boreholes (B1, B2 and B3), ten observations nodes were implemented within the three VZ profiles with a specific color assigned to a given borehole.

B1 B2 B3 Material Sample Interval [m] Obs. Node B1-B2-B3 [m] Material Sample Interval [m] Obs. Node B2 [m] Obs. Node B1-B2-B3 [m] Material Sample Interval [m] Obs. Node B1-B2-B3 [m] 1 S A 0.

Description of the fitting procedure

The RETC software [START_REF] Van Genuchten | The RETC code for quantifying hydraulic functions of unsaturated soils[END_REF] was used to fit the (h) and K(h) curves (Eq. 2 to 7) to the measured water retention and hydraulic conductivity data. For the three boreholes, the hydrodynamic parameters for the three analytical models (MVG, BBC and DBM) were obtained according to the following steps. obtained by using the Rosetta software [START_REF] Schaap | Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions[END_REF] based on the particle size distribution and bulk density measured for each sample [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF]. For the limestone rock samples (R A , R B , R C , R D , R E and R F ), initial values of , , n were obtained using the 𝜃 𝑟 results of a study performed near the study site [START_REF] Amraoui | Projet ESPOL : Modélisation à différentes échelles de l'écoulement et du transport dans le sol et la zone non saturée[END_REF]. Initial values of were 𝜃 𝑠 obtained from experimental data and initial values of K s were obtained from the first experimental values of the K(h) curves. Concerning the DBM model, the initial values of parameters  1 , n 1 and  2 , n 2 were chosen according to the shape of the experimental water retention curve measured for each sample. These sets of initial values were then used as input to the RETC software and all the parameters were fitted (except l which was fixed at 0.5). As recommended by [START_REF] Sisson | An improved analysis of gravity drainage experiments for estimating the unsaturated soil hydraulic functions[END_REF] and [START_REF] Yates | Analysis of Measured, Predicted, and Estimated Hydraulic Conductivity Using the RETC Computer Program[END_REF], the relative weights of hydraulic conductivity data against retention data ranged between 0.1 to 1.0. As also recommended by [START_REF] Isch | A Comparative Study of Water and Bromide Transport in a Bare Loam Soil Using Lysimeters and Field Plots[END_REF], when a correlation higher than 90% was noted between two fitted parameters, the value of one of them was set to that obtained during the first fitting round and a second round was implemented.

The quality of the RETC fitting procedure was evaluated for each sample by the calculation of the coefficient of determination (R²) and the Akaike's Information Criterion (AICc) according to Eq. ( 10) (D'Emilio et al., 2018). When the number of experimental observations are small (for n/K < 40), the computation of the AICc is more appropriate than the standard AIC, given the increased relative penalty for small data sets and models with a high number of fitted parameters [START_REF] Laio | Model selection techniques for the frequency analysis of hydrological extremes[END_REF][START_REF] Pham | A New Criterion for Model Selection[END_REF]. It is worth noting that the lowest value of AICc indicates the best fitting efficiency.

𝐴𝐼𝐶𝑐 = 𝑛 × ln ( 𝑅𝑆𝑆 𝑛 ) + 2𝐾 + 2𝐾 (𝐾 + 1) 𝑛 -𝐾 -1 (10)
with n, the number of experimental observations (-), K, the number of fitted parameters (-) and RSS, the weighted residual sum of square (-).

The dispersivity () values of the soil and powdery limestone samples were obtained from elution experiments conducted on undisturbed columns by [START_REF] Viel | Etude des processus de transport des solutés hors équilibre physique: application à la zone non saturée des calcaires de Beauce[END_REF]. The dispersivity value of the calcareous sand samples was fixed to the mean of the values reported by Vanderborght and Vereecken (2007) for this type of texture. The dispersivity value of the limestone rock samples was taken from [START_REF] Kurotori | Measuring, imaging and modelling solute transport in a microporous limestone[END_REF].

The values obtained for all the parameters were then used to simulate water flow and bromide transport with HYDRUS-1D from 01/01/1966 to 31/12/2020 (55 years).

Procedure implemented for the estimation of the water travel time within the vadose zone

The water travel time (TT) through the VZ was estimated by a virtual tracing experiment, using bromide as an inert and conservative tracer [START_REF] Isch | A Comparative Study of Water and Bromide Transport in a Bare Loam Soil Using Lysimeters and Field Plots[END_REF] in HYDRUS-1D simulations. The bromide input concentration applied at the soil surface (9.7 g/L) was fixed in accordance with the maximum annual limit of nitrogen from livestock manure that can be applied per hectare in Nitrate Vulnerable Zones in the Centre-Val de Loire region (170 kg N/ha) and calculated following (11):

𝐶 0(𝐵𝑟) = 𝐶 0(𝑁) × 𝑀𝑊 (𝐵𝑟) 𝑀𝑊 (𝑁) (11)
with C 0(Br) , the bromide concentration in liquid phase (g/cm 3 ), C 0(N) , the nitrogen concentration in liquid phase (0.0017 g/cm 3 ), MW (Br) , the molecular weight of bromide (79.904 g/mol) and MW (N) , the molecular weight of nitrogen (14.0067 g/mol).

As recommended by [START_REF] Szymkiewicz | Numerical Analysis of Recharge Rates and Contaminant Travel Time in Layered Unsaturated Soils[END_REF], the solute transport simulation started with a "warm-up" period (of 2 years in our case) during which no solute was added to the soil. A single input of bromide was then applied on the 01/01/1968 with the solute concentration added in rainfall water (with R = 10 mm and ETP = 0 mm on the day of input). With the aim of studying TT variability according to the year of bromide input, simulations were repeated with a single input applied on the first day of each year (e.g., 01/01/1969, 01/01/1970, etc.).

As already recommended by some authors [START_REF] Fenton | Consequences of using different soil texture determination methodologies for soil physical quality and unsaturated zone time lag estimates[END_REF][START_REF] Fenton | Time lag: a methodology for the estimation of vertical and horizontal travel and flushing timescales to nitrate threshold concentrations in Irish aquifers[END_REF][START_REF] Szymkiewicz | Numerical Analysis of Recharge Rates and Contaminant Travel Time in Layered Unsaturated Soils[END_REF][START_REF] Vero | A framework for determining unsaturated zone water quality time lags at catchment scale[END_REF][START_REF] Vero | Consequences of varied soil hydraulic and meteorological complexity on unsaturated zone time lag estimates[END_REF][START_REF] Wang | Prediction of the arrival of peak nitrate concentrations at the water table at the regional scale in Great Britain[END_REF], the water travel time within the VZ was estimated using three indicators including (i) the arrival time of a bromide concentration of 1 mg/L at the WTL MAX (TT I ), which was viewed as the initial solute breakthrough (which approximately represents 1% of the bromide peak concentration and 0.01% of the bromide input concentration applied at the soil surface) and consequently reflects initial effects of regulation measures; (ii) the arrival time of bromide peak concentration at the WTL MAX (TT P ) (which approximately represents 1% of the input bromide concentration applied at the soil surface), which indicates the temporal trend of the arrival of the maximum concentration of a conservative solute down to the aquifer; (iii) the time after which a bromide concentration of 1 mg/L is no longer observed at WTL MAX (TT E ), as it represents the maximum residence time of the tracer (total exit) within the VZ and consequently the full effects of regulation measures.

Results & Discussion

Hydraulic properties and hydrodynamic parameters of the VZ samples

Experimental hydraulic properties

The hydraulic properties of a majority of the samples concerned by this work have already been discussed in [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF]. However, considering that the results obtained for all the samples have been further investigated and given that six new samples (S A , S B , S C , P A , R C and R D ) have been considered in this work, their experimental hydraulic properties are briefly discussed below.

The experimental hydraulic properties (water retention and hydraulic conductivity) of the S A and S B soil samples were relatively similar (Figure 3a). The S C sample displayed lower  s and K values between pF 0.0 and 2.0 and a higher water retention between pF 1.5 and 3.0, probably due to the slight increase in clay content observed with depth [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Ould Mohamed | Estimating Long-Term Drainage at a Regional Scale Using a Deterministic Model[END_REF]. The shapes of the water retention and hydraulic conductivity curves of the three fragmented powdery limestone samples (P B , P C and P D ) were close despite a  observed from pF 0.0 to 3.0 for sample P C (Figure 3b). The cryoturbated unit (P A ), which was described by Ould Mohamed and Bruand (1994), displayed higher K for pF < 1.5 and lower K between pF 2.0 and 3.0 than the three other powdery limestone samples (Figure 3b). Its  values were also in good agreement with those reported by [START_REF] Michot | Spatial and temporal monitoring of soil water content with an irrigated corn crop cover using surface electrical resistivity tomography[END_REF] over the range of pF 1.0-4.2. Water retention of calcareous sand interbeds (I A and I B samples) was the lowest among the soft VZ materials for pF > 1.0 (except I A vs P C ) (Figure 3c). I B sample displayed a lower water retention capacity than I A (except at saturation) because of its higher sand content [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF].

As already pointed out by [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF], the hydraulic properties of the hard limestone rock samples were highly heterogeneous (Figure 3d). Their values were the lowest of all the 𝜃 𝑠 VZ materials, except for the upstream weathered rock sample (R A ) which showed a coarse secondary porosity resulting in a strong decrease of  between saturation and pF 1.0. Their 𝜃 𝑠 values were also in good agreement with those obtained by other authors in a study conducted nearby [START_REF] Legchenko | Time-lapse magnetic resonance sounding measurements for numerical modeling of water flow in variably saturated media[END_REF]. The massive rock samples R B and R C displayed a nearly constant  between saturation and pF 3.0 in contrast to the downstream weathered rock samples

(R E and R F ), whose water retention decreased gradually over the same pF range. Compared to the weathered rock (R A , R E and R F ), the K values of the massive rock samples (R B , R C and R D )

were at first much lower near saturation (pF < 1.5) then quite close between pF 2.0 and 3.0.

Hydrodynamic parameters fitted by the models

The hydrodynamic parameters fitted with the DBM model using the RETC software provided the most satisfactory description of the experimental water retention and hydraulic conductivity data for all the VZ samples (Figure 3 As a reminder, when a correlation higher than 90% was noted between two fitted parameters, the value of one of them was set to that obtained during the first fitting round and a second round was implemented.

The compared to those fitted with the DBM model (Table 2). The experimental hydraulic conductivity curves of all the VZ samples were generally well described by the three models which gave good fits for 1.5 < pF < 3.0 (Figure 3 Consequently, the simulations made using the HYDRUS-1D software concerning the influence of the choice of the analytical model on the simulation of water flow and the estimation of TT within the VZ were focused on a comparison between the unimodal (MVG) and bimodal (DBM) expressions.

Table 2:

Parameters fitted with the DBM model using the RETC software and dispersivity () taken from the literature for the VZ materials.

Parameters highlighted in italics were fixed during the fitting procedure. Confidence intervals associated with parameters fitted using RETC are given in parentheses. CPL: cryoturbated powdery limestone, UW: upstream weathered, DW: downstream weathered. The standard deviations associated with the peak concentration and the mean TT I , TT P and TT E are given in parentheses. The years associated with the minimum and maximum TT I , TT P and TT E are given in brackets.

𝜽 𝒓 𝜽 𝒔 α 1 n 1 𝑲 𝒔 w 2 α 2 n 2 R² AICc λ Sample (Sampling) Lithology Depth [m] [cm 3 /cm 3 ] [cm -1 ] [-] [cm/d] [-] [cm -1 ] [-] [-] [-] [cm] S A (Soil

Simulation of water flow and estimation of travel time within the B2 profile

Since a majority of the samples were taken from the B2 profile (Table 2), the influence of the choice of the analytical model (MVG or DBM) on the simulation of water flow and the estimation of TT within the VZ was studied on this profile.

The B2 water content profiles simulated at the date of the drilling (14 March 2017) using the MVG and DBM models were in relatively good agreement with the experimental measurements made on the undisturbed core samples on that same date (R² ≥ 0.892, Figure 4b). The differences in  between the two models were rather small and the main disparities (< 0.05 cm 3 /cm 3 ) were observed within the soil (0-1 m deep), the powdery limestone facies (from 2.0 to 3.5 m deep and from 5.5 to 6.6 m deep) and the saturated zone (below 18 m deep).

Overestimations of experimental  by the models were noted in the powdery limestone material Since the Beauce limestone aquifer is known for its highly permeable soils and its vertical soil water flow regime [START_REF] Bruand | Variabilité de la recharge de la nappe de Beauce: Rôle de l'irrigation et des caractéristiques du sol[END_REF][START_REF] Flipo | Hybrid fitting of a hydrosystem model: Long-term insight into the Beauce aquifer functioning (France)[END_REF] A.3). The mean downward water flux simulated within the soil was higher for the DBM model than for the MVG model (Table A.

3).

This could be attributed to the higher cumulative ET R (actual evapotranspiration) simulated with the MVG model ( 30210 A.3), the bromide transport simulated with the MVG model from the base of the soil down to WTL MAX (14.6 m deep) was always slower than that simulated with the DBM model (Figure 6b-d).

Finally, the mean TT I , TT P and TT E (calculated from all the TT values simulated for each input year over the same period of time) were always higher for the MVG model than for the DBM model whatever the period of time considered (Table 3). The differences in mean TT values between the DBM and MVG models were significant and higher than 5 years for TT I , 9 years for TT P and 14 years for TT E (Table 3). The year of input corresponding to the minimum and maximum TT I , TT P and TT E were relatively close between the two models (Figure B.5 and Table 3). The most significant variations in TT simulated between the two models were noted for the years of input corresponding to a sudden increase in TT for the MVG model and which was not reproduced by the DBM model (e.g., 1978 and 1988 properties, the latter has been adopted for the comparative simulation of water flow and TT performed between three VZ profiles reconstituted with HYDRUS-1D from the three boreholes (B1, B2 and B3).

Simulation of the water flow and estimation of travel time within the vadose zone

(B1, B2 and B3 profiles)

Impact of the heterogeneity of the vadose zone profiles

Simulated  profiles with the DBM model at the date of the drilling (14/03/2017) for the three boreholes were in good agreement with the experimental measurements made on the undisturbed core samples on that same date (0.882 ≤ R² ≤ 0.925, Figure 4). Experimental  obtained for the three boreholes were generally slightly overestimated by the model, notably in P B and P C materials (from 1.5 to 4.0 m deep, depending on the borehole) and in limestone rock materials from 9.0 to 16.0 m deep. The main differences between experimental and simulated  (from 0.05 to 0.10 cm 3 /cm 3 ) within the three boreholes were found (i) in S B , P A and P B materials above 2.0 m deep; (ii) in R B materials around 8.0 m deep, at the transition between the powdery limestone (P D ) and the massive rock (R B ) materials [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF], which displayed high and low  s values respectively (Table 2); (iii) in R E and R F materials between 15.5 and 18.0 m deep, which corresponds to the weathered limestone rock that presented a highly heterogeneous structure with an increasing evolution toward karstification, fracturation and even silicification with depth [START_REF] Mallet | Integrated static and dynamic geophysical and geomechanical data for characterization of transport properties[END_REF].

The water regime within the soil seemed to be all the more dynamic as its thickness was small, as highlighted by the variations of the matric head, which suggested more pronounced drying and wetting cycles in the B2 profile than in B1 and B3 profiles (Figure 7a). The mean values of water content, matric head and water flux simulated for the three boreholes (Table A.4) indicated that the soil of the B3 profile was slightly wetter. It also displayed a lower vertical downward water flux than B1 and B2 profiles, although the three soil profiles were strictly identical down to 0.9 m deep (Table 1). This was attributed to the higher cumulative ET R simulated within the B3 profile (28296 mm), due to the larger thickness of its soil layer S C

(1.2 m, cf. Table 1) compared to B1 (0.4 m) and B2 (0.3 m). Indeed, the soil showed the highest water retention capacity among the VZ materials (Figure 3). This resulted in a visible differentiation of the bromide transport between the B3 and the B2 and B1 profiles with a slightly slower leaching of bromide down to 1.2 m deep in the B3 soil profile (Figure 8a). The differences in soil profiles and the influence of the meteorological conditions led to an arrival time for the peak of bromide ranging from 1.3 to 5.1 years at 1.2 m deep (Figure 8a), which represented a travel speed ranging from 0.23 to 0.94 m/yr.

The mean h values simulated by the DBM model in the VZ materials below 2.0 m deep for the three boreholes were higher than -100 cm (pF < 2.0) (Table A .4). This revealed a global wet water status over the simulation period (Figure 7b-d) and the occurrence of relatively high vertical downward water flux (due to K values close to saturation) within the VZ profiles (Table A.4). The strong lithological similarities between the B1 and B2 profiles down to 3.5 m deep (Table 1) led to a comparable bromide transport within the soil and the upper powdery limestone materials (P A , P B and P C ) of the two profiles (Figure 8a-b). The differentiation in bromide transport simulated between the B1 (the slowest) and B2 (the fastest) profiles from 3.5 m to 8.0 m deep (Figure 8b-d) was mainly caused by the differences in thickness of the powdery limestone materials P C and P D (3.8 m for B1 and 1.6 m for B2, cf. Table 1). The arrival time for the peak of bromide at 3.5 m and 8.0 m deep ranged from 5.0 to 11.0 years (Figure 8b) and from 11.1 to 18.8 years (Figure 8c), respectively. The peak concentration of bromide took between 3.8 and 5.9 years to migrate from 1.2 to 3.5 m deep (travel speed ranging from 0.39 to 0.63 m/yr) and between 5.0 and 7.8 years to migrate from 3.5 to 8.0 m deep (travel speed ranging from 0.58 to 0.89 m/yr).

The mean bromide peak concentration observed at WTL MAX for each VZ profile was between 94.4 (B3) and 104.8 mg/L (B2) (Table 4) which approximately represents 1% of the bromide input concentration (9.7 g/L). The lowest peak concentrations observed for B3 were caused by the higher global dispersivity of the profile (thicker soil and thinner powdery limestone, cf.

Table 1 and Table 2). The highest peak concentrations found for B2 profile were induced by its lower global dispersivity (thinner soil and thicker rock materials, cf. Table 1 andTable 2).

It is worth noting that the higher global dispersivity of the B3 profile resulted in the observation of bromide first concentrations at WTL MAX (TT I ) at times close to B2 and the total exit of the bromide (TT E ) at times close to B1 (Figure B.6 and Table 4). The mean indicators of TT for bromide to reach WTL MAX (14.6 m deep) were of 13.8 years for TT I (9.8 ≤ TT I ≤ 17.5 years), 20.9 years for TT P (14.6 ≤ TT P ≤ 25.3 years) and 31.5 years for TT E (26.6 ≤ TT P ≤ 39.0 years) (Table 4). The peak concentration of bromide was transported from the soil surface down to WTL MAX at a mean travel speed of 0.70 m/yr (min.: 0.58 m/yr and max.: 1.00 m/yr) and took between 5.8 to 8.3 years to migrate through the limestone rock facies from 8.0 m deep down to WTL MAX (travel speed ranging from 0.79 to 1.15 m/yr) (Figure 8d). For a given year of input, TT I and TT P were always the largest for the B1 profile and the smallest for the B2 profile and TT E was always the largest for the B3 profile and the smallest for the B2 profile ( The impact of the lateral variations in the VZ profiles (lithological differences at small scale: less than 15 m, cf. Figure 1) on the variations of the TT was relatively small. This is a consequence of the overall wet status observed over the simulation period within the VZ materials located below the soil (Figure 7b-d andTable A.4) which induced a leaching of bromide mainly impacted by K values close to saturation.

Impact of the meteorological conditions

As previously underlined by some authors [START_REF] Sprenger | Travel times in the vadose zone: Variability in space and time[END_REF][START_REF] Vero | Consequences of varied soil hydraulic and meteorological complexity on unsaturated zone time lag estimates[END_REF], meteorological conditions have a strong influence on the variation of TT through the VZ.

Maximum differences in the TT indicators according to input years for each VZ profile varied from 7.0 (B1) to 7.2 (B3) years for TT I , from 8.2 (B1) to 9.5 (B3) years for TT P and from 7.4 (B2) to 9.7 (B3) years for TT E (Table 4). The variations of the TT observed within a single VZ profile over the simulation period were higher than those observed between the three VZ profiles for a given year of input (Figure B.6 and Table 4). It appeared that, in the case of a globally wet VZ, the meteorological conditions had more influence on the variations of the TT within the VZ than the small scale lateral variations of its profile lithological composition.

Positive values of matric head were noted in several VZ materials located above the least Given the high lithological heterogeneity of the VZ materials, these perched water tables will however probably be less significant at the 3D field scale because of the lateral movement of water into open fractures or karst networks situated within the VZ [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Mallet | Integrated static and dynamic geophysical and geomechanical data for characterization of transport properties[END_REF].

The TT variations were tightly correlated with meteorological conditions and water table level The standard deviations associated with the peak concentration and the mean TT I , TT P and TT E are given in parentheses. The years associated with the minimum and maximum TT I , TT P and TT E are given in brackets.

Limitation of the work

The results obtained and discussed in this study focused on the simulation of the movement of water and bromide in the VZ of the Beauce limestone aquifer at the core scale. It has been well

established that interconnected fractures or karst networks may act as preferential flow paths leading to faster water and solute transport [START_REF] Chen | Water, nitrate and atrazine transfer through the unsaturated zone of the Chalk aquifer in northern France[END_REF][START_REF] Mcling | The Application of Radon for Mapping Open Fracture Networks in a Thin Vadose Zone[END_REF][START_REF] Wood | Unsaturated Flow through a Small Fracture-Matrix Network: Part 1. Experimental Observations[END_REF], notably after heavy rainfall events and thus only under or close to fully saturated conditions [START_REF] Wellings | Movement of water and nitrate in the unsaturated zone of Upper Chalk near Winchester, Hants., England[END_REF]. Some authors suggested that several tens of percent of the total flow in the VZ could be transported through fractures [START_REF] Smith | Water Movement in the Unsaturated Zone of High and Low Permeability Strata by Measuring Natural Tritium[END_REF]. As shown on Figure 9, natural fractures with karstification were observed at the field scale and were found to have been episodically subjected to intense water flow, under or close to fully saturated conditions [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Mallet | Integrated static and dynamic geophysical and geomechanical data for characterization of transport properties[END_REF]. The simulations carried out in this work showed the possible occurrence of perched water tables within the VZ which might generate preferential water flow through these open fractures networks, especially in the case of very rainy years. This would result in a faster solute transport within the fracture porosity which could be simulated by the use of a dual permeability model (Šimůnek et al., 2003). This modeling approach has proven to be useful for taking into account the contribution of preferential flow and representing the water flow and solute transport between fractures and matrix (Aguilar-López et al., 2020; [START_REF] Köhne | Estimating the hydraulic functions of dual permeability models from bulk soil data[END_REF]. Recent work also showed that, even if a dual-porosity model was capable of adequately describing the arrival of the measured initial concentrations and the shape of the breakthrough curves of bromide at the base of a loamy soil, only the dual permeability model used with the inclusion of a mobileimmobile water component (MIM) provided a good ability to identify the arrival of peak concentrations (Varvaris et al., 2021b). The identification and characterization of these open fractures networks and the understanding of their functioning will consequently be crucial for complementing our knowledge of water flow, travel time and contaminant transport through the vulnerable Beauce limestone aquifer.

Summary and Conclusion

In this study, the numerical simulation of water flow and the estimation of water travel time through the highly heterogeneous VZ of a limestone aquifer was undertaken over a period of 55 years. The parameterization of the model (HYDRUS-1D) was based on hydraulic properties laboratory measurements performed on fifteen samples representative of the lithologies encountered within the VZ using unimodal (MVG) and bimodal (DBM) approaches. The results obtained in this work highlighted that the experimental hydraulic properties of the samples were more accurately described by the DBM model. The latter allowed a precise description of the bimodal characteristic of the VZ materials, i.e., macro-and micro-porosity for the soil, powdery limestone and calcareous sand samples, and matrix and (micro-)fissures porosity for the limestone rock. The choice of the analytical model used to describe the VZ hydraulic properties (MVG or DBM) showed a significant influence on the simulation of the water flow and TT within the VZ. The results of the simulations carried out with the two models displayed strong discrepancies regarding the water flow dynamics (matric head and water content variations) through the VZ and showed large differences in TT (occasionally higher than 10 years for TT I and TT P and 15 years for TT E ). These disparities were mainly attributed to the differences observed in the water flow dynamics simulated in the soil and the higher hydraulic conductivity values obtained with the DBM model for the VZ materials located below the soil over the range of matric head (h > -250 cm) observed during the simulation period.

A mean TT I of 13.8 years (min. 9.8 -max. 17.5), TT P of 20.9 years (min 14.6 -max. 25.3) and TT E of 31.5 years (min 26.6 -max. 39.0) were found for the VZ of the Beauce limestone aquifer with simulations carried out using the DBM model. The vertical heterogeneity of the properties of the VZ materials showed a strong influence on the simulation of the water flow and the estimation of TT. Significant differences were found in the travel speed of the peak concentrations of bromide through the VZ materials, which was the fastest in the limestone rock facies (max. of 1.15 m/yr) and the slowest in the powdery limestone facies (max. of 0.63 m/yr). Despite the lithological differences observed between the three VZ profiles, only slight variations were found regarding the water flow dynamics and the TT (max. 4.7 years whatever the indicator) between the three boreholes. This small impact of the lateral heterogeneities was mainly attributed to the global wet status found in the VZ materials located below the soil, which led to hydraulic conductivities always close to saturation and consequently hardly noticeable differences between the three VZ profiles. The meteorological conditions have shown a large influence on the variations of the TT through the VZ (max. 9.7 years whatever the indicator). Together with the rise of the mean annual temperature and ETP over the simulation period, an increase in TT was noted since the late 1970s, becoming even more apparent after the late 1990s and may be a consequence of global warming. In the case of this mostly unconfined aquifer, this could mean that the impact of the regulatory measures that have been put into effect since the beginning of the 1990s (e.g., the Nitrates Directive, European Commission, 1991) for the preservation of groundwater resources may not yet be clearly identifiable today and that it could take progressively longer to identify such impacts in the future.

The results presented in this work were obtained with fifteen samples selected along the whole VZ profile, and represent a first step towards capturing and deciphering in depth the high heterogeneity of the VZ materials and characterizing the complexity of the transfer through the VZ of the Beauce limestone aquifer [START_REF] Aldana | Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer[END_REF][START_REF] Mallet | Integrated static and dynamic geophysical and geomechanical data for characterization of transport properties[END_REF]. It will indeed be necessary to improve our knowledge about the vertical and lateral flow of water and solute through the whole porosity of this highly heterogeneous VZ, by conducting extensive studies at larger scales, combining 3D experimental observations and numerical simulations with dual permeability approaches. To this end, noninvasive hydrogeophysical monitoring tools seem to represent a worthwhile solution for mapping geological facies properties and linking measured geophysical parameters to the hydraulic properties of the VZ materials [START_REF] Binley | The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales: The Emergence of Hydrogeophysics[END_REF][START_REF] Parsekian | Multiscale geophysical imaging of the critical zone[END_REF][START_REF] Romero-Ruiz | A Review of Geophysical Methods for Soil Structure Characterization[END_REF]. When combined with the use of wellestablished or advanced in situ methods in the VZ (i.e., tracer tests and fracture networks characterization) [START_REF] Mcling | The Application of Radon for Mapping Open Fracture Networks in a Thin Vadose Zone[END_REF][START_REF] Sprenger | Travel times in the vadose zone: Variability in space and time[END_REF] or in the saturated zone (pumping, slug or tracer tests) [START_REF] Binet | Water exchange, mixing and transient storage between a saturated karstic conduit and the surrounding aquifer: Groundwater flow modeling and inputs from stable water isotopes[END_REF][START_REF] Mosthaf | Conceptualization of flow and transport in a limestone aquifer by multiple dedicated hydraulic and tracer tests[END_REF], such tools should lead to a refined hydraulic characterization of the materials and a more accurate assessment of the travel time of water and the fate of contaminants through this highly vulnerable limestone aquifer of major importance.

content (cm 3 /cm 3 ) and water flux (cm/d) obtained at the observations nodes defined within the B1, B2 and B3 profiles with the DBM model.

  ) have been drilled in March 2017 from 0 to 20 m deep. The soil formation (0.0-1.5 m deep) was drilled by hydraulic percussion while the other VZ materials (1.5-20.0 m deep) were taken by rotary drilling with a three-annular compartments corer to avoid sample contamination by the drilling fluid (water).

Figure 1 :

 1 Figure 1: Localization of the experimental site, distance between the three cored boreholes

Figure 2 :

 2 Figure 2: Annual mean temperature, rainfall and ETP, and daily water table level observed from 1966 to 2020.

For

  the soil samples (S A , S B and S C ), initial values of parameters , , , n and K s were taken 𝜃 𝑟 𝜃 𝑠 from Ould Mohamed et al. (1997). Experimental K s values have already been obtained in the laboratory (constant-head method) by Ould Mohamed et al. (1997) for soil samples taken between a few tens and hundreds of meters from our experimental site. Since the authors have shown little spatial variability in the hydraulic properties of soils taken within the same pedological layer, the values of K s of the soil samples S A , S B and S C were fixed to the experimental values given by Ould Mohamed et al. (1997) for three of their soil samples taken within the same layers and at the same depths. For the powdery limestone (P A , P B , P C and P D ) and the calcareous sand samples (I A and I B ), initial values of parameters , , , n and K s were 𝜃 𝑟 𝜃 𝑠

  ) compared to those obtained with the MVG (Figure B.1) and BBC (Figure B.2) models. Indeed, R² (≥ 0.938) and AICc (≤ -104.1) values calculated with the DBM model for each sample were always higher and lower than those obtained with the MVG (except R² for R C and R E ) and BBC (except R² for R E ) models, respectively (

  quality of the description of the water retention curves was more impacted by the choice of the analytical model than the hydraulic conductivity curves. The experimental water retention properties of the soft VZ materials, that mostly displayed a clear bimodal behavior (notably S A , S B , P A , P D , I A and I B ), were poorly depicted by the MVG and BBC models, especially for 0.5 < pF < 1.5 where retention was overestimated by the models and for 2.5 < pF < 4.0 where it was underestimated (Figure B.1a-c and Figure B.2a-c). Higher values of parameter  s were obtained with the MVG (

  , Figure B.1 and Figure B.2). The main differences in the fitted K(h) curves among the three models were found in the description of the hydraulic conductivity near saturation (pF < 1). Indeed, the K values obtained by the DBM model near saturation were frequently higher or equal than those given by the MVG and BBC models (except S C for MVG and BBC and I A for BBC). This resulted in K s values always higher or equal for the DBM model compared to the MVG (except P A , I A and R B ) and BBC (except I A ) models (

Figure 3 :

 3 Figure 3: Comparison between experimental and fitted (DBM model) water retention ((h), up) and hydraulic conductivity (K(h), down) curves

Figure 4 :

 4 Figure 4: Comparison between measured water content profiles for B1 (a), B2 (b) and B3 (c) at the date of the drilling (14 March 2017) and those

Figure 5 :

 5 Figure 5: Simulated values of matric head obtained between 1991 and 1994 (dry VZ profile

Figure 6 :

 6 Figure 6: Bromide concentration simulated at 0.9 (a), 3.5 (b), 6.6 (c) and 14.6 m (WTL MAX ) (d) deep by the MVG and DBM models in the B2

P

  B (2 m deep) and in limestone rock materials (R C , R D and R E ) from 10 to 15 m deep. Underestimations of experimental  by the models were observed in massive (R B , 8.5 m deep) and downstream weathered limestone rock materials (R F , 16 m deep).

  , no run off was observed during the whole simulation period (55 years and a total rainfall of 35180 mm). The matric head values simulated by the DBM model in each soil layer displayed more fluctuations than those of the MVG model (Figure 5 and Figure B.3). The soil water status simulated by the DBM model was drier (lower mean values of h) and characterized by more pronounced drying and wetting cycles than for the MVG model (Figure B.3 and Table

  Figure B.2a).As a result, the bromide transport simulated within the soil by the DBM model

  for TT I ; 1981 and 1982 for TT P ) (Figure B.5).The simulation of the water flow and the estimation of the TT made using unimodal(MVG) and bimodal (DBM) models fitted to the same experimental hydraulic properties can show large disparities (9.8 ≤ TT I ≤ 23.1 years, 14.6 ≤ TT P ≤ 34.9 years, 27.3 ≤ TT E ≤ 46.5 years). The DBM model allows to better account for the impact of the macroporosity in the soft VZ materials and the presence of (micro-) fissures (< 1 mm thick) in the hard limestone rock materials compared to the MVG model. In view of all these findings and since the hydrodynamic parameters fitted with the DBM model provided the most satisfactory description of the experimental hydraulic

  Figure B.6). Differences of up to 3.6 (B1-B2, 1979), 3.5 (B1-B2, 1975) and 4.7 years (B3-B2; 1977) were observed regarding TT I , TT P and TT E respectively (Figure B.6).

  permeable materials (P B , P C and R C ) (Figure3b,d and Figure7b-c), and appeared in years showing annual R-ETP balance close or above 0 mm(such as 1978-1984, 1988 or 2001-2002, cf. Figure2). This indicates that perched water tables may have occurred in the VZ profiles, caused by accumulation of water above these weakly permeable materials. The building up of perched water tables were highlighted in the soft VZ materials up to a depth of 0.7 m and in the limestone rock from about 11 m deep up to a depth of 7 m (Figure7c-d and Figure B.4b).

  data that were used as boundary conditions over the simulation period. The relatively slow leaching of bromide observed in the VZ (mean travel speed: 0.70 m/yr) could be partly explained by the low mean annual rainfall (639.6 mm) and the high mean annual ETP (813.7 mm) observed over the simulation period (Figure2) and also by the low hydraulic conductivity (even near saturation) of the powdery limestone (P B , P C and P D ) and the massive rock (R B and R C ) materials (Figure3b,d). The relatively faster leaching of bromide (corresponding to low TT) observed for inputs made between 1968 and 1975 (Figure B.6) could be explained by the meteorological conditions observed between 1977 and 1984, when annual R-ETP balances were close or above 0 mm and when WTL was rising (Figure 2). The global slowdown of the leaching of bromide (corresponding to a general increase in TT) observed for inputs made from the late 1970s (Figure B.6) was strongly influenced by the small annual rainfall and the lowering of the WTL observed from 1989 to 1998, 2003 to 2011 and 2015 to 2020 (Figure 2) and may be a consequence of global warming. Indeed, it is worth noting that from 1989 to 2020, only two years of input showed a positive annual R-ETP balance (1999 and 2001) while an increase in mean annual temperature (of 1.0°C) and ETP (of 98 mm) was noted compared to the 1968-1988 period (Figure 2). As highlighted by the trends observed for TT I since the end of the 1990s (Figure B.6), a further increase in TT for inputs made after 2000 is expected as a response to the global increase in mean annual temperature and ETP observed since 2015 (maximum of 12.8°C and 1044 mm observed in 2020), associated with a strongly negative annual R-ETP balance (Figure 2).

Figure 7 :

 7 Figure 7: Simulated values of matric head obtained between 1991 and 1994 (dry VZ profile with WTL MIN = -22.21 m on 19/08/1992, up) and

Figure 8 :

 8 Figure 8: Bromide concentration simulated at 1.2 (a), 3.5 (b), 8.0 (c) and 14.6 m (WTL MAX ) (d) deep by the DBM model in the three VZ profiles

Figure 9 :

 9 Figure 9: Optical imaging of five new boreholes (B5, B6, B7, B8 and B9 drilled in Spring 2020

Figure B. 2 :

 2 Figure B.2: Comparison between experimental and fitted (BBC model) water retention ((h), up) and hydraulic conductivity (K(h), down)

Figure B. 3 :

 3 Figure B.3: Simulated values of matric head obtained between 1966 and 2020 within the soil

Figure B. 4 :

 4 Figure B.4: Simulated values of matric head obtained between 1966 and 2020 at the

Figure B. 5 :

 5 Figure B.5: TT I (dotted line, down), TT P (solid line) and TT E (dotted line, up) estimated from the simulation of bromide transport within the B2

Figure B. 6 :

 6 Figure B.6: TT I (dotted line, down), TT P (solid line) and TT E (dotted line, up) estimated from the simulation of bromide transport within the B1,

Figure 1 :

 1 Figure 1: Localization of the experimental site, distance between the three cored boreholes (B1, B2 and B3) and position of the soil pit.

Figure 2 :

 2 Figure 2: Annual mean temperature, rainfall and ETP, and daily water table level observed from 1966 to 2020.

Figure 3 :Figure 4 :

 34 Figure 3: Comparison between experimental and fitted (DBM model) water retention ((h), up) and hydraulic conductivity (K(h), down) curves

Figure 5 :

 5 Figure 5: Simulated values of matric head obtained between 1991 and 1994 (dry VZ profile with WTL MIN = -22.21 m on 19/08/1992, up) and

Figure 6 :

 6 Figure 6: Bromide concentration simulated at 0.9 (a), 3.5 (b), 6.6 (c) and 14.6 m (WTL MAX ) (d) deep by the MVG and DBM models in the B2

Figure 7 :

 7 Figure 7: Simulated values of matric head obtained between 1991 and 1994 (dry VZ profile with WTL MIN = -22.21 m on 19/08/1992, up) and

Figure 8 :

 8 Figure 8: Bromide concentration simulated at 1.2 (a), 3.5 (b), 8.0 (c) and 14.6 m (WTL MAX ) (d) deep by the DBM model in the three VZ profiles

Figure 9 :

 9 Figure 9: Optical imaging of five new boreholes (B5, B6, B7, B8 and B9 drilled in Spring 2020 at a few tens meters from B1, B2 and B3)

Figure B. 1 :

 1 Figure B.1: Comparison between experimental and fitted (MVG model) water retention ((h), up) and hydraulic conductivity (K(h), down) curves

Figure B. 2 :

 2 Figure B.2: Comparison between experimental and fitted (BBC model) water retention ((h), up) and hydraulic conductivity (K(h), down) curves

Figure B. 3 :

 3 Figure B.3: Simulated values of matric head obtained between 1966 and 2020 within the soil layers of the B2 profile with the MVG (a) and DBM

Figure B. 4 :

 4 Figure B.4: Simulated values of matric head obtained between 1966 and 2020 at the observation nodes defined for the B2 profile with the MVG

Figure B. 5 :

 5 Figure B.5: TT I (dotted line, down), TT P (solid line) and TT E (dotted line, up) estimated from the simulation of bromide transport within the B2

Figure B. 6 :

 6 Figure B.6: TT I (dotted line, down), TT P (solid line) and TT E (dotted line, up) estimated from the simulation of bromide transport within the B1,

Table 2 ,

 2 Table A.1 and Table A.2).  r values were frequently set to 0 by the RETC software during the fitting procedure (Table 2, Table A.1 and Table A.2). Some high values of correlations (> 90%) have occasionally been noted between K s et  during the fitting made with the MVG

and BBC models and between K s ,   and n 1 and  2 and n 2 for the DBM model (Table 2, Table A.1 and Table A.2).

Table 2 ,

 2 Table A.1 and TableA.2). MVG and BBC models still allow a decent representation of the hydraulic properties of the VZ materials, with R² ≥ 0.866 (except for S C and R A ) and AICc ≤ -85.6 for MVG and R² ≥ 0.849 (except for S C and R C ) and AICc ≤ -75.9 for BBC (TableA.1 and TableA.2). However, the quality of the fit was better for the MVG model as the AICc values computed for each sample were generally lower than those obtained for the BBC model (except for I B and R A ).

	The

Table 3 :

 3 Results obtained for the simulation of bromide transport down to WTL MAX(-14.60 m) for the B2 profile according to the choice of the analytical model: peak concentration and mean, minimum and maximum indicators of travel time (TT I , TT P and TT E ).

		1968-1999		1968-1986		1968-1975
	Analytical	Mean TT I	Min TT I Max TT I	Peak Concentration	Mean TT P	Min TT P Max TT P	Mean TT E	Min TT E Max TT E
	Model
		y	y	mg/L	y	y	y	y
	MVG	18.5 (3.1)	13.1 [1974] 23.1 [1989]	92.8 (4.0)	28.4 (3.2)	24.4 [1979] 34.9 [1986]	44.7 (1.1)	43.4 [1974] 46.5 [1968]
	DBM	13.2 (2.1)	9.8 [1974] 16.9 [1983]	105.9 (5.9)	18.8 (2.2)	14.6 [1974] 22.7 [1985]	29.7 (1.7)	27.3 [1975] 32.2 [1968]

Table 4 :

 4 Results obtained for the bromide transport simulated down to WTL MAX (-14.60 m) by the DBM model for each VZ profile (B1, B2 and B3): peak concentration and mean, minimum and maximum indicators of travel time (TT I , TT P and TT E ).

	1968-2001	1968-1998	1968-1982

  Threeheterogeneous VZ profiles were reconstituted in HYDRUS-1D based on three boreholes drilled in March 2017. Daily rainfall, ETP and water table level data were used as inputs for the implementation of a virtual bromide tracing experiment made over the 55 years. The arrival time of the initial breakthrough (TT I ), the peak (TT P ) and the total exit of the bromide concentration (TT E ) at the maximum water table level were considered as suitable indicators of water travel time within the VZ of the Beauce limestone aquifer.

Acknowledgements

This study was conducted within the framework of the O-ZNS project, which is part of the PIVOTS program (https://plateformes-pivots.eu/o-zns/?lang=en). We gratefully acknowledge the financial support provided to the PIVOTS project by the Centre-Val de Loire region (ARD 2020 program and CPER 2015-2020) and the French Ministry of Higher Education and Research (CPER 2015 -2020 and public service subsidy to BRGM). This operation is co-funded by the European Union. Europe is committed to the Centre-Val de Loire region with the European Regional Development Fund (FEDER). This research work was co-funded by the Labex VOLTAIRE (ANR-10-LABX-100-01). We would like to warmly acknowledge the reviewers and the associate editor for their useful comments which allowed us to considerably improve the content and highlight novelties of our manuscript.

Appendix A

for the vadose zone materials.

Parameters highlighted in italics were fixed during the fitting procedure. Confidence intervals associated with parameters fitted using RETC are given in parentheses.


est destinée au dépôt

Shift of belowground traits during perennial grain growth phase in annual grains crop rotation

Olivier Duchene, Florian Celette, Ana Barreiro, Linda-Maria Dimitrova A range of 'agroecological practices' are currently proposed to increase the sustainability of intensive grain systems. However, little attention has been paid to the use of perennial grains in the crop successions. The use of the perennial grain 'intermediate wheatgrass' (Thinopyrum intermedium (Host) Barkworth & D.R. Dewey) may have the potential to sustain soil fertility through the development of an extensive root system beneficial to a range of soil functions. In the context of cereal grain crop rotation, we compared young stands of intermediate wheatgrass to annual grains during two growing seasons, with the aim to determine how the rooting system of a perennial crop differs from an annual crop along a 1.6 m deep soil profile. Our results indicate that the two-year-old intermediate wheatgrass promotes denser and deeper rooting systems with proportionally higher root distribution and biomass downward the soil profile. In the first growing season, higher root biomass was observed under the perennial grain than under the annual wheat at 0-10 cm depth. In 12 months, between the first and second spring growing seasons, perennial belowground biomass increased by 111%. Consequently, root biomass was larger under the perennial treatment at each soil layer, leading to 5.3 average more tons of roots (dry weight) per hectare compared to the annual grain at the same time. Looking at the root to shoot ratio, perennial roots represented 38.6% of the aerial biomass produced during the first growing season, compared with 24.4% for the annual grain; while reaching 55.3% the second year, compared to 22.9% for the annual grain. Higher level of soil colonization were also observed under the perennial grain between 10 and 160 cm depth, including greater evenness of the roots vertical distribution and improved soil structure between 10 and 20 cm. From the first growing season, the perennial grain also harbors a suite of root traits (specific root length, root diameter, root tissue density) typical of a more resource-conservative strategy, and more belowground-oriented resource allocation. Additionally, soil fungal biomass indicators, including those of arbuscular mycorrhizal fungi, were enhanced under perennial grain, which potentially show improved soil quality under such management. Overall, this study brings evidence that grain agriculture would benefit from the possibility to rely on deeper and long-lived root systems to manage soils. The periodic use of a perennial phase in the crop rotation has the potential to improve soil functioning in the long term, while allowing both forage and grain production.

Keywords: perennial grains, rooting system, root traits, soil microbial indicators, soil quality

  Mårtensson, Grégoire T. Freschet, Christophe David To cite this version: Olivier Duchene, Florian Celette, Ana Barreiro, Linda-Maria Dimitrova Mårtensson, Grégoire T. Freschet, et al.. Shift of belowground traits during perennial grain growth phase in annual grains crop rotation. ESA Congress, 2020, Sevilla, Spain. hal-03653274 1-4 September 2020, Sevilla, Spain 1 SHIFT OF BELOWGROUND TRAITS DURING PERENNIAL GRAIN GROWTH PHASE IN ANNUAL GRAINS CROP ROTATION Olivier Duchene 1 , Florian Celette 1 , Ana Barreiro 2 , Linda-Maria Dimitrova Mårtensson 2 , Grégoire T. Freschet 3 , Christophe David 1 1 ISARA, Agroecology and Environment Research Unit, 23 Rue Jean Baldassini, 69364, Lyon cedex 07, France. olduchene@isara.fr 2 Swedish University of Agricultural Sciences, Department of Biosystems and Technology P.O. Box 103, SE-230 53 Alnarp, Sweden 3 Station d'Ecologie Théorique et Expérimentale, CNRS, Université Toulouse III, 2 route du CNRS, 09200 Moulis, France