Nolwenn Lesparre
email: lesparre@unistra.fr

Sylvain Pasquet
email: sylvain.pasquet@sorbonne-universite.fr

Philippe Ackerer
email: ackerer@unistra.fr

Impacts of hydrofacies geometry designed from seismic refraction tomography on estimated hydrogeophysical variables

Understanding the critical zone processes related to groundwater flows relies on subsurface structure knowledge and its associated parameters. We propose a methodology to draw the patterns of the subsurface critical zone at the catchment scale from seismic refraction data and show its interest for hydrological modelling. The designed patterns define the structure for a physically based distributed hydrological model applied to a mountainous catchment. In that goal, we acquired 10 seismic profiles covering the different geomorphology zones of the studied catchment. We develop a methodology to analyze the geostatistical characteristics of the seismic data and interpolate them over the whole catchment. The applied geostatistical model considers the scale variability of the subsurface structures observed from the seismic data analysis. We use compressional seismic wave velocity thresholds to identify the depth of the soil and saprolite bottom boundaries. Assuming that such porous compartments host the main part of the active aquifer, their patterns are embedded in a distributed hydrological model. We examine the sensitivity of classical hydrological data (piezometric heads) and geophysical data (magnetic resonance soundings) to the applied velocity thresholds used to define the soil and saprolite boundaries. Different sets of hydrogeological parameters are used in order to distinguish general trends or specificities related to the choice of the parameter values. The application of the methodology to an actual catchment illustrates the interest of seismic refraction to constrain the structure of the critical zone subsurface compartments. The sensitivity tests highlight the complementarity of the analyzed hydrogeophysical data sets.

Introduction

Groundwater flow and catchment discharge are strongly controlled by the structure of the critical zone (CZ) underground part and its related hydraulic properties and boundary conditions [START_REF] Cassidy | Combining multi-scale geophysical techniques for robust hydro-structural characterisation in catchments underlain by hard rock in post-glacial regions[END_REF][START_REF] Diek | The effect of spatial soil variation on the hydrology of a semi-arid Rocky Mountains catchment[END_REF][START_REF] Fleckenstein | River-Aquifer Interactions, Geologic Heterogeneity, and Low-Flow Management[END_REF][START_REF] Gabrielli | The role of bedrock groundwater in rainfall-runoff response at hillslope and catchment scales[END_REF]. The bottom limit of the CZ corresponds to the base of the aquifer above which alteration of underground materials typically increases towards the surface [START_REF] Anderson | Physical and Chemical Controls on the Critical Zone[END_REF][START_REF] Brantley | Crossing disciplines and scales to understand the critical zone[END_REF]. From the substratum to the soil surface, the deeper weathered rocks progressively evolve to saprolite and soil, designing the compartments of the underground CZ classical scheme [START_REF] Anderson | Physical and Chemical Controls on the Critical Zone[END_REF]. The porosity and hydraulic conductivity of such compartments increase toward the surface as fractures, weathering, and alteration processes enlarge the porous space and ease water flows [START_REF] Brooks | Hydrological partitioning in the critical zone: Recent advances and opportunities for developing transferable understanding of water cycle dynamics: Critical Zone Hydrology[END_REF]. In the following, we refer to the soil as the superficial compartment made of disaggregated materials, while the saprolite denotes deeper weathered materials with a lower porosity. In mountainous environments, the subsurface is particularly heterogeneous, and the weathered bedrock, saprolite, and soil compartments can show significant thickness variations at the catchment scale [START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF][START_REF] Diek | The effect of spatial soil variation on the hydrology of a semi-arid Rocky Mountains catchment[END_REF][START_REF] Koch | Joint interpretation of hydrological and geophysical data: electrical resistivity tomography results from a process hydrological research site in the Black Forest Mountains, Germany[END_REF]St. Clair et al., 2015). The thickness variability of these CZ underground compartments is related to the history of climate, weathering and erosion processes, regional tectonic forcing, or the occurrence of some metamorphic intrusions. These processes might have different impacts across the catchment due to local topography and lithology [START_REF] Anderson | Physical and Chemical Controls on the Critical Zone[END_REF][START_REF] Holbrook | Links between physical and chemical weathering inferred from a 65-m-deep borehole through Earth's critical zone[END_REF][START_REF] Rempe | A bottom-up control on fresh-bedrock topography under landscapes[END_REF][START_REF] Riebe | Controls on deep critical zone architecture: a historical review and four testable hypotheses: Four Testable Hypotheses about the Deep Critical Zone[END_REF].

The spatial distribution of the subsurface hydraulic properties determines the way infiltrated water drains into storage areas [START_REF] Brooks | Hydrological partitioning in the critical zone: Recent advances and opportunities for developing transferable understanding of water cycle dynamics: Critical Zone Hydrology[END_REF]. It has been shown that the thickness distribution of the underground CZ compartments impacts the watershed's water budget [START_REF] Bertoldi | Impact of Watershed Geomorphic Characteristics on the Energy and Water Budgets[END_REF][START_REF] Lanni | Simulated effect of soil depth and bedrock topography on near-surface hydrologic response and slope stability[END_REF]. Moreover, the hydrogeological facies geometry is a key property for understanding the dynamic of the groundwater from piezometric measurements (piezometers may intercept different layers). Knowing the hydrogeological facies geometry is thus crucial to design hydrogeological models, especially under phreatic conditions [START_REF] Carrera | Inverse problem in hydrogeology[END_REF]. However, the inverse problems that seek the hydraulic parameters of distributed hydrologic models applied at the catchment scale are known to show a strong non-uniqueness [START_REF] Ebel | Physics-based hydrologic-response simulation: Seeing through the fog of equifinality[END_REF]. The value of the hydraulic parameters related to each mesh element of such models has to be determined, while the available measurements might be equally fitted by different sets of properties. To tackle this non-uniqueness issue, spatialized observations providing information on diverse hydraulic properties should be integrated in the inversion process (Zhou et al. 2014). Nevertheless, hidden by nature, the measurement of water storage and flow properties in the underground complex structure is still arduous. Basic but crucial information, such as the interface geometry between the different CZ underground compartments, are often missing [START_REF] Brooks | Hydrological partitioning in the critical zone: Recent advances and opportunities for developing transferable understanding of water cycle dynamics: Critical Zone Hydrology[END_REF]. Recent studies show that the characterization of the underground CZ structure remains challenging [START_REF] Flinchum | Critical Zone structure under a granite ridge inferred from drilling and three-dimensional seismic refraction data[END_REF][START_REF] Gourdol | Exploring the regolith with electrical resistivity tomography in large scale surveys: electrode spacing related issues and possibility[END_REF][START_REF] Kaufmann | Simultaneous multichannel multi-offset ground-penetrating radar measurements for soil characterization[END_REF][START_REF] Pasquet | Catchment-Scale Architecture of the Deep Critical Zone Revealed by Seismic Imaging[END_REF].

Geophysical imaging methods provide an insight into the underground CZ architecture as they furnish a vision of the subsurface geophysical properties with a continuous spatial coverage along acquisition profiles. In particular, seismic refraction tomography (SRT) supplies structural information of the CZ underground part. SRT highlights the spatial variability of the subsurface properties and can be used to distinguish characteristic patterns [START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF][START_REF] Cassidy | Combining multi-scale geophysical techniques for robust hydro-structural characterisation in catchments underlain by hard rock in post-glacial regions[END_REF][START_REF] Dal Bo | Geophysical imaging of regolith in landscapes along a climate and vegetation gradient in the Chilean coastal cordillera[END_REF][START_REF] Olona | Weathering variations in a granitic massif and related geotechnical properties through seismic and electrical resistivity methods[END_REF][START_REF] Huang | Bayesian seismic refraction inversion for critical zone science and near-surface applications[END_REF]. The inversion of SRT data provides a distribution of compressional P waves velocity ( p v ) in the subsurface, which depends mainly on the medium's porosity, density and mineralogy. Weathering processes occurring in the CZ induce a decrease of p v by increasing the degree of fracturation and porosity. Indeed, p v is slower in pores filled by air or water than in the rock matrix (Pasquet et al., 2015;[START_REF] Parsekian | Multiscale geophysical imaging of the critical zone: Geophysical Imaging of the Critical Zone[END_REF]. Moreover, p v is lower in secondary minerals (i.e., clays, oxides) than in parent minerals (i.e., quartz, plagioclase) [START_REF] Olona | Weathering variations in a granitic massif and related geotechnical properties through seismic and electrical resistivity methods[END_REF][START_REF] Parsekian | Multiscale geophysical imaging of the critical zone: Geophysical Imaging of the Critical Zone[END_REF]. SRT is thus well suited to distinguish the spatial variability of the interfaces between the CZ underground compartments [START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF][START_REF] Flinchum | Critical Zone structure under a granite ridge inferred from drilling and three-dimensional seismic refraction data[END_REF][START_REF] Holbrook | Geophysical constraints on deep weathering and water storage potential in the Southern Sierra Critical Zone Observatory[END_REF][START_REF] Olona | Weathering variations in a granitic massif and related geotechnical properties through seismic and electrical resistivity methods[END_REF][START_REF] Olyphant | Topographic correlations with soil and regolith thickness from shallow-seismic refraction constraints across upland hillslopes in the Valles Caldera, New Mexico[END_REF][START_REF] Huang | Bayesian seismic refraction inversion for critical zone science and near-surface applications[END_REF][START_REF] Pasquet | Catchment-Scale Architecture of the Deep Critical Zone Revealed by Seismic Imaging[END_REF]St. Clair et al., 2015).

