Relationship between hydraulic properties and material features in a heterogeneous vadose zone of a vulnerable limestone aquifer

Limestone aquifers constitute important drinking water resources but are vulnerable to groundwater contaminations. They are also known to be often heterogeneous. However, the hydraulic properties of their entire vadose zone and the relationship of these hydraulic properties with their mineralogical and geochemical characteristics remain rarely studied. The hydraulic properties of soft (soil, powdery limestone, and calcareous sand) and hard (limestone rock) materials sampled throughout the entire vadose zone profile of the Beauce limestone aquifer (France) were determined with the multistep outflow method applied using a triaxial system. Physical, mineralogical, and geochemical analysis brought valuable information about factors contributing to the heterogeneity of the hydraulic properties and helped understanding water flow pathways within the vadose zone. The hydraulic properties of the soft materials were strongly related to their physical properties but also to the proportion and the nature of clay minerals. The massive rock presented a few thin microfissures where calcite was replaced by phyllosilicates, thus increasing its water retention capacity. The low hydraulic conductivity of the massive rock could explain the occurrence of perched water tables in the vadose zone. The weathered rock displayed fissures and vugs and was also characterized by the presence of phyllosilicates. The development of a secondary porosity and footprints of water–rock interactions and mass transfers highlighted the transformation of the limestone rock into a more permeable material. This study also pointed out the need to characterize the impact of natural fractures observed at the core scale on the preferential water flow at the field scale.


INTRODUCTION
Limestone aquifers have an overriding interest in many countries, especially in the geographic areas under intensive agriculture, as they constitute important drinking water resources but are also vulnerable to groundwater pollutions (Aisopou et al., 2015). Even though limestone aquifers are regularly considered as a single unit, the geology of such aquifers is, in fact, often highly heterogeneous as the upper units may be crushed or weathered (soft materials), whereas the lower units are generally unaltered and/or fractured (hard materials) (Dar et al., 2017;Williams et al., 2006). The knowledge of water flow and solute transport through limestone aquifers is consequently limited and perceived as especially complex, given the differences in the properties of their geological units (Mosthaf et al., 2018). The hydraulic characterization of these aquifers is commonly carried out in situ and under saturated conditions with slug and pumping tests (Ackerer & Delay, 2010;Audouin & Bodin, 2008) or characterization methods such as advanced tracer tests (Maurice et al., 2012;Mosthaf et al., 2018;Somogyvári et al., 2016). The vadose zone (VZ), which extends from the soil surface down to the aquifer, is considered as a key interface of the critical zone and is of major importance in understanding and managing contamination problems for the protection and preservation of groundwater resources. Although many authors highlighted that the evaluation of the hydraulic properties (water retention and hydraulic conductivity) of the VZ materials is a cornerstone in the study of water flow, heat, and mass transfers in the critical zone (Lekshmi et al., 2014;Nolz, 2016;Robinson et al., 2008;Shahbazi et al., 2020;Singh et al., 2018;Vereecken et al., 2008;Whalley et al., 2013), there is, to our best knowledge, no study focused on the characterization of the saturated and unsaturated hydraulic properties of the whole VZ materials of a limestone aquifer. A wide variety of laboratory methods exist for the determination of the saturated and unsaturated hydraulic properties of VZ materials (Gribb et al., 2004;Masrouri et al., 2008), such as the pressure plate apparatus (Bittelli & Flury, 2009), the centrifuge method (Caputo & Nimmo, 2005), the evaporation method (Schindler & Müller, 2006;Wendroth et al., 1993), the rigid-wall (Chapuis, 2004) and the flexible-wall permeameters (triaxial cell) (Huang et al., 1998;Samingan et al., 2003), the one-step (van Dam et al., 1992) and multistep outflow (MSO) method (Eching et al., 1994;van Dam et al., 1994), each with their advantages and disadvantages (Bordoni et al., 2017;Chapuis, 2012;Cresswell et al., 2008;Schelle et al., 2010;van den Berg et al., 2009).
Some of the works dedicated to the characterization of the VZ try to determine relationships that link physical properties (bulk density, texture, etc.) to hydraulic properties using pedotransfer functions (Assouline & Or, 2013;Patil & Singh, 2016;Vereecken et al., 2010Vereecken et al., , 2016Zhang & Schaap, 2019). However, very few aim at understanding the mechanisms at the origin of these relationships, especially those linking hydraulic properties to material characteristics (physical, mineralogical, and geochemical properties).
In this study, the hydraulic properties of samples of five soft and five hard VZ materials taken from three cored boreholes have been determined using the MSO method

Core Ideas
• Vadose zone (VZ) hydraulic properties have been measured on core samples. • VZ materials displayed high heterogeneity in their properties. • Hydraulic properties of VZ materials are related to their microscopic and macroscopic features. • Mineralogical and geochemical analysis reveal water flow pathways. • Perched water table may occur in the VZ profile.
applied by a triaxial system. Physical, mineralogical (X-ray diffraction), and qualitative geochemical analyses (scanning electron microscopy/energy dispersive X-ray spectroscopy [SEM/EDS]) have also been made in order to investigate the relationships between macroscopic or microscopic features and hydraulic properties with the aim of revealing water flow pathways. In the meantime, numerical simulations have been performed with the HYDRUS-1D model using the measured hydraulic properties to provide key information on water flow throughout the whole VZ of the aquifer.

