Catchment‐Scale Architecture of the Deep Critical Zone Revealed by Seismic Imaging

Weathering and erosion processes are crucial to Critical Zone (CZ) evolution, landscape formation and availability of natural resources. Although many of these processes take place in the deep CZ (∼10–100 m), direct information about its architecture remain scarce. Near‐surface geophysics offers a cost‐effective and minimally intrusive alternative to drilling that provides information about the physical properties of the CZ. We propose a novel workflow combining seismic measurements, petrophysical modelling and geostatistical analysis to characterize the architecture of the deep CZ at the catchment scale, on the volcanic tropical island of Basse‐Terre (Guadeloupe, France). With this original workflow, we are for the first time able to jointly produce maps of the water table and the weathering front across an entire catchment, by means of a single geophysical method. This integrated view of the CZ highlights complex weathering patterns that call for going beyond “simple” hillslope CZ models.

et al., 2019; Rempe & Dietrich, 2014). This leads to the exposure of new primary minerals (e.g., plagioclase and pyroxene) to reactive meteoric waters (Ackerer et al., 2021), thus accelerating chemical weathering of the rock (Brantley & Lebedeva, 2021;Lebedeva & Brantley, 2020). Around the water table, higher concentrations of dissolved O 2 and CO 2 especially promote dissolution and oxidation reactions which consume primary minerals and release mineral nutrients (Gu, Heaney, et al., 2020;Gu, Rempe, et al., 2020). As weathering intensifies, the bedrock material loses much of its mechanical strength and turns into saprolite, a friable layer yet physically intact enough to retain the texture of the parent bedrock (Graham et al., 2010), containing many new, less soluble secondary phases (e.g., clays, oxide). Close to the surface, the weakening of the rock structure through bioturbation processes eventually leads to the formation of soils (Wilford & Thomas, 2013).
Our understanding of the processes controlling the functioning and evolution of the CZ remains obscured by the difficulty of accessing information at depth. The CZ architecture (i.e., the depth of the water table and the main weathering interfaces) is often characterized with boreholes and piezometric wells, but data remain limited by costs, field access, spatial coverage and the destructive nature of such measurements Hubbard & Linde, 2011;Maillot et al., 2019). Minimally invasive surface-based geophysical methods can be used to fill spatial gaps between wells by producing higher lateral resolution and lower cost data (Hubbard & Linde, 2011;Parsekian et al., 2015). Over the past 10 years, an increasing number of studies have used near surface geophysics to image the architecture of the CZ (Befus et al., 2011;Braun et al., 2009;Comas et al., 2019;Eppinger et al., 2021;Leopold et al., 2013;Novitsky et al., 2018;Olona et al., 2010;Orlando et al., 2016;Parsekian et al., 2021;Pasquet, Bodet, Longuevergne, et al., 2015;St. Clair et al., 2015;Wang et al., 2021Wang et al., , 2022Yaede et al., 2015). These methods measure different geophysical parameters (e.g., electrical resistivity, seismic velocity), which are sensitive to a wide range of physical properties of the subsurface materials (e.g., porosity, density, water or clay content). To disentangle the cumulative effects of these parameters, a growing number of studies have used petrophysical relationships to convert these geophysical measurements into quantitative estimates of subsurface properties (Callahan et al., 2020;Flinchum, Holbrook, Rempe, et al., 2018;Gase et al., 2018;Holbrook et al., 2014;Pasquet et al., 2016).
In most cases, the authors use a single geophysical observation and assume dry subsurface conditions to estimate porosity, thus not taking into account variations of water content as a function of depth (Holbrook et al., 2014). To incorporate information about water saturation, they mainly rely on additional geophysical or piezometric data (Boucher et al., 2009;Buchanan & Triantafilis, 2009;Flinchum et al., 2019;Hayes et al., 2019;Linde et al., 2007). However, conducting these extra measurements can be cumbersome, especially when investigating remote field sites, with limited vehicle access, steep slopes and dense vegetation, as it involves carrying a few hundred kilograms of equipment. In these conditions, one often has to choose between two main surveying strategies: (a) a catchment scale characterization that consists in deploying a single geophysical method over several profiles, so as to cover most of the watershed, or (b) a hillslope scale characterization that employs different methods along a single profile, so as to combine several geophysical parameters and better inform subsurface properties. While most of the works mentioned above employ the latter strategy, the few studies involving catchment scale surveys are only able to resolve the CZ structure assuming a constant water table level (Braun et al., 2009;Befus et al., 2011;Flinchum, Holbrook, Rempe, et al., 2018;Leopold et al., 2013).
Here we propose a novel workflow that aims to simultaneously characterize the CZ structure and the water table depth, at the catchment scale, using a single geophysical method. Several recent studies have shown that structure and water content information can be simultaneously inferred from active seismic data Grelle & Guadagno, 2009;Pasquet, Bodet, Longuevergne, et al., 2015). These studies rely on the joint analysis of pressure (P) and shear (S) wave velocities (V P and V S , respectively). These geophysical parameters are sensitive to subsurface mechanical properties which vary according to changes of porosity, bulk density or water content (Pride, 2005). V S is less sensitive to changes in water content than V P (Biot, 1956a(Biot, , 1956b, since shear waves do not propagate in fluids. Combining both velocities in a petrophysical inversion framework thus allows us to better inform the relative effects of lithology and water content on these measurements (Pasquet et al., 2016). Our workflow couples seismic refraction tomography (SRT) and multichannel analysis of surface waves (MASW) to simultaneously estimate V P and V S velocities from a single seismic data set. These velocities are then converted into spatial distributions of subsurface porosity and saturation using a petrophysical inversion framework. We apply this framework on 5 seismic profiles collected in contrasted parts of a small (8 ha) forested watershed on the volcanic tropical island of Basse-Terre (Guadeloupe, France). We then extend this information at the catchment scale using a kriging interpolation. With this original workflow, we provide a unique characterization of the deep CZ architecture. For the first time, we jointly produce maps of the water table and the weathering front across an entire catchment, by means of a single geophysical method.

Site Description
The Quiock Creek watershed is an 8 ha headwater catchment located on the windward side of Basse-Terre Island, the volcanic part of the Guadeloupe archipelago (France) in the Lesser Antilles ( Figure 1). The catchment is monitored by the ObsErA observatory which is part of the OZCAR critical zone research infrastructure  and is dedicated to the study of weathering and erosion processes in the CZ under tropical climates Dessert et al., 2015Dessert et al., , 2020Guérin et al., 2019). Indeed, due the volcanic nature of the rock formations and the tropical climate of Guadeloupe, weathering rates (i.e., the rate of transformation of rock into saprolite and soil) are amongst the highest on the planet (Gaillardet et al., 2011). The Quiock catchment is representative of volcanic tropical landscapes which are known to be hotspots of nutrient production, biological productivity and soil CO 2 consumption by chemical weathering (Dessert et al., 2001;Louvat & Allègre, 1997). Quiock Creek is a small tributary of the Bras David river and is located in the primary tropical rainforest of Guadeloupe National Park with mean annual temperature of 25°C. The elevation of the catchment ranges from 200 to 350 m and includes an active river knickpoint (i.e., sharp convexity in an otherwise concave-up longitudinal river profile). The knickpoint is located at ∼250 m upstream of the confluence with the Bras David river and separates regionally extensive low relief surface upstream of the knickzone, from more deeply incised streams downstream (Sak et al., 2018). The hydrology of the site is strongly influenced by tropical storms and hurricanes, with mean annual precipitation of 3,500 mm/yr, evapotranspiration between 60% and 70%, and runoff of 1,130 mm/yr (Dessert et al., 2020). Water table levels are continuously monitored upstream of the knickpoint with pressure transducers installed in 7 piezometric wells arranged along a 30-m linear transect perpendicular to the stream ( Figure 1) and vary between 0.1 and 4.1 m in the farthest well. The catchment lies within Pleistocene andesitic pyroclastic deposits (Boudon et al., 1988) that are weathered into a very thick (>10 m) regolith profile with ferralitic soils at the surface (Buss et al., 2010). This weathering profile is highly depleted in mobile elements  and is mainly constituted of secondary minerals, with about 66% clay (halloysite and kaolinite) and 28% iron and aluminum hydroxides (magnetite, goethite, maghemite and gibbsite), minor amounts of primary minerals (mostly quartz and cristobalite) making up the rest of the regolith composition. Bulk density measured in auger samples in the upper 5 m is particularly low, on average ∼1 g.cm −3 (Buss et al., 2010).

Seismic Data Acquisition and Processing
Seismic data were first collected in May 2016 along a 188-m-long profile (P5) crossing the Quiock Creek near the piezometric wells ( Figure 1). Four supplemental 142-m-long profiles (P6 to P9) were collected in May 2019, both up and downstream from the original P5 transect. For each profile, P-wave first arrival times were picked manually, then inverted for subsurface P-wave velocity (V P ) structure using the seismic refraction tomography (SRT) code included in the Python geophysical inversion and modelling library pyGIMLI (Rücker et al., 2017). The program starts with an initial 2D model consisting of a velocity field that increases linearly with depth, and then finds an appropriately smooth update to the model that reduces the difference between predicted and observed traveltimes (more details in the Supporting Information S1). These traveltimes are compared for each source-receiver pair to check the quality of the inversion ( Figure S1 in Supporting Information S1). Velocity uncertainties ( Figure S2 in Supporting Information S1) were estimated by running 144 inversions for each profile, using a different set of starting model and regularization parameters for each inversion run (St. Clair et al., 2015;Pasquet et al., 2016). The 144 inverted models are merged to build an average velocity model describing the V P distribution at depth along each profile (Figures 2b-2f).
The seismic data were also processed to perform multichannel analysis of surface-waves (MASW) using the SWIP software package (Pasquet & Bodet, 2017). SWIP uses spatial windowing and spectral stacking techniques (Neducza, 2007;O'Neill et al., 2003) to extract surface-wave dispersion data from the seismic records and retrieve a 2D model of shear-wave velocity (V S ) from a succession of 1D inversions. We specifically apply the novel multiwindow weighted stacking of surface-wave procedure  to extract higher quality dispersion data and improve both the lateral resolution and depth of investigation of the models (more details in the Supporting Information S1). Dispersion curves are picked manually along each profile, and then inverted using the neighborhood algorithm (NA) with the open software package Geopsy (Wathelet et al., 2004). For each extracted dispersion curve, the NA inversion procedure generates 25,000 models which are used to build a 1D misfit-weighted final V S model. The overall quality of these inversions is quantified by computing the residuals between observed and calculated dispersion curves along each profile ( Figure S3 in Supporting Information S1). For each profile, all the consecutive 1D V S models are finally assembled to create a 2D V S section with the depth of investigation estimated from the uncertainties of each NA inversion ( Figure S4 in Supporting Information S1).

Petrophysical Inversion
The inversion strategy relies on a petrophysical model based on Hertz-Mindlin contact theory (Mindlin, 1949) that describes the weathered regolith as a pack of spherical beads. Clayey materials can still be modelled as a random pack of spherical grains, as long as clays are dispersed in the rock (Mavko et al., 2003). Here we assume that in-place weathering results in dispersed clays throughout the medium, as opposed to clastic deposition sequences which would lead to laminated shales. This is consistent with the textural description of highly weathered saprolite in tropical volcanic islands (White et al., 1998). With this model, we can express the bulk elastic properties of the regolith (i.e., bulk and shear moduli) as functions of the elastic properties of the minerals and fluids constituting the medium, and of their relative proportions (i.e., mineralogy, porosity and saturation). Here, we assume that all the beads in the model have an identical mineral composition that corresponds to the average composition of the regolith (66% clay, 28% hydroxides and 6% quartz) observed in direct samples collected at the site (Buss et al., 2010). A complete description of the model and the calibration of its parameters is given in Supporting Information S1. The Hertz-Mindlin petrophysical model is then used in a grid search inversion scheme to look for the best set of porosity (Φ) and water saturation (W) values that minimizes the differences between observed and modelled V P and V S ( Figure S6 in Supporting Information S1), considering errors previously estimated with SRT and MASW inversions. The grid search procedure was performed with porosity and saturation ranging between 0 and 1 in every cell of the five geophysical profiles presenting valid values of both V P and V S . It allowed us to reconstruct 2D sections of porosity ( Figure S7 in Supporting Information S1) and saturation (Figures 2g-2k), and evaluate their uncertainties ( Figure S8 in Supporting Information S1).

Describing the Critical Zone Architecture
Here we propose to use seismic data to image the water table along with the physical interfaces that shape the CZ. From depth to the surface, these interfaces define the following five distinct layers, as described by Riebe et al. (2017). has enough open fractures to allow water to be drained. Reactive meteoric waters can reach fresh rock at depth, allowing significant weathering between fractures and eventually leading to the formation of large corestones of fresh rock surrounded by altered materials. (d) Saprolite corresponds to friable, altered rock that has lost most of its primary minerals, yet keeping the original fabric of the parent rock. (e) Soil is eventually made of disaggregated materials composed mostly of secondary minerals and organic matter. Saprolite and soil are usually grouped as a common unit called regolith , which is separated from the underlying weathered bedrock by the so-called weathering front.
The depths of these interfaces have increasingly been determined by V P (Flinchum et al., 2022). However, when one of these interfaces lies in the vicinity of the water table, V P is also influenced by variations of subsurface water content that can bias the estimation of the corresponding interface depth .
Here we propose to use the results of the petrophysical inversion to map the depth of the water table and thus unbias the V P -based estimation of the CZ interfaces. We define the depth of the water table when water saturation reaches a critical value (W c ) of 0.9, which corresponds to the top of the capillary fringe (de Marsily, 1986). We then consider an unbiased V P -based estimation of a given interface depth as long as it is located below the water table (i.e., in fully saturated conditions). St. Clair et al. (2015) have pointed out that fresh bedrock is usually characterized by V P > 4,000 m/s, whereas fractured bedrock is expected to have V P > 2,700 m/s in volcanic rocks (Adelinet et al., 2018). Several recent studies also described saprolite with V P < 1,200 m/s (Flinchum, Holbrook, Rempe, et al., 2018;Hayes et al., 2019). Soils express a large range of V P which depends on their compaction and saturation levels, but generally show V P < 800 m/s (Santamarina et al., 2005).

Interpolating the Weathering Front and the Water Table
Since only the weathering front and the water table are clearly identified in all the seismic profiles, we focus the following section solely on reconstructing the 3D shape of these two interfaces across the catchment. We first extracted, from the digital elevation model (DEM) (Figure 3a), the spatial coordinates of these two interfaces at each point along the seismic profiles. We also added boundary conditions in the Bras David river to better constrain the interpolations, assuming: (a) a mean regolith thickness of 2 m that coincides with the most lowland values observed along P9 (Figure 2f), and (b) a mean water level of 0.5 m, considering that the aquifer is directly connected to the river. We then used the GSTOOLS python library (Müller et al., 2022) to interpolate the weathering front and the water table across the catchment (more details in the Supporting Information S1). We specifically applied ordinary kriging along a regular grid of 10 × 10 m cells covering the entire watershed in order to generate 3D surfaces of both the depth (i.e., vertical distance under the surface) and the elevation (i.e., vertical distance above sea level) of these interfaces. As shown by Snyder (2008), the water table position is better constrained by combining interpolations of both its depth and elevation. The interpolation of the water table elevation is more sensitive to the main trend associated with the regional hydrological gradient, whereas the interpolation of its depth helps reconstructing local perturbations associated with land surface irregularities. Following Snyder (2008), we used the average of these two interpolations to produce a map of the water table depth (Figure 3c) that incorporates both local and regional information. The same strategy was applied to interpolate the depth of the weathering front across the catchment (Figure 3b).

Characterization of the Critical Zone Architecture
In the upper 2-12 m, V P are mostly <1,200 m/s which is characteristic of clays constituting saprolite in highly weathered terrains such as those on Basse-Terre (Buss et al., 2010). As this 1,200 m/s threshold is always located below the estimated water table (Figure 2), we consider that the estimation of the weathering front based on this threshold is not biased by variations of subsurface water content. The seismic data clearly underline a deepening of the weathering front in the knickzone (Figures 2c-2e), in comparison to upland and lowland areas (Figures 2b  and 2f, respectively). The interpolated depth of the weathering front ( Figure 3b) reveals a similar pattern, with a clear thickening of the regolith (>15 m) in the vicinity of the knickzone. In upland and lowland areas, the  (in blue) and displays the structure of the CZ with specific V P contours describing the weathering front (in orange), and the transition zone between weathered and fractured bedrock (in brown). Colored dots and associated error bars correspond to the average depth of each interface, extracted from geophysical profiles. Empty brown squares with no error bars correspond to the maximum investigation depth at which fractured bedrock is not reached. See Supporting Information S1 for details about error bar estimation.
weathering front is closer to the surface at depths <5 m. This particular organization is summarized in a synthetic cross-section (Figure 3f) where the depth of each interface identified in the seismic transects is represented along the topographic profile of the Quiock stream at their respective location (i.e., at the intersection or at the closest point in the stream). Soil thickness could not be precisely determined due to the large spacing between geophones (>2 m) used for both seismic campaigns, and is therefore not discussed further.
The velocity threshold V P > 2,700 m/s, associated with the transition zone between weathered bedrock and fractured bedrock, is only reached upstream of the knickzone (Figures 2b-2d), at depths of about 40-50 m. Therefore it could not be interpolated across the whole catchment, and is only interpreted along the synthetic cross-section of the stream (Figure 3f). In the knickzone (Figure 2e) and downstream of the knickzone (Figure 2f), V P is always <2,700 m/s over the whole investigated area, which goes as deep as 50 m in P8. Although we cannot image the base of the weathered bedrock in these areas, these results show that weathered bedrock extends at least down to 45 m in the lower end of the knickzone, and down to 30 m in lowland areas. As no obvious deepening of the weathered-to-fractured bedrock transition zone is detected, we can hypothesize that this interface is located at a depth of about 40-50 m across the whole catchment, roughly following the landscape surface topography. Intact bedrock, usually described with velocities of 4,000 m/s, is never reached in the catchment, and is therefore located at depths >50 m. Such a thick weathered zone confirms deep drilling observations made in similar tropical volcanic environment, where saprolite has been observed down to about 40 m (Buss et al., 2013).
Inverted saturation cross sections provide estimates of water table levels that can be compared with field observations. Upstream of the knickzone (Figures 2g-2i), the estimated water table outcrops at the intersections between the seismic profiles and the stream, in which water was flowing during the field campaign. The water table also outcrops in the small tributary crossed at the southern end of profile P8, whereas it goes deeper (∼10 m) in the knickzone (Figure 2j). Downstream of the knickzone, the water table outcrops again widely (Figure 2k) in an area that was saturated during the May 2019 campaign. All these observations are consistent with the expected shape of the free water table in this small tropical catchment (Guérin et al., 2019). Interpolated water table elevations ( Figure S10 in Supporting Information S1) show a W-E trend that roughly follows the surface elevation gradient (Figure 3a) oriented toward the Bras David river. Interpolated water table depths (Figure 3c) highlight three main seepage areas: (a) upstream of the knickzone, before the junction of upper stream branches, (b) just downstream of the knickzone, and (c) near the outlet. The depth to the top of the water table also increases beneath ridges, upstream and in the knickzone (i.e., between profiles P5 and P7). The estimated levels were compared to those observed in the piezometric wells installed in the catchment (Guérin et al., 2019). The observed water table levels were averaged over periods of 1 month centered on both field campaigns. These averaged piezometric levels show good agreement with water table levels interpolated with the average kriging approach (mean absolute percentage error of 29%) (Figure 3e). In comparison, water table levels interpolated solely with elevations or depths show significant deviation from levels observed in piezometric wells (errors of 120% and 57%, respectively) ( Figure  S11 in Supporting Information S1), confirming the robustness of the average kriging procedure.

Implications for Critical Zone Models
The proposed workflow allows us to jointly inform the subsurface structure and water table position at the catchment scale, pointing out that both these parameters are significantly impacted by the knickzone. In upland areas (i.e., southwestern parts of the catchment), far from the knickzone, the water table seems to be classically controlled by the topography of the hillslope organization (Haitjema & Mitchell-Bruker, 2005). There, the water table is the deepest under the hill crest, at the water divide, and then connects with the land surface in the adjacent stream (Figure 3c). The weathering front follows a similar organization, reaching greater depths under the hill crest than at the river ridges (Figure 3b), although always below the water table (Figure 3d). On the contrary, just upstream of, and in the knickzone (i.e., between P5 and P8), the weathering front and the water table appear to deepen, the latter remaining disconnected from the surface. Downstream of the knickzone, in lowland areas (i.e., northeastern parts of the catchment), the water table outcrops in vast saturated areas most likely caused by the resurgence of deep groundwater flowpaths coming from upstream of the knickzone (Figure 3c). In these areas, the regolith becomes thinner, between 0 and 5 m. There, a shallow water table is likely to prevent deep infiltration and rather favor shallow subsurface flowpaths, thus hindering regolith development Weiler et al., 2006). Overall, the difference between the water table and the weathering front ( Figure 3d), appears to be spatially organized along the stream profile (i.e., following the knickzone) rather than transverse to it (i.e., classical hillslope organization).
Most CZ evolution models coupling hydrology, weathering and erosion are intrinsically 2D, considering flowpaths from the hill crest to the river (Brantley & Lebedeva, 2021;Braun et al., 2016;Harman & Cosans, 2019). Here we show that two lateral organizations are structuring the catchment: A classic one, transverse from the hill crest toward the closest stream reach, and a less studied one, longitudinal from the upland to the lowland area through the knickzone. Understanding this particular organization, its impact on groundwater circulations, chemical weathering and erosion activity, requires more advanced models to simulate the complex 3D hydrological and reactive processes. The geophysical insights from our novel workflow could ultimately help constraining these models, both on the structure (i.e., depth of the weathered zone) and dynamics (i.e., water table depth monitoring). The non-destructive geophysical characterization of the vadose zone water content also has interesting implications to study the role of this subsurface moisture reservoir as a water source for the vegetation (Rempe & Dietrich, 2018).
Furthermore, these results allowed us to draw an interpretive cross-section of the CZ architecture along the stream and across the knickzone (Figure 3f). While it clearly shows a deepening of the weathering front in the knickzone, it also suggests that the weathered-to-fractured bedrock interface roughly follows the topography across the catchment, and is thus not influenced by the knickzone. We hypothesize that these different patterns are linked with the combination of processes taking place at two different scales. At the regional scale, topographic and tectonic stresses control the opening of fractures at depth and can lead to a depth-to-bedrock interface that roughly follows topography (St. Clair et al., 2015). Alternatively, local scale features such as hillslope orientation (Befus et al., 2011), vegetation cover (Oeser et al., 2018) or knickpoints (Comas et al., 2019) can affect weathering processes in such a way that the weathering front is locally shifted. Validating such hypothesis will require additional observation of these patterns across other catchments in different geological and climatic contexts.

Limitations and Challenges of the Workflow
A close examination of the petrophysical inversion residuals reveals that velocity distributions observed along P6 and P9 are not well reproduced at depth ( Figure S6 in Supporting Information S1). This discrepancy is most likely due to the inability of the Hertz-Mindlin petrophysical model to correctly represent seismic velocities in high-velocity and less-weathered materials (Flinchum, Holbrook, Grana, et al., 2018). Yet, this does not impact our water table level estimates, as these are located above the poorly resolved areas. Another limitation of the modelling approach is the relative lack of studies describing the elastic parameters of secondary minerals that constitute the regolith profile in the Quiock catchment. For instance, we were only able to find a single publication reporting elastic parameters for iron oxides and hydroxides (Chicot et al., 2011). In future works, we suggest to collect regolith samples at multiple locations along the seismic lines to estimate the bulk elastic parameters of the dry frame along with porosity and water content (Heap et al., 2021), so as to further constrain the petrophysical model.
The quality and robustness of the interpolation is sensitive to the number and density of geophysical data points collected throughout the catchment. In this study we were only able to record five seismic profiles with data points unevenly distributed. As a result, the elevations of the weathering front and the water table estimated in the outermost parts of the catchment are rather extrapolated than interpolated and thus only follow the main elevation trends and depth patterns. We assume that the general trend of the elevation gradient used to extrapolate both weathering front and water table levels only remains valid within the catchment, and thus do not display the kriging results beyond the watershed (Figure 3). Improving the density and coverage of geophysical data points across the catchment remains challenging in rugged and densely vegetated landscapes. Drilling additional piezometric wells is costly and strictly regulated by Guadeloupe National Park policy. Deploying extra seismic profiles across the catchment would help filling those gaps, yet requiring an improved methodology to optimize acquisition time and spatial coverage. In the current state, it takes about one or two days to collect a single seismic profile with the field conditions experienced in the Quiock watershed, mainly due to the heavy weight of the seismic equipment (∼300 kg). Recent technological advances have brought light-weight, autonomous and low-cost seismic sensors (nodes) on the market (Wilcox, 2020), opening new avenues for the seismic characterization of the CZ architecture.

Conclusions
We present a novel workflow combining seismic measurements, petrophysical modelling and geostatistical analysis, that aims at providing an extensive characterization of the CZ architecture at the catchment scale. This workflow relies on the use of a single geophysical method, so as to reduce equipment load in the field and ease deployment across an entire catchment, while still enabling the estimation of multiple geophysical parameters in order to jointly inform both the subsurface structure (i.e., depth of the main interfaces) and the water table level. This study specifically advocates for the use of seismic methods as they allow for the joint estimation of both pressure and shear wave velocities (V P and V S , respectively). These velocities are combined into a petrophysical inversion framework to estimate subsurface saturation values and highlight the position of the water table, then used to extract the depths of the weathering front and the fractured bedrock. Spatial variations of these parameters are then inferred across the entire catchment using a kriging interpolation. We apply this approach in a forested watershed representative of tropical volcanic landscapes in the island of Basse-Terre (Guadeloupe, France). The estimated water table levels are consistent with theoretical predictions and field observations. Both the weathering front and the water table appear to be impacted by a knickpoint located half-way along the stream, and deepen in its vicinity. Conversely, the top of the fractured bedrock, when shallow enough to be detected, seems to remain parallel to the topography at depths of about 45 m. These different patterns could be the results of competing weathering processes taking place at different scales. This integrated view of the CZ architecture also highlights two main spatial weathering patterns across the catchment: one transverse, along the hillslope, and one longitudinal, along the stream, strongly impacted by the knickzone. These findings call for going beyond "simple" hillslope representations of the CZ when studying hydrological and weathering processes in such complex environments.

Data Availability Statement
The seismic data used in this study are available on the H+ hydrogeophysical database which stores the geophysical data collected on the critical zone observatories of the OZCAR network. The data set can be accessed via https://doi.org/10.26169/hplus.obsera_seismic_data_2019. for lending part of the geophysical equipment, and E. Lajeunesse and the Obsera staff for running and maintaining the observatory's instrumentation and database. The authors are also grateful to A. Arènes, L. Derry, J. Druhan, M. Eichner, N. Fernandez, C. Le Traon and F. Levy for invaluable assistance during field work. The authors give thanks to D. Grana, C. Riebe and E. Gayer for the fruitful conversations which considerably helped improving this manuscript. The authors finally acknowledge A. Parsekian and one anonymous colleague for their insightful comments during peer-review. This work was supported by the CRITEX ANR-11-EQPX-0011 project and NSF EAR Award No. 1251952.