In this study, we assume that the subsurface description obtained by SRT holds also for hydraulic properties, i.e. the layers distinguished with SRT present a homogeneous porosity or hydraulic conductivity at the catchment scale. From this assumption, we assess the impact of the underground geometry based on SRT on variables dependent on the groundwater storage estimated from hydrological modelling. We use parameter values obtained from previous studies in the Strengbach catchment [START_REF] Belfort | Vadose Zone Modeling in a Small Forested Catchment: Impact of Water Pressure Head Sampling Frequency on 1D-Model Calibration[END_REF]Lesparre et al., 2020) to perform this analysis under field site conditions. The following methodology is applied:

1.

Measured SRT profiles are analysed in a geostatistical framework using filtering, truncated power value variograms due to the change in measurement scale with depth, and 250 three dimensional (3D) velocity fields are generated over the catchment;

2.

A threshold velocity is prescribed to estimate the layers' thicknesses for the 250 velocity fields;

3. Numerical simulations are performed for the 250 geometries using a hydrological physically-based model of the water catchment and uniform hydraulic parameters for each layer. Model outputs of interest are piezometric heads, i.e. levels of the saturated zone as the groundwater is unconfined and Magnetic Resonance Soundings (MRS), which include information in the water content of the unsaturated zone.

4.

The impact of the geometry and additional uncertainties on threshold velocity values and hydraulic parameters is evaluated by a simplified sensitivity analysis of MRS and piezometric heads. We focused the analysis on spatialized data as they better show how the data sensitivity to the tested conditions depends on the local context (i.e. steep slopes leading to drainage or flat regions favouring storage).

The originality of this paper lies in the framework developed here to build the hydrological model geometry from SRT data and assess how variables informing on the groundwater state estimated from that model are sensitive to the geometry uncertainty. The field site and the treatment of the SRT data are described in section 2 and 3, respectively. The construction of the 3D geometries from the SRT data and the geostatistical analysis are explained and discussed in section 4. Then, the hydrological model NIHM and its related outputs are exposed in section 5. Finally, we evaluate how the simulated water circulations are impacted by the model interface geometry through an analysis of measurable data sensitivity in section 6.

Field context

Studied site

The Strengbach watershed is located in the Vosges Mountains (Northeast France) and covers an area of about 0.8 km² (Fig. 1).

The elevation ranges between 883 and 1146 m, and the topography is rugged with incised slopes that can reach up to 30°. The catchment is divided into two hillsides with different morphology and meteorology influenced by their respective orientation.

The southeast slopes are gentler; the temperature is usually lower and associated with higher precipitation than the northwestern slopes.

The subsurface can be described using the classical CZ scheme with a degree of weathering and fracturation that increases from depth towards the surface [START_REF] Brantley | Designing a network of critical zone observatories to explore the living skin of the terrestrial Earth[END_REF]. Most of the catchment lies on a Hercynian granitic bedrock, but microgranite and gneiss constitute the protolith of the southern and northern crests, respectively [START_REF] Gh'mari | Etude minéralogique, pétrophysique et géochimique de la dynamique d'altération d'un granite soumis aux dépôts atmosphériques acides (bassin versant du Strengbach, Vosges, France) : Mécanismes, bilans et modélisations[END_REF]. That hard-rock level may be locally fractured and is overlaid by weathered bedrock made of chemically altered and fractured rocks. When sufficiently altered, this weathered bedrock turns into saprolite, forming a sandy coarse-grain matrix containing gravels and pebbles (Fichter et al., 1998a). Then, the soil composes the uppermost layer. The physical properties of each of these two layers and their respective thicknesses vary spatially over the catchment, as expected by El Gh' Mari (1995) and confirmed by a recent hydrogeophysical study (Lesparre et al., 2020a). Some catchment regions, such as crests, slopes and valley bottom might present distinct soil and saprolite thicknesses and varying porosity distributions. The exposure and the inclination could also impact thickness and porosity distributions, allowing observations of differences from one hillside to another. The weathering history of the hillsides also differs as hydrothermal circulations altered the northern slope 180 My ago (Fichter et al., 1998b).

Meteorological and hydrological observations

The Strengbach catchment is a well-studied research site that hosts numerous scientific investigations spanning various key questions concerning the functioning and the vulnerability of the CZ [START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF]. Permanent measurement stations have continuously monitored the meteorological and environmental conditions since 1986. These long-term observations are managed by the Observatoire Hydro-Géochimique de l'Environnement (OHGE, http://ohge.unistra.fr; CNRS/University of Strasbourg), which is part of the French network of CZ observatories (OZCAR; [START_REF] Gaillardet | OZCAR: The French network of critical zone observatories[END_REF]. OHGE also provides convenient facilities for punctual scientific experiments [START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF].

The meteorological forcing is monitored at two stations: one placed on the northern crest of the catchment and the other settled near the outlet (Fig. 1). Both stations record rainfall, temperature and relative humidity. The upper station also monitors global radiation, wind speed and snow thickness. Seven rain gauges provide regular measurements covering the catchment to infer rainfall spatial variability.

The stream flow rate is continuously monitored at the catchment outlet with an H-flume (RS station). A second flume records the flow rate upstream (RAZS station). The underground structure of the catchment was investigated by drilling nine boreholes with depths from 15 to 120 m. Three boreholes were cored to provide direct insight into the deep CZ structure. Most of the boreholes intercept fractures in the bedrock, such that monitored water levels do not necessarily display hydraulic pressure heads corresponding to the water flow dynamics in the shallow porous medium. Recently, 10 piezometers were drilled to depths lower than 7 m and have recorded data since September 2020. Those piezometers are spatially distributed over the catchment with a few of them installed near the original boreholes.

Hydrogeological knowledge

A combined analysis of the catchment pedology and MRS data covering the catchment has shown a relatively flat region of water storage upstream the creek main spring [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF]Lesparre et al., 2020a). A previous analysis of MRS measurements rendered a qualitative depiction of the subsurface water volume distribution in the catchment [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF][START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF]. The map of the water volume concealed in the weathered layer shows significant variability strongly correlated with the pedologic zonation. Low water contents are suggested on the northern crest by MRS measurements, but clayey materials covering that region might prohibit the detection of subsurface water [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF]. On the northern hillside, the shallow subsurface is made of fissured/fractured granites that intense hydrothermal circulations have altered in the past. There, the estimated water volume is intermediate and seems to feed perennial flow over time [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF]. The southern hillside, which appears to be less weathered (Fichter et al., 1998b), shows lower water content with a drier vadose zone less prone to infiltration or better drained than the northern hillside.

Finally, water content is higher underneath the wetland in the downstream part of the catchment and under the flat colluvium zone, which most likely corresponds to the thickest porous subarea of the catchment [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF]. MRS signals recorded in that zone (called zone 2 in Lesparre et al., 2020a, see Fig. 1) show higher amplitudes than other acquisition locations, indicating a higher water content at depth. This can be explained by a thicker water bearing unit in that zone, as suggested by the results of the hydrological modeling (Lesparre et al., 2020a).

3 Seismic refraction data

Seismic velocity profiles

Ten SRT profiles, covering a total length of 2 km, were acquired in June 2018 and August 2019. Their locations were chosen to cover specific areas of the catchment, such as the valley bottom, the crests, the region upstream of the creek spring and both hillsides (Fig. 1). The surveys were designed to explore how the underground part of the CZ evolves in these different regions, which were previously distinguished by a joint analysis of pedological and MRS data collected across the catchment [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF]Lesparre et al., 2020a). Seismic data were collected along profiles of different lengths with 24-channel seismic recorders (Geometrics) and 14-Hz vertical-component geophones. The inter-distance between geophones was fixed to 2 m, and the sources were distant by 8-10 m. The source signal was generated using a 5 kg sledgehammer swung on a metal plate.

For each shot, the seismic wave propagation was recorded with 72-144 geophones, depending on the profile as summarized in Table S1.

First arrival times were picked manually on each trace gathered by recorded shots, when the signal-to-noise ratio is high enough to confidently identify first breaks. The observed travel times were then used to build the subsurface P-wave velocity structure ( p v ) by solving an inverse problem with the pyGIMLi refraction tomography inversion module (Rücker et al., 2017). In pyGIMLi, the inversion domain corresponds to a triangular mesh with cells of constant velocity through which rays are traced using a shortest-path algorithm [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF][START_REF] Moser | Shortest path calculation of seismic rays[END_REF]. The velocity in each mesh cell is estimated using a generalized Gauss-Newton inversion framework. The inversion is iterative and starts with an initial model consisting of a velocity field that increases linearly with depth from [250 -750] m/s at surface to [2000 -5000] m/s in depth (Table S2). The velocity field is then smoothly updated at each iteration in order to reach the closest match between predicted and observed travel times.

Inversions were performed with 144 combinations of starting models and regularization parameters (Table S2) in order to explore the possible solutions and estimate the uncertainty of the velocity distribution along each profile [START_REF] Pasquet | Geophysical imaging of shallow degassing in a Yellowstone hydrothermal system[END_REF].

A selection is then applied to keep only the results of inversions performed with a set of parameters that obtained a root mean square error < 2.5 ms and a root mean square error weighted by the variance  the travel time measurement error. We applied a systematic error of 5% on each picked travel time, setting obs  lower and upper bounds at 0.3 and 3 ms, respectively. Among the 144 combinations of starting models and regularization parameters, 104 inversion results fulfill these requirements for all 10 profiles. The mean and the standard deviation of p v are then computed for each pixel of the SRT profiles from the 104 selected models (Fig. S1 andS2). The standard deviation distribution provides an estimate of the p v variations' likelihood.