Description of the Beauce aquifer system
The Beauce aquifer system is located in the center of France and extends over 9,700 km 2 (de Frutos Cachorro et al., 2017) between the Seine river (northeastern part) and the Loire river (southwestern part) ( Figure 1). This mostly unconfined aquifer is one of the largest groundwater reservoirs in France and is mainly composed of Cenozoic limestone from upper Oligocene to lower Miocene whose origin is mostly lacustrine (Ménillet & Edwards, 2000). Its thickness and its topography ranges from 10 to 200 m and from 70 up to 190 m, respectively (Flipo et al., 2012). The land use consists essentially of agriculture (74%), 50% of which is irrigated (Lejars et al., 2012). Over the past decades, human activities and mostly intensive agriculture have affected the groundwater quality with nitrates and pesticides, notably atrazine and bentazon, at concentrations that still exceed regulatory limits (DDT, 2016).
In this context, the O-ZNS (Observatory of transfers in the VZ) project aims to understand and quantify mass and heat transfers throughout the heterogeneous VZ of the Beauce limestone aquifer, by integrating observations over a wide range of spatial (from nanometer to meter) and temporal (from minutes to decades) scales by means of a large access well (depth = 20 m, diameter = 4 m) surrounded by several Vadose Zone Journal F I G U R E 1 Representation of the Beauce aquifer system (modified from Flipo et al., 2012). The red dot indicates the location of study site F I G U R E 2 Location of the experimental site and distance between the three cored boreholes (B1, B2, and B3) boreholes in order to combine broad characterization and focused specific monitoring techniques.

Location of the study site, climate, and water table level data
The experimental field site is situated at Villamblain, 30 km northwest of the city of Orléans (France) (coordinates: X = 48˚1′5.131′′; Y = 1˚34′55.333′′) (Figure 2). The climate of the study site is continental-temperate with an annual average temperature of 10.5˚C.
Meteorological data were collected from the Bricy weather station located about 20 km east of the study site. From 1966 to 2017, the mean annual rainfall was 642.7 mm with a minimum of 413.3 mm observed in 1990 and a maximum of 929.8 mm observed in 2001. The mean annual Penman-Monteith potential evapotranspiration (ETP) (Monteith, 1965)

2.3
In situ sampling of the VZ Three cored boreholes (B1, B2, and B3) separated by about 10 m one from the other were drilled in March 2017 down to the depth of 20 m in order to sample the entire VZ of the Beauce aquifer ( Figure 2). Core samples were drilled into polycarbonate liners (96-mm i.d.) with a three annular compartments corer to avoid sample contamination by the drilling fluid. The first 2 m of the VZ were drilled by hydraulic percussion and the 18 m below were sampled by rotary drilling using water as drilling fluid. The lithologies encountered throughout the VZ (0-20 m deep) were described by visual examination of the undisturbed core samples. A soil with a thickness between 1.0 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 (Duval & Isambert, 1992;Michot et al., 2003) according to the French soil reference system (Baize & Girard, 2009) or a Hypereutric Cambisol according to the World Reference Base for Soils (IUSS Working Group WRB, 2015) and has a silty loam texture according to USDA (2012) classification. The soil lays on a highly heterogeneous Miocene lacustrine fragmented powdery limestone facies (referred to as "powdery limestone") that was cryoturbated in its upper part during the Quaternary (Michot et al., 2003;Ould Mohamed & Bruand, 1994) and which also contains calcareous sand interbeds. The thickness of this layer is from 5.2 to 6.7 m. Finally, the last main stratigraphic facies is a 12.2-to 13.4-m-thick massive, fractured, or weathered hard limestone rock (Pithiviers limestone) (Schnebelen et al., 1999).

Hydraulic properties measurements
Undisturbed samples were re-cored in the laboratory at different depths based on the lithological description made on the cored samples. The sample diameter was adapted to the experimental device (70-mm diam.) and most of the samples had a thickness of 30 mm. Rock samples were ranging from 25 to 50 mm in thickness due to difficulties in sampling this massive fractured/weathered material. Three samples from the soil formation (S A , S B , and S C ), three samples from the powdery limestone facies with different textures (P A , P B , and P C ) and two samples from the calcareous sand interbeds (I A and I B ) were re-cored. Since rock materials presented more variations in their appearance (fissuring, fracturing, upstream/downstream weathering), five representative rock samples (R A , R B , R C , R D , and R E ) were selected to capture the impact of these variations on hydraulic properties. All samples used for hydraulic property measurement were undisturbed except for I B (sand deposit), which was repacked to the observed field bulk density. In order to estimate the unsaturated hydraulic properties of the VZ materials, the MSO was applied using a flexible-wall permeameter (triaxial system) ( Figure 3) because of the possibility (a) to simultaneously determine the water retention and hydraulic conductivity curves of soft and hard materials without the need to change the method; and (b) to reproduce the field conditions generated by the lithostatic pressure of the overlying layers by applying a confining stress to the sample.
The measurement system consists of a triaxial cell of 230-mm diam. and 530-mm height with a maximum confinement pressure of 3.5 MPa (GDS Instruments). A high-flow ceramic plate with a 100-kPa air-entry value was bonded into the pedestal so that large unsaturated water pressure gradients can be maintained within the sample when unsaturated experiments were carried out. The sample was isolated from the cell via a flexible membrane, which allowed the application of a confinement pressure to the sample. The pedestal was connected to four 2-MPa pressure-volume controllers (PVCs) via five high-pressure Swagelok valves. Three PVCs (Back-, Base-and Cell-PVC) were designed to control water pressures and one PVC (Gas-PVC) controlled air pressure (Figure 3).
The MSO experiments were carried out using the triaxial system and by following the same experimental procedure applied by Eching et al. (1994). This technique allows measuring the hydraulic properties of porous samples by applying an isotropic consolidation (Huang et al., 1998) and varying the matric potential by steps (Schelle et al., 2010). Unsaturated hydraulic properties were measured within the pF value range 0.5-3.0.
The liquid phase was a solution of deionized water with 0.01 M of CaCl 2 and three drops of 10% diluted NaOCl per liter. The solution was deaerated before the experiment. The gas (nonwetting) phase was N 2 with 99.9% purity. The temperature of the room was continuously monitored with a 107-series temperature probe (Campbell Scientific) in order to register fluctuations and correct for viscosity of the aqueous solution.
The samples were saturated upwards during 24-48 h. The axis translation (Fredlund & Rahardjo, 1993) technique was applied in order to dissolve the remaining air bubbles and to prevent the cavitation of water in the system (Huang et al., 1998). Its principle is based on the upwards translation of the origin of reference for the pore water pressure from the standard atmospheric condition to a larger pressure so as to promote the dissolution of the eventual air phase initially trapped into the sample (Vanapalli et al., 2008). Sample saturation was checked through the calculation of Skempton's coefficient (Fredlund & Rahardjo, 1993), which is defined as the ratio between the increment of pore water pressure and the increment of confinement pressure when the sample is subjected to an increment of confinement pressure (Makhnenko & Labuz, 2013). Once the saturation step was validated, the measurement of the saturated hydraulic conductivity (K s ) was made using the standard constant-head method according to the French norm NF X 30-443 (AFNOR, 2014). Under steady state conditions, the flow rate was monitored until both upstream and downstream flow rates were permanently equal during at least 24 h. Subsequently, the equilibrium condition (zero flux) was imposed by applying a zero pressure gradient through the sample before the unsaturated condition experiment begins. After reaching hydraulic equilibrium, the first pressure increment was applied with the Gas-PVC (Figure 3). The matric head and the cumulative outflow volume were continuously monitored at the base of the ceramic plate. For each pressure step, the net normal stress σ c (kPa) and the matric head h (cm) were calculated by where σ 3 is the cell pressure (kPa), U′ w is the pore water pressure head at the material-plate interface (kPa), and U a is the pore air pressure (kPa). Air pressure was applied through the top cap ( Figure 3) in steps to values of 10, 30, 60, 100, 200, 330, 600, and 900 cm for the calcareous sand interbeds and rock samples. For the powdery limestone samples, the steps were of 20, 50, 100, 200, 400, 600, and 900 cm. The criterion considered to switch to the next pressure step was an outflow value lower than 50 mm 3 h −1 (Eching et al., 1994). The volumetric water content at the end of the outflow experiment was determined by oven drying the sample at 105˚C for 24 h. The saturated volumetric water content (θ s ), in addition to the volumetric water content values at each pressure step, were then back calculated using the oven-dry mass of the sample and the measured cumulative outflow volumes at each pressure step. The water retention curves were then obtained for each sample within the pF value range 0.0-3.0. The measurement of the water retention curve was extended to the pF range 4.0-6.0 by using the WP4C Dewpoint Potentiometer (METER Group). For direct calculation of the unsaturated hydraulic conductivity (K [cm d −1 ]), it was assumed that the water pressure head in the samples varies linearly with depth (Eching et al., 1994). Calculation of K from outflow data was based on the method of Gardner (1956). Additional details can be found in Aldana (2019).

Bulk physical properties measurements
Particle size distribution of unconsolidated materials was determined using the French procedures NF P 94-057 (AFNOR, 1992) for sedimentation analysis and NF P 94-056 (AFNOR, 1996a) for sieving analysis. Specific surface area was measured with the surface analyzer NOVA 2200e (Quantachrome Instruments) via sorption isotherm measurements. Bulk density values (ρ b ) were calculated by dividing the oven dry weight of the samples by their bulk volume. Bulk volumes were calculated from the dimensions of the cylindrical cores.

Mineralogical and geochemical analyses
Mineralogical composition was determined using X-ray diffraction (XRD), both on bulk materials and on their clay fraction. For the bulk mineralogy determination, random powders (<50 μm) were prepared and analyzed using the Debye-Sherrer method (Gravereau, 2012). After carbonate removal, <5-μm clay fraction was characterized using oriented mounts. Diffraction was obtained using the Bragg-Brentano geometry (Gravereau, 2012). Nondestructive and qualitative chemical analysis of the hard limestone rock samples was performed using backscattered electrons detection with an energy dispersive X-ray spectrometer (EDS) integrated to a scanning electron microscope (SEM) (Merlin Compact ZEISS). Textural analysis was performed with an Everhart-Thornley secondary electron detector.
Carbonate content of the samples was determined by calcimetry according to the French procedure NF P 94-048 (AFNOR, 1996b).

Simulation of the water flow in the VZ
The measured hydraulic properties (water retention and hydraulic conductivity data) were used to estimate the hydraulic parameters of each sample using the RETC software . A numerical simulation of the water flow within the VZ of the Beauce limestone aquifer was then undertaken over a period of 52 yr (1 Jan. 1965-31 Dec. 2017) with the HYDRUS-1D software (Šimůnek et al., 2016).