Each seismic profile inversion result is extracted to build horizontal maps of the average p v distributions at different depths (Fig. 2). Each profile was flattened, so the depth of each point corresponds to its orthogonal distance to the surface. The standard representations of the seismic profiles (distance vs. elevation) are also given in Fig. S1. p v varies globally between 400 m/s and 4500 m/s and increases progressively downward. Above a depth of 3 m, p v values are globally homogeneous and remain below 700 m/s (Fig. 2). At a depth of 3 m, profiles are heterogeneous with p v varying in between 700 and 2000 m/s. At depths of 5 and 8 m, profile 1 and large parts of profiles 2 and 3 show low p v values with a discrepancy of 1000 m/s compared to the other profiles. At a depth of 24 m, p v is again homogeneous with values above 2700 m/s for all profiles.

3.2

Underground structure along the SRT profiles

We explain the progressive increase of p v downward by a decreasing of weathered materials with depth, as observed in other sites lying on crystalline or rhyolitic bedrocks [START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF]Holbrook et al., 2014;[START_REF] Olyphant | Topographic correlations with soil and regolith thickness from shallow-seismic refraction constraints across upland hillslopes in the Valles Caldera, New Mexico[END_REF]. p v variations observed from a profile to another at a 5 m depth suggest that the thickness of the weathered medium varies in different areas of the catchment. Results obtained along profile 1 show that the region upstream the main spring presents a thicker weathered zone compared to the rest of the catchment. The same conclusion was previously deduced from MRS measurements showing a region with higher water content [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF]. In Lesparre et al. (2020a), MRS data estimated by the hydrological model NIHM (described below) were fitted to field measurements in order to calibrate the thickness and the porosity of the model. This calibration showed that a thicker weathered zone was required in that same area upstream of the main spring. The SRT data confirm the occurrence of that deeper weathered zone that is not limited around the MRS acquisition station but extends all along SRT profile 1. Our results also reveal that a thicker weathered region is susceptible to occur at the bottom of steep slopes as shown in profiles 2 and 3 (Fig. S1 andS2). Alternatively, weathered materials located in the valley bottom may be relatively thinner than other regions. Discrepancies are noticed from one slope to another, notably along the third profile, but no particular trend can be extracted to distinguish the north-and south-facing slopes.

In the Strengbach catchment, the underground porous material is described by two layers: the soil and the saprolite [START_REF] Gh'mari | Etude minéralogique, pétrophysique et géochimique de la dynamique d'altération d'un granite soumis aux dépôts atmosphériques acides (bassin versant du Strengbach, Vosges, France) : Mécanismes, bilans et modélisations[END_REF]Fichter et al., 1998a;Lesparre et al., 2020a). From the literature, only a few studies have explored the choice of a velocity threshold to delimit the saprolite upper and lower interfaces in such hard-rock contexts. [START_REF] Begonha | Weathering of the Oporto granite: geotechnical and physical properties[END_REF] measured ultrasonic velocities on saprolite and weathered granite samples from Oporto (Portugal). They showed that porosity is the most influential property on the seismic velocity when studying the influence of weathering. Their analysis of 167 samples concluded that the velocity threshold between saprolite and moderately weathered granite is around 2000 m/s. Several field SRT measurements above crystalline bedrocks have confirmed this threshold value by comparing the profiles with pits, borehole logs or images acquired with other geophysical methods [START_REF] Olona | Weathering variations in a granitic massif and related geotechnical properties through seismic and electrical resistivity methods[END_REF][START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF]Holbrook et al. 2014).

Other studies allocated the saprolite bottom interface at the depth where p v exceeds either 1100 m/s, 1200 m/s or 1400 m/s [START_REF] Flinchum | Critical Zone structure under a granite ridge inferred from drilling and three-dimensional seismic refraction data[END_REF][START_REF] Holbrook | Links between physical and chemical weathering inferred from a 65-m-deep borehole through Earth's critical zone[END_REF]. The range of p v in soil is less discussed because SRT is not always efficient in providing information with a fine-enough resolution to study such a thin layer. The resolution depends on the inter-distance between geophones, and for studies exploring the protolith upper interface, long inter-distances between geophones are preferred. Moreover, ultrasonic measurements on soil samples raise issues concerning preserving the in-situ conditions of the medium analyzed. In a similar crystalline context, [START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF] performed SRT using a 1-m spacing between geophones to delimit soil < 0.5 m thick. They estimated that 700 m/s p v  corresponded to the interface between these disaggregated materials and saprolite.

On the Strengbach catchment, different boreholes and pits were excavated to study the soil properties, the structure of the shallow underground CZ, and the erosion processes [START_REF] Ackerer | Regolith evolution on the millennial timescale from combined U-Th-Ra isotopes and in situ cosmogenic 10Be analysis in a weathering profile (Strengbach catchment, France)[END_REF][START_REF] Belfort | Vadose Zone Modeling in a Small Forested Catchment: Impact of Water Pressure Head Sampling Frequency on 1D-Model Calibration[END_REF]. Unfortunately, the pits are distant by more than 100 m from the SRT profiles. We had to consider the steep slopes and the density of the vegetation when designing the layout of SRT surveys. Thus, we initiate our analysis by only examining p v variations along the profiles that are distant by less than 50 m from a borehole to provide an order of magnitude of the p v thresholds at the soil and saprolite interfaces. In that goal, we consider p v values corresponding to the interfaces depth of the soil and saprolite identified when drilling the boreholes (Table 1). The soil thickness is not precisely estimated from the boreholes drilling as it is generally thin at the drilling locations (i.e., around 0.5 m and never above 1 m thick). The p v threshold of the soil bottom interface varies in [410; 720] m/s along profiles 9, 13 and 3, which are distant by less than 35 m from the F1, Pz3 and Pz10b boreholes, respectively (Table 1). The saprolite bottom interface is estimated as the depth where the drilling tool had to be changed as it was penetrating a much less weathered rock. This interface is estimated at a depth of 4.5 m in the Pz3 borehole located close to profile 13 (20 m). For that borehole, the p v saprolite threshold varies in [1480; 2245] m/s (Table 1). The correspondence between the F1 borehole and its closest SRT profile gives a much lower velocity range in [900; 1030] m/s. This lower range can be explained by the comparison between a local measurement of the saprolite bottom location from the drilling, while on the SRT profiles the resolution is of a few meters so local heterogeneities are smoothed. All the more, there is more distance (35 m) between the F1 borehole and its nearest profile compared to the other boreholes and their respective neighboring profiles. The F8 borehole that is close to the 3rd profile is excluded from our analysis since the borehole is located in the valley bottom where the SRT profile shows a strong heterogeneity and, therefore, a much wider velocity range (Fig. S1 and Table 1).

To assess the variability of the soil (saprolite) compartment thickness along each profile, we apply a p v threshold of 700 m/s (2000 m/s) (Fig. S3). We note that the average thickness of 3 m soil in zone 3 is twice as high as in the other zones (Fig. S4).

The average thickness estimated for the saprolite is around 3.5 m in zones 1 and 4, while it reaches 8 m in zone 3 and 12 m in zone 2 (Fig. S4).

4

Construction of 3D p v models from SRT

SRT data filtering

Geostatistical tools are applied to interpolate p v in order to construct 3D p v blocks that could help in defining the geometry of the hydrological model covering the whole catchment. As mentioned above, velocity trends are observed with depth due to weathering processes related to changes in porous material properties along the profiles (Fig. 2). Besides, at a given depth we observe strong p v variations at the catchment scale leading to strong variations of the soil and saprolite thicknesses obtained with a fixed p v threshold (Fig. S3). Since p v maps show non stationary significant variations, SRT data have to be filtered to remove these trends and perform the geostatistical analyses. In that goal, the water catchment is partitioned in zones. The zonation and the filtering of the p v values are performed in four steps:

1.

We analyze the SRT profiles that show strong lateral variations to distinguish eventual locations of zones' boundaries crossing SRT profiles.

2.

Since the SRT profiles do not cover the whole catchment, we add constraints on zones delimitation by considering the soil surface slope (Fig. 3a) and altitude (Fig. 3b). We chose these two variables because we assume that the evolution of the porous material is linked to erosion and weathering processes, which both depend on slope and altitude [START_REF] Riebe | Controls on deep critical zone architecture: a historical review and four testable hypotheses: Four Testable Hypotheses about the Deep Critical Zone[END_REF].

The validity of a such hypothesis is confirmed by [START_REF] Uhlemann | Surface parameters and bedrock properties covary across a mountainous watershed: Insights from machine learning and geophysics[END_REF]. Slope and altitude thresholds are defined from the analysis of a digital elevation model (DEM) characterized by a 0.5 m lateral resolution. The slopes are computed after applying a 40 40 m  rectangular filter to remove the effects of the small-scale asperities of the topography. The thresholds are determined so the zonation is consistent with the lateral variations observed on the seismic profiles (Fig. 3c). We favor a limited number of four zones for having enough data in each zone to compute reliable statistics.

3.

For each zone i , an average velocity at a given depth () i p vz  is computed (Fig. 4). Close to the surface (depth < 2 m), p v distributions are similar from one zone to another (Fig. 4 and S5). Deeper,

4.

Each SRT profile is split according to the zonation (Fig. 3c) and for each sub-profile corresponding to zone i , the residual is computed using

    log10 log10 ( ) i pp w v v z     .
The logarithm of the velocity is used because its distribution is closer to a Gaussian distribution than the velocity itself.

  log10 ( )

i p vz  represents the average log- velocity at a depth z in zone i .