General assumptions
The simulation of water flow in the VZ was performed 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 (Richards, 1931): 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), and K the hydraulic conductivity (cm d −1 ). We used van Genuchten's expression (van Genuchten, 1980) to describe the water retention curve (Equation 4) and the application of the statistical pore connection model established by Mualem (1976) to the van Genuchten model to predict the hydraulic conductivity (Equation 5). with θ r and θ s being the residual and saturated volumetric water content (cm 3 cm −3 ), respectively, α being 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.
with s being the saturated hydraulic conductivity (cm d −1 ), e being the effective saturation (-), and l being a pore connectivity parameter (-), which was fixed at 0.5 (Mualem, 1976).

2.7.2
Representation of the VZ profile Since the majority of the samples (except P B and I B , cf. Table 1) were re-cored from B2, the study of the water flow in the VZ was focalized on this profile, which was reconstituted in HYDRUS-1D. The lithological heterogeneities along this VZ profile were reproduced based on visual descriptions of the undisturbed cored samples. A 23-m-deep profile composed of 13 different materials was created for B2 (Table 1). Each material corresponds to a sample whose hydraulic properties have been determined in the laboratory.

Initial and boundary conditions
As the WTL at the time of the start of the simulation was −19.84 m, the initial matric head profile was defined following Equation 6 with h i = −100 cm for −18.84 < z < 0.00 m and h i varying linearly from −100 to +316 cm from z = −18.84 m to z = −23.00 m.
with 0 being the time of the start of the simulation (d), and ℎ the matric head at 0 at depth z (m). At the soil surface boundary, a water flux was imposed as the upper boundary condition by the daily ETP and rainfall data.
The daily variations of the WTL were used as the lower boundary condition: with ℎ b being the matric head observed at z max (m), and z max = −23.00 m.

Parameterization of the model
The RETC software  was used to fit the θ(h) and K(h) curves (Equations 4 and 5) to the measured water retention and hydraulic conductivity data. The hydraulic parameters of the samples were obtained according to the following steps. For the powdery limestone (P A , P B , and P C ) and the calcareous sand interbeds samples (I A and I B ), initial values of θ r , θ s , α, n and K s were obtained by using the Rosetta software (Schaap et al., 2001) based on the particle size distribution and bulk density measured for each sample. For the hard rock samples, initial values of θ r , α, and n were obtained using the results of a study made near the experimental field site (Amraoui et al., 2017), θ s was obtained using experimental values, and K s was obtained using the first experimental value of the K(h) curves. These initial values were then used as an input for the RETC software. As rec-ommended by Sisson and van Genuchten (1991) and Yates et al. (1992), the relative weights of hydraulic conductivity data against retention data ranged from 0.1 to 1.0. The values optimized with RETC for parameters θ r , θ s , α, n and K s were then used to simulate water flow with HYDRUS-1D from 1 Jan. 1966 to 31 Dec. 2017 (52 yr).

Lithological description and physical properties of the VZ materials
After coring the whole VZ profiles, all materials have been observed and described, from the surface down to a depth of 20 m. We have observed both vertical and horizontal spatial variability (i.e., within a single borehole and between them [B1, B2, and B3]). Four main lithologies were identified along the VZ profiles, including three soft and one hard materials. The first soft material was a brown soil having a fine-grained texture and a plastic consistency ( Figure 4). The first 20 cm corresponded to the ploughed layer except for B1, which was cored in a noncultivated field. The undisturbed soil horizon was found below the ploughed layer and extends from 0.8 (B2) to 1.6 m (B3) deep. The undisturbed soil had a silty loam texture according to USDA (2012) classification with higher clay, silt, and fine sand content and lower proportion of calcite-rich gravel (<3.5%) than the other VZ materials (Table 2). Since compaction caused by hydraulic percussion from coring had modified the soil samples, the results of their hydraulic properties measurements will not be presented.
The second soft material observed in the VZ represents a fragmented powdery limestone. The thickness of this material ranged from 3.4 to 5.9 m. The powdery limestone samples had a loam to silty loam texture according to USDA (2012) classification, with an average clay fraction of 19% (Table 2). This lithology had a crumbly consistency containing cemented aggregates and large-sized particles in the upper part of the interval (between 1 and 6 m deep). The lower end of this layer (between 6 and 7 m deep) had a more silty texture and displayed a slightly plastic consistency ( Figure 5, Table 2).
The third soft material was referred to as calcareous sand interbeds. This material is from 0.5 to 1.5 m thick and was observed between 3.5 and 6.1 m deep. This lithology primarily contains fine to large sand particles and gravels (Table 2), presents a high macroporosity (i.e., clearly visible to the naked eye) (Figure 6), and has a sandy loam to sand texture according to USDA (2012) classification.
The fourth type of lithology observed in the VZ profile was a hard calcite-rich limestone rock having a microcrystalline texture with high carbonate content (>90%) and larger values of bulk density with respect to the soft materials  Note. Particle size distribution and specific surface area analysis were performed only for soft materials since the hard rock materials displayed a well-cemented microcrystalline texture. The analysis of the fine particles (<50 μm) for sample I B was not possible due to insufficient quantity. The value (12.5%) represents the mean bulk fine fraction (clay to fine sand). ρ b , bulk density; SSA, specific surface area; N/A, not available.