The result of the procedure is presented in Fig. 5a and 5b for profile 2. The trend with depth and the contrast in velocity at the interface between two zones at a distance of 170 m can be seen in figure 5a. After filtering, the residuals do not show any vertical trend but some minor differences still remain at the interface between zones 3 and 4 (green and blue lines above the profiles Fig. 5b).

Geostatistical modelling of the seismic P velocities

In preliminary tests, horizontal and vertical variograms were estimated without considering the zonation of the catchment. In the XY plane, the variogram shows a horizontal coherency (blue line, Fig. S6), but no vertical correlation arises (red line, Fig. S6). Therefore, variograms for horizontal slices of 0.5 m are computed from the surface down to a depth of 25 m.

We chose the truncated power value (TPV) model to fit each experimental variogram because the support volume of SRT measurements increases with distance between two geophones (Di [START_REF] Federico | Scaling of random fields by means of truncated power variograms and associated spectra[END_REF][START_REF] Neuman | On the geostatistical characterization of hierarchical media: characterization of hierarchical media[END_REF]. The TPV model filters out random fields with an integral scale larger than u  and lower than l  (Di Federico and Neuman 1997; [START_REF] Neuman | On the geostatistical characterization of hierarchical media: characterization of hierarchical media[END_REF]. u

 is assimilated to the dimension of the sampling scale -here the catchment size -while l  refers to the data support-in our case the SRT resolution [START_REF] Heße | Generating random fields with a truncated power-law variogram: A comparison of several numerical methods[END_REF][START_REF] Neuman | On the geostatistical characterization of hierarchical media: characterization of hierarchical media[END_REF]. The TPV variogram ( , , )

l u s n n  is defined as: 0 ( , , ) ( , ) ( , ) l u l u s n n c s n s n       , (1)
with s the lag distance, ( , )

2 2 2 ( , ) ( ) [1 exp( ( ) ( ( ) ) (1 , ( ) )] 4 4 4 H m m m m m s n n sn sn H sn            , ( 2
)
where  represents the gamma function, the variance The random fields constitute 3D blocks of 25 m depth and are created with the Geostatistical Software Library (GSLIB; [START_REF] Deutsch | Geostatistical Software Library and Users Guide[END_REF] updated with additional libraries to compute the TPV Gaussian law [START_REF] Neuman | On the geostatistical characterization of hierarchical media: characterization of hierarchical media[END_REF]. GSLIB is a collection of geostatistical programs developed to build variograms, apply kriging and generate stochastic simulations [START_REF] Deutsch | Geostatistical Software Library and Users Guide[END_REF]. The quality of the simulations was checked by looking at the distribution of the simulated residuals (Gaussian distribution with zero mean and prescribed variance) and by computing the variograms of the generated fields (see Fig. 6). The simulations were also verified by removing one by one each SRT profile to compare the distribution of the generated velocity with the removed one. Vertical cross-sections of p v parallel to profile 2 extracted from the generated 3D blocks are illustrated in Fig. S7 together with a map showing their respective locations.

Underground structure of the whole catchment

With p v thresholds of 700 m/s (2000 m/s), we obtain the distribution of the soil (saprolite) thicknesses on the whole catchment from the 3D p v blocks (Fig. 8). The average and standard deviation of the soil and saprolite thicknesses, computed from the 250 geostatistical models, reproduce the zonation division (Fig. 8). As expected, zones 1 and 4 share similar characteristics with soil and saprolite thicknesses of 1.4 0.5 m  and 4 1 m 

, respectively (Fig. 8). In zone 2, the soil thickness increases In the following, we use the SRT data to define the thickness of the aquifer layers used in a hydrological model. We examine then how the estimated thickness uncertainty influences some of the models' outputs: piezometric and MRS data distributed over the catchment. We chose these two variables because one is representative of the saturated zone (piezometric level) while MRS also includes information of the water content in the unsaturated zone. Both simulated data types are estimated at the same location to allow a comparison of their sensitivities. They are located at the same place where field MRS data were acquired. The location zone of those stations, their respective distance with their closest SRT profile and their Topographic Wetness Index (TWI, defined in the appendix) are summarized in Table 2.

The uncertainty on the layers' thicknesses is related to the uncertainty of the SRT data, their conversion in velocities p v , the interpolation of p v over the whole catchment and to the unknown p v threshold values used to define the interfaces between layers. Uncertainties related to the SRT data inversion and to the p v interpolation have been handled in the geostatistical 13 framework described above. The selected p v threshold values correspond to likely values encountered in the literature (Begonha and Bragga, 2002;[START_REF] Olona | Weathering variations in a granitic massif and related geotechnical properties through seismic and electrical resistivity methods[END_REF][START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF]Holbrook et al., 2014) and are in the value ranges estimated when comparing the SRT profiles with the field observations (Table 1). We investigate the impact of the soil bottom location by testing p v threshold values of 500, 700 and 900 m/s, keeping a fixed p v threshold at 2000 m/s to define the saprolite interface (Fig. 9a). Alternatively, we look for the influence of the saprolite bottom interface depth with p v threshold values of 1500, 2000 and 2500 m/s, the soil bottom location being defined with a 700 m/s p v threshold (Fig. 9b). The choice of those values is justified by the bibliographic analysis described in section 3.2.

From the obtained geometries, we estimate the average thickness under each MRS or piezometric stations for each applied p v threshold (Fig. 9). The generated fields correctly reproduce thicker soil under stations located in zone 3 (stations 3 and 7, Fig. 9a) and thicker saprolite in zones 2 and 3 (stations 5, 8, 22; 3 and 7, Fig. 9b) with respect to zone 1 and 4 hosting the other stations. In zones 1, 2, and 3, the soil thickness difference is higher than 1 m when comparing interfaces corresponding to distinct p v threshold values (Fig. 9a). In zone 4, that difference is less than 1 m. The soil thickness standard deviation is globally in the same order of magnitude as the average thickness difference between distinct p v thresholds. The thickness difference between p v thresholds in the saprolite is larger, with an estimated thickness difference higher than 3 m in zones 2 and 3 (stations 3, 5, 7, 8, 22; Fig. 9b). This is slightly above the standard deviation values of the thickness lower than 2 m in such zones.

5 Hydrological model and outputs

The Normally Integrated Hydrological Model -NIHM

The Normally Integrated Hydrological Model (NIHM) is a physically-based model that computes water flows by coupling processes occurring at the surface (1D stream flow and 2D surface flow) and in the subsurface compartments of a water catchment. Meteorological forcing data such as precipitations, evapotranspiration and temperatures are required NIHM inputs.

We describe below the main characteristics of NIHM. A detailed description of the model and its numerical aspects are provided in [START_REF] Pan | A coupled stream flow and depth-integrated subsurface flow model for catchment hydrology[END_REF] and [START_REF] Jeannot | A low-dimensional integrated subsurface hydrological model coupled with 2-D overland flow_ Application to a restored fluvial hydrosystem (Upper Rhine River -France)[END_REF].

The surface flow (1D and 2D) is computed through a simplified formulation of the St-Venant equations, the diffusive wave model, neglecting the inertial effects [START_REF] Panday | A fully coupled physically-based spatiallydistributed model for evaluating surface/subsurface flow[END_REF]. Henderson (1966) consider inertia terms to be negligible in most cases and Ahn et al. (1993) argues that such a simplification induces errors between 5% and 10% that can be treated as negligible in comparison with uncertainties on the meteorological forcing or on the hydrological data. For our application, the option that manages the diffuse 2D surface run-off and exfiltration is switched off as such processes have never been evidenced at the Strengbach catchment. The soil covering the catchment is generally sandy, so it favors rapid infiltration even over steep slopes [START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF].

The diffusive wave formulation writes:

  1/2 5/3 2/3 1 () () r r L r s GM r r A h A h q h h t x x h nP h x                           . ( 3
)
The flow cross-sectional area A [L2] and the wetted perimeter P [L] both depend on the stream geometry. The Gauckler-

Manning coefficient GM n [T/L1/3] is
    , if , ( ) if r s rr h t h z ht z h z       x x x (4)
where h [L] is the groundwater head and r z [L] the riverbed elevation. Initial conditions are defined by initial values of the free surface elevation. Boundary conditions are of Dirichlet or Neuman type. At the outlet, it is assumed that the head gradient is equal to the river bed slope (flow parallel to the river bed also called zero depth gradient).

In the subsurface compartment, we assume that the water flux perpendicular to the substratum is negligible compared to the water flux parallel to the substratum. In other words, we assume that the head is constant along the perpendicular to the substratum. Following this assumption, the 3D Richards' equation is integrated (averaged) over that direction to obtain a 2D flow model. This workaround allows a significant reduction of the meshing effort, the required memory space and the computational cost while preserving the main physics of the flows [START_REF] Weill | A low-dimensional subsurface model for saturated and unsaturated flow processes: ability to address heterogeneity[END_REF][START_REF] Jeannot | A low-dimensional integrated subsurface hydrological model coupled with 2-D overland flow_ Application to a restored fluvial hydrosystem (Upper Rhine River -France)[END_REF]. Comparisons with other hydrological models on benchmarks have shown that this assumption is valid [START_REF] Pan | A coupled stream flow and depth-integrated subsurface flow model for catchment hydrology[END_REF][START_REF] Jeannot | A low-dimensional integrated subsurface hydrological model coupled with 2-D overland flow_ Application to a restored fluvial hydrosystem (Upper Rhine River -France)[END_REF][START_REF] Weill | A low-dimensional subsurface model for saturated and unsaturated flow processes: ability to address heterogeneity[END_REF].

The mathematical model of the subsurface compartment writes:

                    0 ( , ) , , , , ( , 0)
( ) ( , ) ( , ) 0, , ( , ) 0, rs D D s N N s h S t ht h t f t h t h t t hh h t h t t h t q t t                                       x T x x x x x x x x x x x T x u x x (5)
and

            () s w s w bb s w z z wb z z z s s r z z z h h dz S h S z z h h dz dz k h dz                  T K K K (6)
where  [-] is the water content, S [-] the storativity and T the transmissivity tensor [L2T-1], the latter depending on the groundwater head.

  1 0 () () 1 0 r e sr S                      (7)
and the Mualem model [START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF] for the relative hydraulic conductivity r k :   s their influence on the stored groundwater is lesser. The three dimensional distribution of the water content can be computed by NIHM through post-processing, using the constant head assumption (since the head is assumed to be constant perpendicular to the substratum) and ( 7). Water contents can then be used to estimate MRS signals at given stations as described in the next sub-section.

  2 1/ 1 1 0 () 1.0 0 ee re s SS K kS K                  (8)
The equations are solved with a fully implicit non-conforming finite element method that allows a high flexibility of the discretization and ensures continuity of the normal component of the velocity from one element to the adjacent one. Although the subsurface flow model is 2D, it requires an explicit description of the parameters in three dimensions. The computation of the integrals in ( 5) is based on the elevation and slope of the aquifer's substratum. In this paper, this geometry is estimated through seismic refraction data. The medium is divided in two vertical layers: the soil and the saprolite. The thickness of those compartments vary from a mesh to another but we consider that s  and s K are homogeneous in each layer.

The model has already been applied to the Strengbach catchment and showed its capacity to reproduce the behavior of the catchment flows [START_REF] Pan | A coupled stream flow and depth-integrated subsurface flow model for catchment hydrology[END_REF]. NIHM was also used to constrain the distribution of the flow lines in the Strengbach catchment [START_REF] Ackerer | Crossing hydrological and geochemical modeling to understand the spatiotemporal variability of water chemistry in a headwater catchment (Strengbach, France)[END_REF] and to explore the variability of the water transit times through the watershed [START_REF] Weill | Variability of Water Transit Time Distributions at the Strengbach Catchment (Vosges Mountains, France) Inferred Through Integrated Hydrological Modeling and Particle Tracking Algorithms[END_REF]. The comparison between observed MRS data and NIHM deduced MRS estimates was performed on the Strengbach catchment for conditioning NIHM's thickness and s  (Lesparre et al. 2020a).

The equations defining the groundwater flows show that key hydraulic variables such as the transmissivity T and the water content  correspond to the integration over the porous media thickness of the hydraulic parameters () Kh and () h  , respectively as stated in (6). Thus, to solve the inverse problem seeking the hydrological model parameters, misestimating the thickness of the hydrological model underground compartments would inherently lead to a wrong assessment of the hydraulic parameters. The porous media thickness might then be considered as a sought parameter or at least as a prior information associated with an uncertainty. All the more, measurable data sensitive to T and  should be completed with data directly related to the porous media thickness to tackle the porous media thickness correlation with () Kh and () h  in (6).

MRS data estimate

MRS is a non-invasive geophysical method that is classically used to estimate the underground water content in the saturated and unsaturated zones of the subsurface [START_REF] Legchenko | Magnetic resonance sounding applied to aquifer characterization[END_REF][START_REF] Costabel | Noninvasive Estimation of Water Retention Parameters by Observing the Capillary Fringe with Magnetic Resonance Sounding[END_REF][START_REF] Mazzilli | Contribution of magnetic resonance soundings for characterizing water storage in the unsaturated zone of karst aquifers[END_REF]. Thirtytwo MRS measurements were performed on 23 different stations covering the Strengbach catchment during two campaigns in April and May 2013. Data were acquired with a Numis plus device system from IRIS instruments using eight-shaped square loops. This data set was fully described in Lesparre et al. (2020b). A first analysis of the MRS measurements described the subsurface water content distribution over the catchment [START_REF] Boucher | MRS characterization of a mountain hard rock aquifer: the Strengbach catchment[END_REF][START_REF] Pierret | The Strengbach catchment: A multidisciplinary environmental sentry for 30 years[END_REF]. A subset of the data acquired at 16 stations was then used as a posterior information to select subsurface parameters of NIHM applied on the Strengbach catchment (Lesparre et al., 2020a). Here, we estimate MRS synthetic data from NIHM simulations. The MRS signal envelope ( , )

V q t decays with time t during the sounding for a pulse moment q . It can be written as follows [START_REF] Legchenko | A review of the basic principles for proton magnetic resonance sounding measurements[END_REF]:

  * 2 ( , ) ( , ) ( ) exp / ( ) dz z V q t q z z t T z       (9)
where ( , ) qz  represents the kernel function of the MRS vertical sensitivity and depends on the geometry of the acquisition system and the amplitude of the injected pulse q . ( , ) qz  is influenced by environmental conditions such as the geomagnetic field amplitude, the Larmor frequency and the electrical resistivity of the subsurface [START_REF] Legchenko | A review of the basic principles for proton magnetic resonance sounding measurements[END_REF]. The values of the parameters used for the computation of ( , ) qz  are given in Lesparre et al. (2020b). The shape of ( , ) qz T the apparent value of the relaxation time estimated for each pulse (see Lesparre et al., 2020a). Then, we use the () z  values provided by NIHM to compute values of ( , )


V q t with ( 7) and investigate how they evolve with the tested geometries and parameters' sets.

6 Impacts of layer thicknesses on hydrology variables

Test case setup

The influence of the soil and saprolite thicknesses on hydrological variables is analysed using two outputs: piezometric heads linked to the saturated thickness and water content (through MRS) related to the water stored in the saturated and unsaturated media. This influence is quantified by a simplified sensitivity analysis that consists in running the hydrological model NIHM for each 250 simulated velocity fields with the following input parameters: distinct velocity thresholds to define the layers' thicknesses and different sets of hydraulic parameters. We focus our investigations on testing the impact of the hydraulic conductivity s K , the saturated water content s  and air pressure entry  . Preliminary tests showed that the considered outputs (MRS data and piezometric heads) are mainly sensitive to those hydraulic parameters together with the thickness of the underground layers.

In a first step, we fix the set of hydraulic parameters and investigate combinations of soil and saprolite p v thresholds. We 3. The values given to each parameter are defined considering a previous study of the Strengbach vadose zone [START_REF] Belfort | Vadose Zone Modeling in a Small Forested Catchment: Impact of Water Pressure Head Sampling Frequency on 1D-Model Calibration[END_REF].

We note that in similar granitic catchment contexts, porosity values (that we relate to s  ) as high as 50% and 60% have been estimated in the shallow region (Holbrook et al., 2014[START_REF] Holbrook | Links between physical and chemical weathering inferred from a 65-m-deep borehole through Earth's critical zone[END_REF].

Simulations are run with the meteorological forcing measured on the Strengbach catchment from June 1, 2012, to May 31, 2013, as this period covers the MRS measurement campaign. We analyze data estimated at a same date, 19th of April 2013, so we can compare data related to a same meteorological forcing history. This date corresponds to a relatively low water level and only a few artesian locations might be observed. Artesian events might indeed happen depending on the applied parameters, the p v thresholds and the station location. Because NIHM is not designed to simulate these situations properly, we prefer to focus the data sensitivity analysis to an average flow period to limit the occurence of such events and so variations of the water table can still occur.

The head levels are converted to water table depths (WTD). For MRS data, we focus the analysis on the signal simulated for the pulse that shows the largest variability when compared to the other pulses applied on the field. A high water content in the underground induces a high MRS signal amplitude and a water level close to the surface corresponds to a low WTD value.

Results are first presented on 3D plots that represent the projection of the simulations on three planes: thicknesses of both layers (horizontal plane) and MRS signal or WTD values in function of the 2 layers' thickness (vertical planes; Fig. 10 and11). When exploring the influence of the parameters' set, data on the horizontal plane are in grey since the soil and saprolite thicknesses vary with the same distribution for the three studied sets (Fig. 11). We discuss data estimated at stations 5 and 6 which are representative of the main results. Results of all stations are given in supplementary materials (Fig. S8 to S11).

Stations 5 and 6 differ in soil thickness (less than 3 m for station 6, less than 5 m for station 5) and in saprolite thickness (between 1 m and 12 m for station 6, and 5.5 m to 16 m for station 5). Therefore, the total thickness of the aquifer below those stations is not similar (Fig. 9). Station 5 represents zones where the topography favors water storage (high TWI) whereas the topography is propitious to water drainage around station 6 (low TWI).

We then estimate the coefficient of determination of a linear regression 2 R between the estimated data and the soil or saprolite thicknesses for all stations to describe how their specific location influences the data sensitivity to the thicknesses variations (Fig. 12 and13).

2 R highlights a linear relationship between the estimated data and the layer thickness when it is close to 1.

However, a coefficient significantly different from 1 does not mean that the data are not dependent on the layer thickness.

Stations 9, 13, and 14, are located less than 10 m from a seismic profile (Table 2), therefore measured p v values strongly constrain the soil and saprolite thicknesses that are accurately estimated for given velocity thresholds (Fig. 9, S8 to S11). Those narrow variation ranges hinders analyzing the correlation between the soil and saprolite thicknesses and the estimated data so we do not include such stations in the 2 R analysis.

Contrarily to other geophysical methods, MRS is directly sensitive to the underground water content as no petrophysical relationship is required to estimate the MRS signal from water contents estimated by a hydrological model. However, the signal measured on the field is impacted by the instrument dead time, the pulse length and the presence of bounded water cannot be detected. In the analysis applied to synthetic estimates, we did not consider such aspects that influence MRS measurements in addition to the hydraulic parameters' values. They should be taken under consideration in the analysis of real MRS data.

Groundwater variations with respect to the porous medium thickness

For a given set of parameters (e.g., set B in Table 3), we investigate the influence of the p v thresholds on the MRS and WTD values. Note that the p v threshold of the soil layer influences the saprolite thickness: the lower the soil threshold, the thicker the saprolite layer for a same p v threshold of the saprolite layer. Results clearly show the important effect of the station location on the MRS amplitude which varies in [10-100] nV at station 6 and [100-300] nV at station 5 (Fig. 10a andb). WTD values are also strongly impacted as they vary between 1 and 15 m at station 6 and between 1 and 3 m at station 5 (Fig. 10c andd). The thicker underground medium under station 5 and its position on a region favoring storage (high TWI) explain its higher MRS and lower WTD values. The sensitivity of the data to p v is clearly different for these two stations. At station 6, the MRS signal is proportional to the soil thickness for small saprolite thickness (less than 2 m, brown dots Fig. 10a). For higher saprolite thicknesses, the MRS signal is lower and linearly dependent on the saprolite thickness (purple dots Fig. 10a).

The WTD is linearly dependent on the thickness of the saprolite layer, since the WTD is mostly below the soil layer (Fig. 10c).

At station 5, MRS estimates obtained with all p v thresholds show a linear trend with the soil thickness, but none of them show such a trend with the saprolite thickness (Fig. 10b). WTD values do not show any linear dependence on the soil or saprolite thicknesses (Fig. 10d). However, a thicker soil layer is related to a deeper WTD and high MRS values (green dots, Fig. 10d). A thicker soil provides more space to store water leading to a stronger MRS amplitude, but it also increases the transmissivity that might favor drainage and thus reduce the water level. Fig. 10 also highlights non-uniqueness of MRS and WTD with respect to the geometry for a given hydrological parameter set. In particular, at station 5, a given value of WTD can be obtained by different combinations of the layers' thicknesses. It is less true for MRS at the same station where the number of possible combinations is lower due to the correlation with the soil thickness. This clearly shows the interest of using different kinds of measured variables to better constrain the model.

A global overview of the correlations that may exist between layers' thicknesses and MRS or WTD is provided in Fig. 12. For almost all stations, when the correlation with the soil (resp. saprolite) thickness for MSR or WTD is significant, the estimated data are not linearly correlated for saprolite (resp. soil). On average, MRS with 2 R values above 0.5 (Fig. 12a) are more linearly dependent on the soil thickness than WTD for which 2 R values remain mostly below 0.5 (Fig. 12b). WTD is more controlled by the saprolite thickness as 2 R values above 0.5 are observed (Fig. 12d). This can be explained by the fact that WTD depicts the water level of the saturated medium that might remain in the saprolite under dry conditions, while MRS depends on the water content variations in both the saturated and unsaturated media. In general, MRS and WTD better correlate with the soil thickness when the saprolite is thinner (brown lines Fig. 12). On the contrary, both data types better correlate with the saprolite thickness when it is thicker (purple lines Fig. 12). A thicker saprolite hinders the presence of water in the soil as the water level might be lower, and also since it increases the transmissivity and thus favors drainage. Therefore, the influence of the soil thickness on the estimates is annihilated. In contrast, a thin saprolite is more likely saturated by a higher water level and a reduced transmissivity so its thickness influence on the estimates diminishes. For both data types, stations with a low TWI are generally better correlated with the saprolite thickness than stations with a high TWI. A low TWI indicates a region favorable to drainage and thus to a low water level for a given aquifer bottom, therefore the groundwater level is more likely present in the saprolite.

Stations 3 and 7 show lower 2 R values compared with stations characterized with a similar TWI, in particular for MRS (WTD)

values compared with the saprolite (soil) thickness (Fig. 12b andc). Those stations located in zones 3 present thicker soil and saprolite (Table 2 and Fig. 9). The high WTD under those stations does not help exploring the influence of the soil thickness (Fig. S9). WTD and MRS at stations 5, 8, and 22 seem to be independent from the saprolite thickness. Those stations associated to a relatively high TWI are located in zone 2 (Table 2) that is relatively flat and thus propitious for water storage. The WTD underneath those stations is close to the surface and varies in a range of 1 m or less indicating that the WTD is not strongly influenced by the variability of the layers' thickness (Fig. S9). MRS is still strongly correlated to the soil thickness under stations 5 and 8 as MRS depends on the water content in the unsaturated medium and the water table is between 1 and 2 m below the surface at those stations (Fig. S8). However, at station 22 with the highest TWI value, the water table is very close to the surface, when not in artesian conditions (Fig. S9), MRS or WTD cannot be influenced by the underground medium thickness for our parameter sets.

Groundwater variations with respect to the hydraulic parameters

We investigate now the effects of the hydraulic parameters for a given set of p v thresholds of 700 m/s for the soil and 2000 m/s for the saprolite (Fig. 11). Thicknesses variations of the soil and saprolite are thus reduced since they are only related to the generation of the 250 geostatistical models. Despite that, we note that the range of variations of both signals are similar as in the previous test. Again, non-uniqueness occurs as different parameter sets give the same MRS or WTD values for given p v thresholds. However, the relationship between both data types and layers' thicknesses is parameter set dependent. This is clearly shown for MRS and soil thickness at station 5 (Fig. 11b) and for WTD and saprolite thickness at station 6 (Fig. 11c).

At station 6, we observe lowest WTD values for the parameter set C that corresponds also to the highest MRS signal reflecting high saturated conditions (blue dots Fig. 11 a andc). The parameter set C has the lowest saprolite s  and so a smaller storage capacity compared to the two other parameters' sets. This small storage capacity leads to a higher water level in the medium below the station. Set C corresponds also to saprolite layers with the lowest s K that further induces a slower drainage and thus might better maintain the groundwater under the station. All the more, the set C shows the highest s  of the soil which provides a larger space to store water in the unsaturated zone and induce a higher MRS signal.

At station 5, WTD values corresponding to parameter A (red dots) are slightly higher than values estimated with the other parameters (Fig. 11d). However, if the parameter set influences the trend between MRS estimates versus the soil thickness, we do not distinguish a clear impact of the paramater set on the MRS signal amplitude (Fig. 11b) as observed at station 6. This means that above a given aquifer thickness, variations of the WTD due to distinct sets of parameters influence less the MRS signal than variations of the soil thickness.

Influence of the parameter sets on the correlations that may exist between layers' thicknesses and MRS or WTD for all stations is illustrated in Fig. 13. In general, MRS and WTD better correlate with the soil (resp. saprolite) thickness for parameter set C (resp. A). Thus, low values of s  and s K in the saprolite associated with a high s  in the soil that characterizes parameter set C lead to a better correlation of the estimated data to the soil. Such parameters favor a water table closer to the surface and a higher water storage, both allowing a stronger influence of the soil thickness on the WTD and MRS. In contrast, high values of s  and s K in the saprolite of the parameter set A induce a better drainage and thus lead to a lower water level. In that case, WTD and MRS are more sensitive to the saprolite thickness. Stations 3, 7, 8, 5 and 22 show a general lower sensitivity to the medium thickness as mentioned in section 6.2. Here again, this peculiar behavior can be explained by a thicker medium, with higher TWI values for stations 8, 5 and 22.

The complementarity between the piezometric heads and the MRS is again emphasized as both signals are differently influenced by the parameters tested. It is yet difficult from the results we obtain to distinguish what is the respective influence of those parameters on the synthetic data.

Insights from the test case

The structure of the CZ underground compartments can be deduced from geophysical images of the subsurface obtained classically along profiles. By analyzing the geostatistical variations of such images, it is then possible to provide an insight of the subsurface geometry beyond the geophysical profiles. Assuming that patterns of hydrological facies can be deduced from variations of the geophysical properties, it is then possible to build the geometry of distributed hydrological models at the catchment scale from a few geophysical profiles. Here we use SRT data to define the interface geometry between layers with distinct hydrological properties [START_REF] Befus | Seismic constraints on critical zone architecture, Boulder Creek watershed, Front Range, Colorado[END_REF][START_REF] Flinchum | Critical Zone structure under a granite ridge inferred from drilling and three-dimensional seismic refraction data[END_REF][START_REF] Pasquet | Catchment-Scale Architecture of the Deep Critical Zone Revealed by Seismic Imaging[END_REF]. The interpolation methodology used to design the geometry of a hydrological model from such data could be applied to other catchments displaying a heterogeneous thickness distribution of its subsurface compartments. Similarly, Electrical Resistivity Tomography (ERT) could be used in catchments characterized by relatively homogeneous and/or layered compartments. Yet the occurrence of a variable clay component at the catchment scale, as observed on the Strengbach, further complicates the interpolation of ERT data. Furthermore, the relation between the geophysical properties and the subsurface structure requires field measurements to sustain the generalization of boundary delineation between the distinct hydrofacies. Our basic sensitivity analysis explores how measurable data as WTD and MRS are impacted by the variability of the surface compartments' thicknesses related to the geostatistical uncertainty and to the choice of the p v thresholds. We show that this sensitivity is influenced by the values of some key hydraulic parameters (namely the water content at saturation and the hydraulic conductivity at saturation).

It appears from this analysis that MRS is differently sensitive to the subsurface properties than WTD. In particular, MRS shows a better sensitivity to the shallowest subsurface compartment under higher drainage conditions. The drainage capacity of the subsurface can be related to 1) the topography as slopes favor drainage; 2) the aquifer thickness since a thicker aquifer leads to a higher transmissivity for a same hydraulic conductivity and 3) the hydraulic parameters as obviously a higher hydraulic conductivity induces a better drainage. Conspicuously, in well drained areas where the WTD does not reach the soil, its measurement does not supply information on the properties of that superficial medium. In contrast, MRS detects the quantity of unbounded water and we show here that MRS presents globally a higher sensitivity to the shallowest compartment thickness, the soil, in comparison with WTD. MRS acquisition protocols adapted to sound the superficial media should then supply complementary information to WTD with the aim of constraining the hydraulic parameter values [START_REF] Costabel | Noninvasive Estimation of Water Retention Parameters by Observing the Capillary Fringe with Magnetic Resonance Sounding[END_REF]. MRS could then be pertinent to calibrate hydrologic models at the catchment scale notably when acquired along draining slopes and in places characterized with a deeper porous medium. In areas where it is crucial to understand the conditions of recharge, MRS data would then be of great interest.

The test case shows the influence of the hydraulic parameters on the MRS and WTD sensitivity to the subsurface compartments' thicknesses. A global sensitivity analysis should help discriminate the impact of the compartment thicknesses and different hydraulic parameter values on those data. We focused our analysis at a given time, while the variations of the water content in the underground evolve depending on the forcing conditions. It would therefore be interesting to extend the sensitivity analysis to a longer time range in order to assess the impact of the hydrological regime on the data sensitivity.

Conclusions

In this paper we propose a methodology to build the pattern of the three dimensional underground heterogeneity from geostatistical analysis of seismic profiles acquired on the Strengbach catchment. The computation of preliminary variograms shows no vertical coherency of the seismic data, allowing a depth-by-depth analysis. The properties of the experimental variograms reflect the data uncertainty variations with depth, the spatial resolution of the SRT, and the dimension of the underground structures. The porous soil and saprolite compartments are assumed to drive most groundwater flow supplying the catchment outlet studied here. The thicknesses of those layers are deduced by defining p v thresholds from field observations and considering previous studies in similar contexts. The study shows that the average soil and saprolite thicknesses are thinner on the catchment crests, upper slopes, and the valley bottom close to the outlet. At the bottom of steep slopes, the largest soil thicknesses occurred together with high saprolite thickness. In a flat area upstream the creek's main spring, the soil is also relatively thick, and the saprolite appears to be the thickest.

Increasing the p v threshold globally shifts the soil and saprolite compartments' bottom limits downward. Thus, an increase in the p v threshold is equivalent to a transmissivity rise in the different layers (Eq. 5) and an increase of the storage capacity.

This tends to lower the groundwater level and induces higher WTD and lower MRS values for a given set of hydraulic parameters. The sensitivity of the WTD and MRS signal to the porous medium thickness is also depending on the set of hydraulic parameters. For instance, low hydraulic conductivity and porosity of the saprolite favor shallower groundwater levels and higher signal sensitivity to the soil thickness. Beyond the valuable information supplied by SRT on the Strengbach catchment underground structure, this paper also shows the double dependence of data influenced by the water quantity (i.e.

WTD and MRS) to both the hydraulic parameters and the thickness of the porous media. Thus, the model geometry knowledge is crucial to reduce the non-uniqueness of the hydrological inverse problem that would fit such data. SRT measurements should be completed with field observations in pits or on outcrops so they could constrain efficiently the hydrological inverse problem.

The tests applied here demonstrate that piezometric heads and MRS signals display different underground structure sensitivity even when collocated. Such a complementarity is very encouraging for setting up future experiments. Data presently recorded with piezometers could be constructively completed with repeated MRS acquisitions sensitive to the medium porosity. The proposed methodology opens the way for applying hydro-geophysical measurements to constrain underground CZ structures (using SRT) and their hydraulic properties (with piezometers and MRS). This demonstrative application could be easily translated into other watersheds where MRS measurements have been or could be acquired for constraining their hydraulic parameters. The design of the SRT profiles distribution should investigate the different underground morphology susceptible to occur on the catchment. This study's field-based synthetic exploration invites a quantitative global sensitivity analysis to deepen the understanding of the respective impact on the different data types of the hydraulic parameters and their eventual combined effects.

Appendix A

The topographic wetness index (TWI) helps distinguishing the capacity of a station to store or drain the groundwater depending on the geometry of the topography. The TWI depends on the upstream contributing area per unit width orthogonal to the flow direction ( a ) and on the local slope ( b ), and is defined as [START_REF] Beven | A physically based, variable contributing area model of basin hydrology[END_REF]:

TWI ln tan( )

a b     . (A1)
A low TWI value indicates a region suitable to drainage while higher TWI values correspond to areas favoring water storage.

We compute TWI values at each MRS station (Fig. 1, Table 2) to classify the obtained results and sustain the data sensitivity interpretation. The sensitivity might indeed be influenced by the spatial configuration of the measurement stations that strengthens a groundwater drainage or storage behavior. 3).

  and estimated travel times, respectively and obs

  to a depth of about 7 m. In the remaining two zones, p v increases faster in zone 2 down to a depth of 7 m, where p v starts to increase faster in zone 3 instead.

c

  corresponds to the nugget and is directly determined by the variance of the 104 seismic results obtained for each profile with 0 s  . TPV models can be characterized either by a Gaussian or an exponential variogram. In our case the Gaussian TPV variogram better fits the experimental variogram and writes:

  2



   is the Hurst coefficient[START_REF] Hurst | Long-term storage capacity of reservoirs[END_REF] and C is a constant. m represents either the index u or l , of the upper or lower wave number, respectively.One theoretical variogram is estimated in each 0.5 m horizontal layer and we analyze the evolution of each variogram characteristics with depth. In that goal, we compute the values of p s that correspond respectively to the abscissa and ordinate of the point where the variograms reach a plateau (yellow stars, Fig.6). to the variance of the variogram. p s corresponds to the projected lag distance where the theoretical variogram reaches () p s  and is related to the correlation length of the TPV variogram. p s is constant until a depth of 7.5, where the variable jumps abruptly before decaying progressively (Fig.7a). m, then it decreases with depth (Fig.7b). The nugget, 0 c , shows a strong decrease between the surface and a depth of 1 m, where it stabilizes until a depth 5 m before it increases with depth (Fig.7c).

c

  variations are influenced by the acquisition geometry of the SRT data. Since the sensors are installed on the surface, the resolution is more accurate in the shallow medium in between 2 and 6 m depth. Smaller targets can be detected near surface so smaller p s values are observed. This better accuracy is confirmed by the lowest 0 c values, and the largest () p s  reflecting the medium heterogeneity. The regularization process used during the SRT inversion involves smoothing the p v distribution. The less-constrained deeper region is depicted by more laterally extended (higher p s values) and blurred targets limited resolution of SRT in the very shallow media explains the low () p s  values in the medium close to the surface. It is impossible to resolve targets with a smaller size than the distance between the geophones. The depth of 3 m at which () p s  is maximum is similar to the geophones inter-distance (i.e., 2 m). This explains as well the higher values of 0 c in the first meter below surface compared to the underlying region.Beyond the acquisition geometry and the characteristics of SRT images related to the inversion process, can be explained by the structure of the underground medium. The p s abrupt jump could be related to the transition where the medium becomes more coherent. In the shallow region, the strongly weathered medium is composed of materials presenting smaller characteristic sizes than in the deeper part. Furthermore, to the presence of roots and pebbles with various dimensions in the shallow region that could induce a strong heterogeneity in the medium.The geostatistical fields are generated following the theoretical TPV model fitted at each depth, and each generated geostatistical field reproduces the variable  corresponding to the normalization of w  . The white noise is added to the residual w to take care of the uncertainty on p v . is estimated from the 104 different velocity tomography computed with distinct inversion configurations. has a Gaussian distribution with zero mean and a variance equal to the variance of the

  The deduced structure in each zone can then be used to delimit the compartment interfaces in the hydrological model NIHM. The zonation methodology we apply induces sharp contrasts of soil and saprolite thicknesses along the zones' boundaries. They reflect the strong lateral p v variations observed notably along the SRT profiles 2 and 3.

  r k is the relative hydraulic conductivity, K [LT-1] and s K [LT-1] represent the hydraulic conductivity tensor and the hydraulic conductivity tensor at saturation respectively. For our application, we consider that those tensors are isotropic, so they are reduced to the scalar values K and s K , respectively. b z [L] is the substratum elevation, w z [L] the groundwater free surface elevation and s z [L] the soil surface elevation. In (5), f [LT-1] is the sink-source term including groundwater and the last term describes the exchange with the river.  is the model domain; boundaries  that correspond to Dirichlet and Neumann conditions, respectively. u is the unit vector normal to the boundary, counted positive outward. conditions defined over the domain and s  is the simulated period. For each element of the catchment model and at each observation time, NIHM provides the water pressure hz   [L] and estimates of  and K based on the van Genuchten model for the water retention (van Genuchten, 1980):

  where e S[-] is the effective water saturation, r  [-] and s  [-] the residual and saturated volumetric water content respectively.  [L-1] (air entry pressure) and  [-] are the Mualem van Genuchten shape parameters, and 1

  assess soil p v threshold values of 500, 700 and 900 m/s for a fixed p v threshold at 2000 m/s in the saprolite. We also examine saprolite p v threshold values of 1500, 2000 and 2500 m/s with a soil p v threshold fixed at 700 m/s. Thus, we test five combinations of p v thresholds, each shifting the soil and saprolite thickness patterns and influencing the global porous volume of the CZ underground compartments as well as their transmissivity. In a second step, we prescribe the p v thresholds to 700 m/s for the soil and 2000 m/s for the saprolite and apply three different sets of hydraulic parameters detailed in Table

Figure 1 :

 1 Figure 1: Map of the Strengbach catchment. The seismic profiles are indicated by red lines, and the flow measurement station at the outlet is represented by the blue triangle. The colored stars show the TWI, defined in Eq. (A1), computed at each MRS station. The white triangles represent the 2D mesh of the hydrological model underground compartment, with the grey zone 2 that was identified as a storage area in Lesparre et al. (2020a). The cyan lines represent the 1D mesh of the hydrological model surface compartment that includes flows in the creek but also on the forestry roads. 870 An inset shows the location of the Strengbach catchment in the Northeast of France.

Figure 2 :

 2 Figure 2: Spatial distribution on the Strengbach catchment of the p v extracted at different depths from the SRT

Figure 3 :

 3 Figure 3: Analysis of the Strengbach catchment topography to delimit zonation in which the p v presents geostationary

Figure

  Figure 4: Average p v as a function of the depth in each zone. The shaded areas represent the average p v more or less

Figure 6 :

 6 Figure 6: Experimental variograms (red lines) estimated from the detrended variable w (red line). The theoretical

Figure 7 :

 7 Figure 7: Characteristics of the theoretical variograms as a function of depth: p s (a), () p

Figure 8 :

 8 Figure 8: Statistical characteristics of the lower boundary of the soil (a, b) and the saprolite (c, d). The averages (a, c) 895 and the standard deviations (b, d) are estimated from the generation of 250 geostatical models following a Gaussian truncated power value geostatical model. The black dots represent the locations of the SRT profiles; the black stars correspond to the MRS station locations.

Figure 9 :

 9 Figure 9: Variation of the soil and saprolite thicknesses below each piezometric and MRS station. The thicknesses 900 plotted represent the average estimated from the 250 generated fields, and the error bars correspond to the thicknesses' standard deviations. Stations are ordered with a crescent TWI.

Figure 10 :

 10 Figure 10: Distribution of the MRS and WTD values as a function of the soil and saprolite thicknesses (labeled th.soil and th.saprolite) below measurement stations 6 and 5. Data are estimated the 19th of April 2013 for different velocity 905 thresholds and the fixed set of parameter B (see Table3).

Table 1 : Depth of the Bottom Interfaces Estimated at the Boreholes and the Corresponding p v Ranges of the Closest Part of the Seismic Profiles at such Depths. The Saprolite Bottom Interface Is Not Intercepted by Pz10b. 860

 1

	Borehole name	Closest profile number	Minimum distance to the closest profile (m)	Bottom interface depth (m) Soil Saprolite	Corresponding p v range (m/s) Soil Saprolite
	F1	9	35	0.5	1.5	480; 720	900; 1030
	Pz3	13	13	1	4.5	410; 630	1480; 2245
	Pz10b	3	9	0.5	-	560; 650	-

Table 3 : Parameters Applied in Each of the Subsurface Compartments for the Different Sets of Simulation Runs

 3

	Parameter		s  (%)		s K (m/s)		 (m-1)
	Medium	Soil	Saprolite	Soil	Saprolite	Soil	Saprolite
	Set A	0.1875	0.08	10-4.5	10-4.5	0.575	1.525
	Set B	0.325	0.06	10-4	10-5	1.05	1.05
	Set C	0.4625	0.04	10-3.5	10-5.5	1.525	0.575

Acknowledgements.

The authors warmly thank Marie-Claire Pierret, Hélène Jund, Marie-Anne Churka, Quentin Chaffaut, Flore Rembert, Sylvain Weill, Benjamin Belfort, Matthias Oursin, Jérôme Vergne, and Sylvain Benarioumlil for their help in the acquisition of the seismic data presented hereThe digital elevation model was obtained from aerial LiDAR images recorded in 2018 by the Helimap system. Computing time was provided by the HPC-UdS.

Financial support.

Meteorological and flow rate data collection was funded by the Observatoire Hydro-Géochimique de l'Environnement (OHGE), which is financially supported by CNRS/INSU France and the University of Strasbourg. OHGE is part of the OZCAR research infrastructure (http://www.ozcar-ri.org). . The CRITEX ANR-11-EQPX-0011 project provided the instruments used for the seismic data acquisition, and the ANR HYDROCRIZSTO-15-CE01-0010-02 project funded the field campaigns.

Data availability

The whole seismic data set is available on the H+ database (http://hplus.ore.fr/en/; [START_REF] Pasquet | Seismic profiles measured on the Strengbach catchment (OHGE observation site), France [Data set[END_REF], which stores the geophysical data collected on the CZ observatories of the OZCAR network. The fortran libraries developed to run the NIHM code are available upon request to the authors. The python library used to analyze the SRT data can be downloaded at 25 www.pygimli.org [START_REF] Rücker | pyGIMLi: An open-source library for modelling and inversion in geophysics[END_REF]. The geostatistical software library, GSLIB, is distributed at www.gslib.com (Deutsch and Journel, 1992).

Competing interests. The authors declare that they have no conflict of interest. 3.


Philippe M Oger

T Radiotolerans

T Siculi

Complete Genome Sequence of 11 type-species from the Thermococcus genus of hyperthermophilic and piezophilic Archaea

genome sequences of the type strains of the species Thermococcus barossii, T. celer, T.

P6.

Thermococcus celer were isolated from surface hydrothermal vent systems in the Mediterranean [START_REF] Zillig | The archaebacterium Thermococcus 49 celer represents, a novel genus within the thermophilic branch of the archaebacteria[END_REF].

Thermococcus sp. Strain P6 was isolated in 1997 by Pr. M. Thomm (University of Kiel, Germany) in 1997 from an hydrothermal site off Palaeochori Bay (Milos, Aegean Sea, Greece). T. thioreducens was isolated from the Rainbow hydrothermal vent system on the Mid-Atlantic ridge, 2600m below the ocean surface (2).All other species were isolated from hydrothermal vents systems in the Pacific ocean (3-9). All species exhibit very similar growth characteristics, with optimal conditions for growth ranging from 83 and 88°C (T), pH 5.8 and pH 7 and salinity between 2 and 3%.

DNA was extracted from cells grown in TRM medium to late exponential phase. Whole-genome shotgun sequencing was carried out by a combination of IonTorrent PGM (800 bp fragment library, 318 Chip, HiQ chemistry) and Illumina HiSeq run (500bp library, 100bp paired reads), which generated at least 583k reads with a mean 389 bp length (IonTorrent) and 13 million paired 100 bp reads after QC check and trimming (10-11) for each species. The IonTorrent data were assembled de novo with Newbler 2.8 (12). The paired-end Illumina data was assembled de novo with velvet [START_REF] Zerbino | Velvet: algorithms for de novo short read assembly using de 83 Bruijn graphs[END_REF]. The pooled IonTorrent and Illumina data was assembled de novo with MIRA (14) and Newbler. The hybrid data assembly yielded betwen 1 and 6 contigs which were connected by PCR.

Chromosome sizes range from 1.52Mb in Thermococcus sp. strain P6 to 2.07Mb in Thermococcus 21 thioreducens. which is similar to known genomes from this family. Only two species, T. peptonophilus 22 and T. profundus, harbor an extrachromosomal element. The T. profundus plasmid belongs to the pTBMP1 family of plasmid described in the piezophilic archaeon T. barophilus (15). To the exception of 24 T. chitonophagus, the GC% ranges between 51.7 and 56% Similar to previously sequenced 25 Thermococcales, all genomes harbor a single copy of the 16S-23S operon, two copies of 5S operons and 26 45 to 48 tRNA genes, large CRISPR loci associated with cas genes (cas1, cas2, cas3, cas4, cas5, cas6, 27 csm2), indicating multiple contacts with mobile elements of viral or plasmidic origin. Sinteny analyzes 28 confirm the plasticity of Thermococcales genomes. The gene content is consistent with an anaerobic and 29 heterotrophic lifestyle. The genome of T. barossi encodes the capacity for hydrogenogenic growth on 30 carbon monoxide, as identified previously in T. onnurienus and T. barophilus. 31

The genome sequence of T. peptonophilus (GC%, core genome phylogeny) confirms previous 32 reports that it may have been misplaced into the Thermococcus genus, and that it should be placed in the 33 Pyrococcus genus instead. Strain P6 is most closely related to Thermococcus celer strain Vu 13 (1), as 34 indicated by in silico DNA-DNA hybridization of their genomes (predicted value = 28.90% using the 35 genome-to-genome distance calculator GGDC2.0 (16). This genomic distance suggests that strain P6 36 defines a novel species inside the Thermococcus genus. 37

Nucleotide sequence accession number. These whole Genome Shotgun project has been deposited in 38

GenBank under the accession numbers CP014750, CP014751, CP014854, CP014855, CP014862, 39 CP014863, CP015101, CP015102, CP015103, CP015104, CP015105, CP015106, and CP015193. 40 41

Fonctionelle de Lyon by Dr. Sandrine Hughes and Benjamin Gillet.

ACKNOWLEDGEMENTS 42

This work was supported in part by the Agence Nationale de la Recherche (ANR-10-BLAN-43

Living deep) and a grant from the Fond Recherche of the ENS Lyon to PO. Ion Torrent