Saturated hydraulic properties of the VZ materials
The saturated hydraulic properties of the VZ materials were generally consistent with their bulk physical properties. We observed that the K s values of the powdery limestone and calcareous sand interbeds materials were generally decreasing as particle size decreases whilst the opposite was noted for the values of θ s (Figure 8, Tables 2 and 3). We also observed higher average values of θ s and K s for calcareous sand interbeds than for powdery limestone and rock samples ( Figure 8, Table 3). Differences by a factor of 10 or more were observed for K s between some samples within the same lithology (e.g., P A and P B for powdery limestone and R C and R D for rock samples) (Figure 8). Ould Mohamed (1995) reported average θ s values of 0.35 cm 3 cm −3 for the powdery limestone materials sampled between 0.65 and 1.70 m deep. These values were in good agreement with those observed for P A samples, which was taken within the same depth interval (Tables 2 and 3). Rock samples presented widely spread values of K s with sample R E having the largest among all samples and sample R C having the smallest one (Figure 8). The values of θ s also differ significantly from one sample to another within the same lithology. The θ s values of the rock samples ranged generally from 0.13 to 0.17 cm 3 cm −3 and were relatively low with respect to the soft materials (except for R A , cf. Table 3). The θ s values of the rock samples were in good agreement with those obtained by Legchenko et al. (2020) in a study conducted near the experimental site. An artifact related to the measurement of K s was observed with respect to the values of K(h). The K s values were always lower than the unsaturated values of K close to saturation (Figure 8). This artifact is probably due to the measurement protocol [saturation of the sample, reverse flow direction between K s and K(h) measurements, confinement pressure applied on the sample, cf. Chapuis, 2012] and is currently under investigation. Measured K s values should therefore be considered with caution as they are probably underestimated. This effect was more pronounced for calcareous sand interbeds and rock samples than for powdery limestone samples (Figure 8).

Unsaturated hydraulic properties of the VZ materials
The shapes of the experimental water retention curves of powdery limestone samples were similar despite a small difference between P B and the two other samples (Figure 8b). Powdery limestone samples displayed θ values ranging from 0.25 to 0.36 cm 3 cm −3 for the applied matric potentials 0.5 < pF < 3.0 (Figure 8b). These values were in good agreement with those reported by Michot et al. (2003) for the powdery limestone horizon found just below the soil between 0.88 and 0.98 m deep. In the dry range (4.0 < pF < 6.0), all curves showed similar trends and were in good agreement with their textures ( Table 2). The shape of K(h) curves was quite  different among the samples with a large decrease of K above pF 1.0 for P A and above pF 1.8 for P B . Sample P C showed a more gradual decrease of K in the measurement range when pF increases (Figure 8b). The average water retention capacity of calcareous sand interbeds was the lowest among the soft materials, except at saturation (Figure 8c). The K(h) values for calcareous sand interbeds samples were close throughout the full range of pF.
The K values at pF 1.0 were higher than those of the other VZ materials and then decrease drastically when pF increases (Figure 8c).
The hydraulic properties of the rock samples were highly heterogeneous (Figure 8d). Sample R A showed a fast decrease of θ between saturation and pF 1.0, probably induced by its secondary porosity. Then, water retention capacity decreased gradually between pF 1.0 and 5.0 (Figure 8d). Samples R B and R C had a nearly constant water content in the range between saturation and pF 2.5 where no or little drainage occurred which means that they hardly desaturated. On the contrary, water retention of samples R D and R E decreased gradually over the same range of pF (Figure 8d). The water content for most of the rock samples at the end of the unsaturated experi-ment (for pF around 3.0) was close to 0.06 cm 3 cm −3 . Regarding K(h), all samples displayed a progressive decrease of K as pF increases except for sample R C , whose K values remained nearly constant through the 0.0-3.0 pF range (Figure 8d).
The θ r , θ s , α, n, and s parameters fitted with RETC based on the experimental results provided a satisfactory description of the experimental water retention and hydraulic conductivity curves (R 2 > .915, except for P B and R C ) (Figure 8, Table 4). The trends in the values of the parameters θ s and s fitted with RETC for the VZ materials were consistent with the experimental observations made above, although the water retention curves for samples I A , I B , and R A could have been more accurately reproduced near saturation by using a dual-porosity model (Figure 8, Tables 3 and 4).

Relationships between the hydraulic properties and the physical, mineralogical, and geochemical characteristics of the VZ materials
In this section, we explore whether the hydraulic properties of the VZ materials could be related to their physical, T A B L E 4 Residual (θ r ) and saturated (θ s ) water content, fitting parameters (α and n), and saturated hydraulic conductivity (K s ) fitted using the RETC software for the vadose zone materials mineralogical, and geochemical characteristics. Because distinct characterization methods were applied to soft and hard samples, these two types of materials will be discussed separately.

Soft materials (soil, powdery limestone, and calcareous sand interbeds samples)
In addition to calcite, we clearly observed that palygorskite (a nonswelling clay) and other secondary clay minerals were present in the soft materials (Figure 9a and 9b). As shown by glycolation, smectite (a swelling clay) was present in the soil sample in larger proportions with respect to the powdery limestone materials (Figure 9b). The presence of smectite and palygorskite is supported by the specific surface area analysis that showed that the largest values correspond to soil (Table 2), as also observed by Ould Mohamed et al. (1997). The soil samples showed low carbonate content (<18%, cf. Table 2) due to weathering (Michot et al., 2003). The presence of kaolinite and illite (nonswelling clays) (Figure 9b) is a sign of the more intense weathering of soil (Velde & Meunier, 2008) compared with the other VZ materials. When linked to the hydraulic properties of the soil samples given by Ould Mohamed et al. (1997) and the information given by Michot et al. (2003), these physical and mineralogical results could explain the highest water retention capacity observed for the soil in comparison with the other VZ materials (Figure 8).
The lower water retention capacity of powdery limestone compared with the soil samples (Figures 8a and 8b) was related to the lower clay fraction and the lower content in smectite of this material (Table 2, Figure 9b). The hydraulic conductivity values of the powdery limestone samples were relatively low given their loam to silty loam texture according to USDA (2012) classification. The fact that the samples P A and P B were less silty than sample P C (Table 2) could explain their lower K(h) values for pF > 2 (Figure 8b).
Regarding the hydraulic properties of the calcareous sand interbeds samples, θ(h) curves displayed a rapid decrease within the range of saturation to pF 1.0 and mean K values near saturation were the largest among all soft materials (Figure 8c). These results were expected since calcareous sand interbeds had large particle sizes (Table 2). We further observed a larger decrease of the θ values near saturation (pF < 1.0) for I B , which was in good agreement with its texture since the latter displayed a coarse-grained texture while I A had a finer-grained texture (Figure 8c, Table 2).

Hard materials (limestone rock samples)
Compared with the soft materials (Figures 8a, 8b, and 8c), the limestone rock samples displayed a higher heterogeneity in their hydraulic properties (Figure 8d). This variability was related to the natural composition of the rock and to the geometry of its pores. The relatively low experimental θ s values found for the rock samples (except R A , cf. Table 3) are due to the diagenetic processes of compaction and cementation which yielded the rock. As long as the rock has no secondary porosity developed by weathering processes, its water retention and hydraulic conductivity are expected to be very low (Tiab & Donaldson, 2004). Samples R B and R C showed strong differences in their hydraulic properties compared with samples R A , R D , and R E (Figure 8d). Observed through the naked eyed, R B and R C samples seemed hardly permeable at first sight ( Figure 10a) and could therefore be considered as a massive rock which is mainly composed of calcite (Figure 11a). However, the SEM/EDS analysis highlighted the presence of few microfissures with an average thickness of 15 μm (Figure 10b). We also observed that calcite was replaced by a 2:1 Al-rich phyllosilicates (palygorskite) (Figures 10b and 11b) in some areas where water was probably stored as a result of adsorption and capillary forces (Tiab & Donaldson, 2004). These types of phyllosilicates are typically observed in partially saturated pores where capillary water exists and solute transport is controlled by chemical diffusion (Meunier, 2005). The presence of swelling clays such as smectites (Figure 11b) and the absence of large pores (>15 μm) could explain the nearly constant θ values measured between saturation and pF 2.5 for samples R B and R C and why no desaturation was observed (Figure 8d). The presence of clays in the rock matrix also showed that the massive rock was permeable even if no F I G U R E 9 Diffractograms obtained by X-ray diffraction (XRD) for (a) the bulk material (<50 μm, random powder) and (b) the clay fraction after decarbonatation (<5 μm, oriented mounts) of one soil (S B ) and one powdery limestone (P B ) sample. EG, ethylene-glycol treatment F I G U R E 1 1 Diffractograms obtained by XRD for (a) the bulk material (<50 μm, random powder) and (b) the clay fraction after decarbonatation (<5 μm, oriented mounts) of a massive rock (R C ), an upstream (R A ) and a downstream (R D ) weathered rock sample. EG, ethylene-glycol treatment macroporosity (fractures, vugs) was observed. As a consequence of physical and/or chemical reactions, the migration of fine particles (clay minerals), notably illite, smectite, and kaolinite may also have caused a reduction of the permeability of the rock materials induced by the accumulation of those minerals in pore throats (Aksu et al., 2015;Russell et al., 2017;Wilson et al., 2014).
Observed through the naked eye, R A , R D , and R E samples showed intense weathering highlighted by the presence of fissures (<1 mm) ( Figure 12a) and vugs (<5 mm) (Figure 13a) at their surfaces and by the replacement of calcite by secondary minerals (Figures 11b, 12b, and 13c). These weathered rock samples have developed a secondary porosity and displayed higher K values near saturation (pF range 0.5-1.5) and lower (or close) K values within the pF value range 2.0-3.0 than the massive (microfissured) samples (R B and R C ) (Figure 8d).
The sample R A was localized in the VZ profile within the sequence of powdery limestone and calcareous sand interbeds materials, and therefore also above the massive rock facies (Table 2), and could consequently be considered as an example of upstream weathered rock. We observed that the calcite (represented in blue in Figure 12b) was replaced by phyllosilicates (represented in red-brown in Figure 12b) in the matrix of the rock which highlighted that water has been flowing through the sample. According to the diffractograms, palygorskite is the main secondary mineral in the rock matrix,  (Figure 11a and 11b). The value of ρ b (Table 2) and the presence of palygorskite (Figure 11b) were consistent with the high experimental value of θ s observed for R A ( Table 3). The near total absence of smectite (Figure 11b) could explain the large decrease in water retention capacity of the R A sample within the pF value range of saturation to 1.0 (Figure 8d).
Within the VZ profile, the samples R D and R E were localized just above or in the zone of water table fluctuation (Table 2) and could therefore be considered as examples of downstream weathered rock. We observed voids spread all over the surface of the two rock samples (Figure 13a), and a chemical mapping of these vugs for the sample R D showed that the material has undergone dissolution-recrystallization F I G U R E 1 4 Analysis of a natural macrofracture and associated secondary minerals observed in a rock sample localized at a depth of 6.6 m in B2: (a) highlight of the macrofracture, and (b) chemical analysis of the macrofracture profile using scanning electron microscopy (SEM, energy dispersive X-ray spectroscopy [EDS] detector); blue: Ca, red: Fe episodes, which produced a secondary porosity (Putnis, 2002) and yielded the rock into a more permeable material through weathering processes (Velde & Meunier, 2008) (Figures 13b  and 13c). Some authors indeed observed that weathering has a major impact in enhancing permeability in Cenozoic limestone aquifers, notably through the enlargement of interconnected vugs (Cunningham et al., 2009;Worthington et al., 2016). The red round-shape ring showed the precipitation of goethite (Figure 13c), which is probably related to an episode of drying of the aquifer. Calcite has also recrystallized in the vugs (dissolution voids) over the oxide ring (Figure 13c), which highlighted that water has been flowing through this network after oxide deposition. As illustrated in Figure 11b and Figure 13c, the matrix of the downstream weathered rock and the dissolution voids also showed the presence of phyllosilicates (swelling and nonswelling clays), which replaced calcite. The presence of kaolinite and smectite, along with iron oxides (goethite), highlights recent water-rock interactions (Velde & Meunier, 2008). The lower water retention capacity of the upstream (R A ) and downstream weathered rock samples (R D and R E ) compared with the massive rock samples (R B and R C ) (Figure 8d) could also be explained by their smaller proportions in smectites (Figure 11b). Natural (macro-) fractures (>>1 mm) were observed during the visual description of core samples (Figure 14a). At the core scale, evidence was found that these fractures had been exposed to water flow as demonstrated by the presence of iron oxides (goethite) next to the macrofracture. Only calcite and goethite were observed next to the fractures, and no replacement of calcite by phyllosilicates was observed (Figure 14b). This indicates that intensive leaching of water has occurred in the open fractures (Meunier, 2005). These fractures have been active at some stage in the history of the Beauce limestone aquifer, leading to much larger values of K than those measured on the samples at the laboratory scale. However, this phenomenon could occur only under or close to fully saturated conditions.

Application to the simulation of water flow through the VZ profile
The simulated θ profile for B2 at the date of the drilling (14 Mar. 2017) was in good agreement with the experimental measurements made on the undisturbed core samples at the same date (R 2 = .900, cf. Figure 15a). Overestimations of experimental θ by the model were noted in P A material (2 m deep) and, to a lesser extent, in R D material from 11.0 to 15.0 m deep, whereas underestimations of experimental θ by the model were noted in the R B material (8.5 m deep) ( Figure 15a, Table 1).
Simulated θ profiles observed at the minimum WTL (−22.45 m on 26 Aug. 1992), corresponding to a dry VZ profile, and at the maximum WTL (−14.84 m on 18 May 2001), corresponding to a wet VZ profile, showed that the most significant changes in θ (>0.10 cm 3 cm −3 ) were observed in the powdery limestone P A (between 1.0 and 3.5 m deep) and calcareous sand interbeds (I A and I B , between 3.5 and 5.5 m deep) materials with differences up to 0.23 cm 3 cm −3 between the two profiles ( Figure 15b, Table 1). Inversely, slight variations of θ (<0.05 cm 3 cm −3 ) were found in P C and in all rock materials until the water table was reached and could be explained by the nearly constant θ values measured between saturation and pF 3.0 in these materials (Figures 8b and 8d).
Positive values of matric head were noted in several VZ materials located above the least permeable materials (P C and R C , cf. Table 4) ( Figure 16) and were correlated with years showing annual rainfall-ETP balance close or above 0 cm (1978-1984, 1988, 2001-2002). This indicates that perched water tables may have occurred in the VZ profiles, caused by accumulation of water above the weakly permeable materials. F I G U R E 1 6 Simulated matric head observed at the observation nodes for the vadose zone B2 profile rock material R C taken into account in the representation of the VZ profile for B2 (0.2 m, cf. Table 1), its very low K values, which remained nearly constant over the range of saturation to pF 2.5 (Figure 8d), caused considerable fluctuations of the matric head and water status in the materials situated above (cf. R A , I B , P C , and R B in Figure 16). The materials located below R C showed less variations of their matric head and a drier water status over the simulation period (cf. R D in Figure 16). These results may suggest that the water flow in the downstream weathered rock materials (R D and R E ) is generally slow, and equivalent to those in the massive rock (R C ) unless the water table reaches up these materials and brings them to saturation. However, given the high lithological heterogeneity of the rock materials (Figure 7) and the presence of (macro-)fractures (>1 mm) along the VZ profile (Figure 14), these observations and the occurrence and thickness of the perched water tables would probably be less significant at the three-dimensional field scale, as the water could and probably will be drained laterally by preferential flows through fractures into the more permeable materials.

Summary of the results
We have explored the relationship between hydraulic properties and physical, mineralogical, and geochemical features within a heterogeneous VZ of a limestone aquifer. Three boreholes were drilled in March 2017, and four types of lithology were identified along the VZ profiles, including three soft (soil, powdery limestone, and calcareous sand interbeds) and one hard (limestone rock) materials. Based on all the results obtained through this study, a conceptual scheme summarizing the properties of the VZ materials obtained at the laboratory scale and the impact of these properties on water flow at the borehole scale has been produced (Figure 17). At the top of the VZ, the soil presented the highest water retention capacity, which could be explained by the presence of swelling clay (smectite) and palygorskite, and high variations of water content and matric head, which indicated that water flow was strongly influenced by meteorological patterns (rainfall and ETP).
The fragmented powdery limestone observed beneath the soil displayed low hydraulic conductivity, which suggests that slow water flow takes place along these materials. The lower part showed smaller variations of water content than the upper part, which could be explained by its higher water retention capacity.
The calcareous sand interbeds showed significant changes in their water content due to their high saturated water content and low water retention capacity. The high hydraulic conductivity of this material near saturation is responsible for fast water flow for a wet water status, whereas its low hydraulic conductivity when drying results in slow water flow for a dry water status.
The properties of the hard limestone rock materials were highly heterogeneous. The upstream weathered rock showed fissures (<1 mm) and the development of a secondary porosity resulting from the replacement of calcite by phyllosilicates (mainly palygorskite) in its matrix. The absence of smectite could explained the low water retention capacity of this rock facies. The hydraulic properties of the upstream weathered rock indicated that the water flow can be fast near saturation and slow under dry conditions. The massive rock facies, which seemed hardly permeable at first sight, actually showed the presence of microfissures (15 μm thick) around which calcite was replaced by secondary clay minerals with the presence of (a) palygorskite, which could suggest that water was stored by capillary forces and adsorption processes and that solute transport through these microfissures was driven by chemical diffusion; (b) smectite, which could indicate that this material was permeable even if macropores were not observed and could explain its water retention capacity. The very low hydraulic conductivity of this material may cause accumulation of water in the facies located above and thus lead to the occurrence of perched water tables in the VZ profile. The downstream weathered rock, which was localized just above or in the zone of water table fluctuation, showed vugs (<5 mm) caused by dissolution processes, with precipitation of iron oxides (goethite) in it and replacement of calcite by phyllosilicates which emphasized water-rock interactions. This facies has also undergone dissolution-recrystallization episodes, which could have enhanced its permeability through the enlargement of vugs. The downstream weathered rock located in the VZ profile showed slight variations of water content and an overall dry water status as long as the water table remains far below. In view of its hydraulic properties measured at the laboratory (sample) scale, the water flow along this facies is highly affected by the presence of the weakly permeable massive rock located above, which could restrict downward water flow and promote perched water tables and lateral flow towards large fractures.

CONCLUSION
Valuable information about the link between hydraulic properties and microscopic and macroscopic features of the geological formations of the Beauce limestone aquifer were gathered throughout this multidisciplinary study. Insufficiently studied so far, the hydraulic properties of the lithologies encountered within the VZ of this aquifer showed strong contrast and stressed the high heterogeneity of the water retention and hydraulic conductivity properties of both soft and hard materials. The petrographic investigations implemented through the mineralogical and geochemical analysis provided key evidence on the influence of the presence or absence of different types of secondary minerals and the impact of weathering processes on water flow within the VZ. Finally, the soft and hard VZ materials displayed highly heterogeneous hydraulic properties and mineralogical composition, which seems to have a considerable influence on water retention and water flow at the borehole scale. It also seems that the fissures and vugs developed by weathering processes represent a more or less interconnected network of pores through which water can flow, but which may drain quite rapidly when drying. The weakly permeable massive rock facies, which displayed higher water retention capacity and lower hydraulic conductivity than the weathered rock, could lead to the occurrence of perched water tables within the VZ and seems to affect significantly water flow through the materials situated above and below. After this study, ongoing work is focused on describing the impact of the materials heterogeneity and climate on water flow and solute travel time through this highly heterogeneous VZ.
The deciphering of the multitude of coupled processes that govern water flow and solute transport within such a heterogeneous VZ is complex and cannot be performed in an exhaustive way with the examination of only a dozen samples at the laboratory scale. It is also worth noting that the presence of natural fractures, found at the field scale and that had been subjected to intense leaching of water at some point in the past, may lead to preferential water flow resulting in much higher values of hydraulic conductivity than those measured on the samples at the laboratory scale. These observations underline the necessity to carry out further comprehensive in situ studies through the instrumentation of the VZ materials with innovative hydrogeological and geophysical sensors, which will allow data acquisition over a long time and help simulating one-or three-dimensional water flow and solute transport at larger scales. This will allow a better understanding of transport processes in the VZ of the vulnerable Beauce limestone aquifer and help sustainable management and preservation of the water resource.

A C K N O W L E D G M E N T S
This study was conducted within the framework of the O-ZNS project which is part of PIVOTS project (https://plateformespivots.eu/o-zns/?lang=en). We gratefully acknowledge the financial support provided to the PIVOTS project by the Région Centre-Val de Loire (ARD 2020 program and CPER 2015−2020) and the French Ministry of Higher Education and Research (CPER 2015−2020 and public service subsidy to BRGM). This operation is co-funded by European Union. Europe is committed to the Centre-Val de Loire region with the European Regional Development Fund. This research work was co-funded by the Labex VOLTAIRE (ANR-10-LABX-100-01). We also gratefully acknowledge Stephane Ruy (UMR 1114 EMMAH INRA-UAPV) for having allowed us to access his laboratory and use the WP4C Dewpoint Potentiometer.