A parametric sensitivity analysis for prioritizing regolith knowledge needs for modeling water transfers in the West African critical zone

Hard rock aquifers (HRAs) in West Africa (WA) are located within a thick regolith layer. The representation of thick tropical regolith in integrated hydrological models lacks consensus on aquifer geometries and parameter ranges. Our main objective was to determine the knowledge requirements on saturated hydraulic conductivity (Ks) to model the critical zone (CZ) of HRAs in WA. A parametric sensitivity analysis with a focus on the representation of the Ks heterogeneity of the regolith was conducted with a critical zone model (Parflow‐CLM [Community Land Model]) of the Upper Ouémé catchment in Benin (14,000 km2) at a 1‐ × 1‐km2 resolution. The impact of parameter changes in the near subsurface (0.3‐to‐5‐m depth) and in the deeper regolith aquifer (24‐ and 48‐m maximum depth) was assessed in five modeling experiments. Streamflow was largely dependent on Ks and on clay distribution in the near subsurface and less on the properties of the deeper subsurface. Groundwater table depths and amplitudes were controlled by vegetation and topography as observed on instrumented hillslopes and for Ks within the literature range. Experiments with higher Ks suggested a Ks threshold where dynamics become less determined by one‐dimensional vertical and more determined by lateral processes. Such heterogeneity impacts from smaller scales need to be accounted for when hydrological models are upscaled to larger domains (1‐ × 1‐km2 resolution or coarser). Our findings highlight the need for a new conceptual approach to represent clay distribution in order to develop catchment‐scale CZ models of HRAs in WA that capture the observed processes.


INTRODUCTION
A large part (∼40%) of sub-Saharan Africa is underlain by crystalline hard rock (Vouillamoz, Lawson, et al., 2015;Vouillamoz, Tossa, et al., 2015) where 46% of the population lives. (BRGM, 2020;CIESIN, 2018). The crystalline hard rock is overlain by a thick weathering layer, called regolith Lachassagne et al., 2011). This geological structure has a large impact on the hydrological cycle and water access in sub-Saharan Africa, especially in areas with sub-humid climate (Cuthbert et al., 2019;Kotchoni et al., 2019;MacDonald et al., 2021). Recent research shows, that groundwater (GW) recharge is likely to be high in these regions, whereas storage is low (Cuthbert et al., 2019;Mac-Donald et al., 2012 but significant and provides drinking water to the population (Vouillamoz, Lawson, et al., 2015). However, borehole yields remain low Lawson, 2019), which makes the region vulnerable to both global changes (climate and the fast observed land cover change; Cotillon & Tappan, 2016) and increased demand. About 40% of the population in sub-Saharan Africa has no access to clean and safe drinking water (WHO/UNICEF, 2017), and the population is increasing at the fastest rate in the world (FAO, 2018). Reaching out a wider and sustainable access to drinking water calls for urgently improving the understanding of the characteristics of regolith aquifers, and their connectivity to soil moisture, vegetation, and surface water (SW). Integrated, physical critical zone models (CZMs) represent the interactions between the components of the water cycle (water table [WT], evapotranspiration (ET), soil moisture, and streamflow [Q]) and are therefore groundbreaking for the design of future operational models supporting integrated water resource management. The multiscale heterogeneity of the regolith body requires a distinct representation of the subsurface in CZMs, which is a challenge that hinders CZM development in data-scarce regions such as West Africa (WA). Availability of data on subsurface properties in WA remains sparse and limited to hardly available nationwide datasets (drill logs, soil profiles), and sporadic research data (Alle et al., 2018;Giertz & Diekkrüger, 2003;Kamagaté et al., 2007;Legchenko et al., 2002;Peugeot et al., 2011;Tirogo et al., 2016;Vouillamoz et al., 2014Vouillamoz, Tossa, et al., 2015;Wubda et al., 2017). For permeability, there are currently only a few datasets available for larger scales (e.g., Global Hydrogeology Maps [GLHYMPSE], Gleeson et al., 2014). Furthermore, current global datasets for depth to bedrock (Hengl et al., 2014;Pelletier et al., 2016;Shangguan et al., 2016) represent the regolith body as a continuum and generally extend to a few meters below the surface, whereas hydrogeological observations revealed that regolith aquifers in WA are stratified and can extend to several tenths of meters in the subsurface Lachassagne et al., 2011). Also, the most recent mapping approaches of recharge in Africa (MacDonald et al., 2021) and of global GW well data (Jasechko & Perrone, 2021) that we are aware of do not take into account the regolith structure. All in all, global-scale estimates of regolith thickness are not detailed enough to be accurate and relevant at the catchment scale (Alley & Konikow, 2015;de Graaf et al., 2020;Foster & Maxwell, 2018). Therefore, the representation of the spatial distribution of aquifer parameters at catchment scale remains sparse and global datasets of aquifer parameters are not detailed enough to be applied for modeling at the catchment scale in WA. Until now, there has been a lack of consensus on how to reduce the uncertainty introduced by data scarcity in integrated models regarding regolith properties. In addition, it is not yet clear on which part of the subsurface parameterization efforts should be focused due to a lack of knowledge on sensitivity of CZM towards subsurface parameters. Hence, it remains unclear how to design conceptual schemes of the regolith body for CZM applications.
In this study, we adopt a parametric sensitivity analysis framework focused on saturated hydraulic conductivity (K s ) to prioritize the knowledge needs regarding K s distributions (both vertical and lateral), within a CZM development perspective. We use the integrated hydrodynamic model ParFlow-CLM (where CLM stands for Community Land Model; Maxwell et al., 2016;PARFLOW-CLM, 2018) as a CZM on the Upper Oueme River in Benin (total area of 14,000 km 2 ), a representative hard-rock catchment for the Soudanian climate region in WA. For the regolith layer, we compile the range of values of the hydrodynamic properties of the regolith (K s and aquifer thickness) available from existing datasets in WA based on a literature review. The model performance was first evaluated by comparing simulated values to observations of streamflow, ET, and GW table depths and fluctuations from the AMMA-CATCH

Core Ideas
• Simulated water balance components for a catchment in West Africa were confirmed by observations. • Subsurface (0.3-to-5-m depth) exerts stronger control on streamflow than deeper regolith. • K s magnitude determines transition from topography-to recharge-controlled water table dynamics. • We identified a limit of K s where the impact of one-dimensional processes on water table dynamics ceases. • A high-permeability fissured zone at the bottom shows little impact on the simulations.
(African Monsoon Multidisciplinary Analysis-Couplage de l'Atmosphère Tropicale et du Cycle eco-Hydrologique) observatory . Model sensitivity was tested by running the model with different parameter settings for Ks, aquifer depths for a bimodal aquifer with a deep weathered layer, and a unimodal regolith aquifer. The amount of experiments was dominantly limited by the high computation time. The results of the five model experiments were evaluated focusing on the impact of different hydrogeological parameter settings for the near subsurface and deeper aquifer on simulated daily streamflow and WT depth (WTD) and amplitude (WTA) and the water balance components. The experiment results show the sensitivity of estimated streamflow, spatial distribution of WTD, and WTA spatial distribution to changes of K s distribution and layer depth.

Study area
The study region is the Upper Ouémé Catchment in Benin which covers a surface of about 14,000 km 2 (8.993-10.001˚N, 1.320−3.594˚E). It is representative for the conditions in the Sudanian zone (Figure 1a), an ecoclimatic zone characterized by subtropical climate with humid conditions in the south and semiarid conditions towards the north (mean annual precipitation between 700-1,400 mm (USGS, 2020a(USGS, , 2020b. A unimodal rainy season between April and October and a dry season of 4-6 mo between November and March is typical for the case study region. In the Upper Ouémé, mean temperature is 25˚C, average rainfall is 1,200 mm yr −1 and mean potential ET is 1,500 mm yr −1 . This area encompasses mostly ferrugineous soils, widespread in WA (Figure 1a, Supplemental Figure S2). Elevations vary roughly between 200 and 700 m with raising elevations in the western and eastern parts of the basin. Present rivers turn around inselbergs and hills covered with ancient duricrust. An overview on the location of the study area, vegetation classes, and the measurement stations used in this study is given in Figure 1b. In this study, we used streamflow data from two subcatchments at Beterou (10,131 km 2 ) and Cote_238 (3,126 km 2 ) located close to the catchment outlet.
Hydrological processes in this catchment have been studied using observation data and virtual model experiments for several years (Giertz, 2004;Giertz & Diekkrüger, 2003;Hector et al., 2015Hector et al., , 2018. Characteristic for the study region is direct GW recharge via percolation of precipitation, in the humid season only (Cuthbert et al., 2019;Kotchoni et al., 2019). The water budget of the Ouémé catchment reacts strongly to interannual climate variability (Giertz et al., 2006;Giertz & Diekkrüger, 2003;Peugeot et al., 2011;. Evapotranspiration is the most important discharge mechanism (75−90% of annual rainfall). It lowers the permanent WT in dry season through deep-rooted trees and riparian forests (Getirana et al., 2017;Richard et al., 2013). Streamflow is the second discharge component (10-15% of annual rainfall). Streams are ephemeral and fed by perched aquifers discharging into inland valleys, which are seasonally waterlogged headlands (Hector et al., 2018;. For the Donga catchment (586 km 2 ) situated in the Upper Ouémé,  found that baseflow from perched aquifers is the dominant component of streamflow (GW discharge from perched aquifers amounts to 10-12% of annual rainfall between [2003][2004][2005][2006]. The proportion of herbaceous-like areas in our vegetation dataset for both catchments Beterou and Cote_238 is relatively similar with a difference of only 5% (for Beterou about 17.6% against 23.6% for Cote_238). The landcover map used in this study is derived from (USGS, 2020b) for the year 2000.

Model setup
Parflow-CLM is a numerical model that has been designed to model hydrological processes in the critical zone in three spatial dimensions (Ashby & Falgout, 1996;J. E. Jones & Woodward, 2001;Kollet & Maxwell, 2006;Kuffour et al., 2020;Maxwell, 2013). For variably saturated and unsaturated flow, the pressure in the Richard's equation is solved for each grid cell in three dimensions. Subsurface properties conform to the van Genuchten approach (van Genuchten, 1980). Subsurface geometry is adjusted by a terrain-following grid (Maxwell, 2013;Maxwell et al., 2016). Any surface cell with a positive pressure head is classified as an overland flow cell (no routing scheme). This overland flow is accumulated within the catchment flowing to the lowest neighbor according to the input slopes and using the kinematic wave equation (Kollet & Maxwell, 2006;Maxwell, 2013) . Richards and kinematic wave equations are fully coupled when solving water transfers. Parflow is fully coupled to CLM, the latter solving the energy balance at the surface. The CLM is composed of a vegetation layer, in which each surface grid cell is given a particular vegetation class. It is coupled through the user-defined top cells (six in this study) from the surface to allow for temperature exchanges and root water uptakes. Latent heat fluxes are calculated by CLM and applied as a sink in the Richards equation, according to a root density function.
In this study, we used forcings from the years 2005-2006. Rainfall data are derived from rain gauge data interpo-lated with a best Lagrangian krigged method (Vischel et al., 2011). Other meteorological forcings (air temperature, specific humidity, surface pressure, and wind speed) included are from the European Centre for Medium-Range Weather Forecasts (ECMWF). Values for long-wave and short-wave radiation are satellite based data with a 30-min timestep and a 0.05˚resolution (LAND surface analysis Satellite Application Facility [LAND-SAF]) from the ALMIP2 experiment (ALMIP2, 2018; Boone et al., 2009;Getirana et al., 2017). Vegetation parameterization input includes leaf area index (LAI), stem area index, displacement heights, and roughness lengths for the two vegetation classes as in Hector et al. (2018). An exponential shape for the cumulative root fraction is assumed for root distribution of the herbaceous cover. For the deep-rooted Sudanian trees, the distribution is assumed uniform for 6 m below the surface. The model is forced at a half-hour time step. The land cover was designed based on the land cover map produced for the whole WA at a 2-km resolution for the year 2000 by the USGS database (USGS, 2020b). To simplify, we merged land cover classes associated with woody cover (i.e., woodlands, gallery and riparian forests, forests, and savannah, which represented respectively 3.6, 6.7, 0.6, and 65.9% of the surface) as a single class named "trees," and the remaining as "herbaceous," which was mostly composed of cultivated areas (22.6% of the surface). The same vegetation properties as in Hector et al. (2018) were then assigned to each of these classes.
As an initial condition, the WT was set to 8 m below the ground surface with hydrostatic conditions everywhere. This was found to be a reasonable average depth to minimize spin-up time after early experiments. We applied a no-flow boundary condition at the lateral borders of the domain. The Ouémé ParFlow model consists in a 1-× 1-km 2 horizontal resolution. The ParFlow-CLM model was run for the years 2005 and 2006. In the spin-up phase, the model was run at least five times using the input data over 2005-2006 until a water balance equilibrium was reached and trends and dynamics in WTDs and overland flow could be subtracted. An equilibrium was considered reached when the annual storage change fell below 1.2% of the incoming annual precipitation and WTs looked stable.

Model experimental design
One 2-yr model run using 32 cores needed about 5 h of wallclock time (about 160 CPU hours, ∼748 g CO 2 equivalent) on the Dahu machine of the Grenoble Alpes Recherche Infrastructure de Calcul Intensif et de Données (GRICAD) cluster of the Université Grenoble-Alpes (Berthoud et al., 2020). Limited by the enormous numerical efforts and data storage necessary for a parameter estimation with a 10-or 100-fold amount of runs per parameter a parametric sensitivity analysis was preferred (Foster & Maxwell, 2018). We test the sensitivity of the simulated hydrological variables to different parameter values and spatial distribution for K s (being one of the most important parameters controlling water transfers). The variability of K s was assumed to appear on a logarithmic scale.
We designed the experiments according to the current known hydrogeological conceptualization of hard rock aquifers (HRAs), which is built on the structure of the weathering profile of the regolith body (Bianchi et al., 2020;Chilton & Foster, 1995;Dewandel et al., 2011;Lachassagne et al., 2011) and on a literature review (see Supplemental Sections S1, S2, and S3). The profile of the subsurface can be divided into four compartments with different properties and hydrodynamic behaviors. In the upper soil, a relatively thin layer of ferruginous soils with lateritic material is found at the surface in places where it has not been cleared by erosion (Table 1, H1) The surface horizon also corresponds to the root zone of the annual herbaceous layer.
The heterogeneous content of clay in the near subsurface (Table 1, H2) leads to local decrease of K s (details on pedology can be found in Supplemental Section S2). Below appears the saprolite, which is composed of porous, unconsolidatedlike weathering material with limited hydraulic conductivity (Table 1, H3) and starting from 20-to-30-m depth a transition zone towards the crystalline matrix called saprock or fissured layer (Table 1, H4) is characteristic Lachassagne et al., 2009).
The different Horizons and their respective depth are shown in Table 1. The grid cells thicknesses D z for the horizons was refined near interfaces to better represent the pressure gradients (particularly between H1 and H2).
The heterogeneity of the regolith aquifers in the deeper subsurface in WA at the regional scale, is twofold: (a) bimodal vertical heterogeneity with an unconsolidated weathering zone overlying a fissured zone and (b) lateral heterogeneity controlled by substratum features and weathering history Lachassagne et al., 2011). In order to represent this heterogeneity, we used different subsurface configurations. From a traditional point of view on modeling, the representation of the structure of the deeper aquifer would be of greater importance for a hydrogeologist, whereas the representation of the weathering surface would be for a soil scientist or hydrologist.
With the model setup as shown in Table 1, model experiments using different parameter settings for K s were carried out (see Table 2). In total five sensitivity experiments with four model runs each were designed, of which three experiments had K s variations in the aquifer layer (H3 and H4) (Exp. 1, 2, and 3) and two experiments had K S variations at the lower border towards the aquifer (H2) (Exp. 4 and 5). Orders of magnitude of K s for the aquifer layer (H3, H4) were chosen according to a literature review (see Supplemental Section S1). The values are in the same order of magnitude as in (Bianchi et al., 2020). The different experiments are summarized in Table 2 and described as ▪ Experiment 1: This experiment was conducted in order to explore the sensitivity of the model to different constant values of K s in H3. This experiment uses scenarios for constant K s with the geometric mean, upper limits, and lower limits of the geometric mean K s from the literature review and one extra with 6 × 10 −1 m h −1 to explore higher K s . The latter allows us to assess if we miss local heterogeneities (e.g., large fractures) at our model resolution (because K s usually increases with control volume). ▪ Experiment 2: To determine the effect of horizontal heterogeneity of K s considering a single aquifer layer a stochastic dispersion of K s based on the values tested in Exp. 1 was assumed. The stochastical distribution of K s is calculated based on the turning bands algorithm (Maxwell et al., 2016) given the correlation length and the SD as in (Maxwell & Kollet, 2008). Conform with the grid resolution and domain size, correlation lengths from 1-50 km were used. The SD was approximated as the mean of the upper and lower limits of the K s geometric mean [log( smax ) + log( smin )]∕2, as no values were available from observations. ▪ Experiment 3: The influence of vertical heterogeneity in modeling was tested with a bimodal aquifer representation (saprolite layer H3+ fissured layer H4). The saprolite (ZA) refers to the upper layer of the aquifer that is important for water storage. The setup with ZA and the fissured layer (ZFA) consists of four horizons and nine vertical grid cells (Table 1). Horizons H1, H2, and H3 are similar in thickness (the first seven cells in Table 1). Total depth is z = 48 m. The saprock is placed in the bottom layer H4 and limited to 24 m in depth. K sH4 was varied using a factor of amplification of 10 and 100 in comparison with the benchmark K s value in H3. These scenarios were chosen in order to represent the vertical profile of transmissivities according to the hydrogeological conceptual model (Lachassagne et al., 2009. Transmissivity is assumed to be higher in the fissured layer relative to the saprolite layer. ▪ Experiment 4: In the near-subsurface layer (H2), a range of constant values of K s was used to explore the sensitivity of the model to different values of K s in H2. After some preliminary tests, the K s range has been narrowed around low K s values where the model has shown to be highly sensitive. ▪ Experiment 5: The lateral redistribution of water in or above the near subsurface layer has been shown to be highly dependent on the extent of clayey areas or inland valleys (Giertz et al., 2006;Hector et al., 2015Hector et al., , 2018Masiyandima et al., 2003). According to the abovementioned literature review, we study the model sensitivity to the extent of a low-permeability layer in the subsurface, hypothesizing that its extension is a function of local topography. We differentiated cells closest to and farthest from the river according to a drainage area map and a threshold. The cells closest to the river were supposed to behave like inland valleys with low permeability in H2 corresponding to clay, whereas the cells furthest from the river were given a higher permeability in H2. Cells assigned as clay cells with K s values corresponding to 1 x 10 −4 m h −1 have been mapped with four different scenarios: no clay cells at all or cells spread upstream from the lowlands, using a drainage area map: cells draining more than 2, 8, and 72 cells and corresponding to 48, 24, and 10%, respectively, of the total number of cells have been assigned as clay cells.
For instance, the subsurface clay extent of 10 and 48% are shown in Figure 1b.
One benchmark configuration (Exp. 1b) has been considered in this study, with the saprolite only (with a total depth of z = 24 m), and no fissured zone (no H4). It consists of three horizons and nine vertical grid cells with variable thicknesses (Table 1) and is used in all modeling experiments except Exp. 3. The benchmark also considers near subsurface clay mapping in H2 using 24% of the lowest cells according to the drainage area map (and a threshold of eight cells being drained) associated with clay properties. This has been found to be the simulation closest to the real conditions according to our knowledge on K s and clay distribution and with the best model performance (with regard to streamflow, ET, and WTD and WTA).

Experiment metrics and evaluation data
Streamflow, GW level and ET data of the catchment area used for validation of the simulation results can be downloaded from the AMMA-CATCH observatory (AMMA- CATCH, 1990CATCH, , 1996aCATCH, , 1996bCATCH, , 2003CATCH, , 2005. The simulations were first compared with observations from measurement stations Cote_238 and Beterou (Figure 1b). The objective was to identify if the model was able to reproduce observations, given the choice of realistic parameters according to the literature review. The Kling-Gupta efficiency (KGE) criterion (Gupta et al., 2009) was used to evaluate the fit of modeled to observed streamflow hydrographs. A KGE of 1 implies a perfect fit. 24% of clay cells Impact of varying clay content in the near surface below river cells

48% of clay
Note. Experiments 1-3 were conducted with a focus on properties of the aquifer layer H3 and H4 (in bold); Experiments 4 and 5 were conducted with a focus on properties of the near subsurface layer H2 (in italics Evapotranspiration was additionally compared with observations from a large aperture scintillometer in the Ara catchment (Guyot et al., 2012) and against measurements from an eddy covariance flux tower at a study site covered by herbaceous vegetation in Nalohou (herbaceous hillslope in the following) and covered by trees in Bellefoungou (woody hillslope in the following. Locations are shown in Figures 1c and 1d (Mamadou et al., 2016). The research piezometers data used in this study are also from the two study sites of Nalohou and Bellefoungou (Figures 1c and 1d), which respectively represent an herbaceous and a woody hillslope.
A comparison of simulated and observed GW states and behavior raises several challenges, such as scale issues (local observations vs. 1-km × 1-km average of simulated WT) and withdrawal affected wells used to observe water resources. We therefore evaluate simulation results in their ability to represent elementary GW redistribution processes at the hillslope scale, which seems to be the representative elementary volume for GW transfers in the area (Richard et al., 2013). This is achieved by comparing the relationship between mean WTD and yearly WTA (difference between the yearly maximum and the yearly minimum) observed along both a woody ( Figure 1c) and a herbaceous (Figure 1d) hillslope and simulated for model cells with corresponding vegetation. Hereby, aquifer properties (i.e., specific yield [Sy], K s ) are supposed to remain homogeneous. The observations represent heterogeneity at the scale of few hundreds of meters, whereas the simulations represent heterogeneity at a kilometric resolution, but both hold the fundamental hillslope-scale redistribution processes. The rationale behind the WTD vs. WTA relationship is that it carries both the steady state behavior of the WT through the distribution of mean WTD, and the transient behavior (both seasonal recharge and discharge) through the distribution of yearly WTA.
To assess the sensitivity of the model to the different scenarios with different hydrogeologic configurations, the benchmark case was used as a reference for further evaluation of the simulated streamflow and GW table of the different experiments. The experiment results were compared regarding simulated streamflow phase and onset visible in the hydrographs. Additionally, we intercompared experiments by evaluating variations in proxies of GW redistribution processesnamely, WTD and WTA-using maps. Groundwater table depth is an important indicator for the productivity of an aquifer. The closer (the deeper) the GW table is with respect to the surface, the thicker (thinner) the saturated thickness of the aquifer, the higher (lower) aquifer productivity and sustainability (Bianchi et al., 2020;Vouillamoz et al., 2014Vouillamoz, Tossa, et al., 2015).
Furthermore, the yearly water budget gives an overview on the water fluxes and their representative quantities in the watershed. Based on streamflow, precipitation (P), ET, and water storage change (WSC), the water budget is calculated using Equation 1. Water storage change is the total variation of water storage in the soil water column and surface ponding water for the simulated time period. Q in is zero (no-flow boundary conditions) and Q out is the simulated streamflow at the surface: Water budget differences are evaluated to intercompare simulation results and to find out what terms are mistaken or similar from one simulation to another.

Comparison with observations
The  Figure 3 shows daily observed and simulated ET for herbaceous and woody savannah vegetation (Figure 3a and 3b, respectively) for available observations (eddy covariance and a large aperture scintillometer data) and simulated herbaceous and woody savannah pixels. This allows an evaluation of simulated ET with respect to different vegetation context. Measurements for 2005 were not available, but we added further available data, which allowed appreciating annual cycles for both vegetation and the interannual variability. Note that the scintillometer path covered a mix of crops and woody savannah (Guyot et al., 2009). It was plotted on purpose in both F I G U R E 3 Evapotranspiration (ET) simulation results for the Benchmark Case 1b for grid cells corresponding to trees and herbaceous cover and observations from large aperture scintillometer (LAS) data in the Ara watershed, which captures trees and herbaceous cover and from eddy covariance towers (EC) at Nalohou (herbaceous) and Belefoungou (trees). Simulated ET is the daily spatial mean and the 1st and 99th percentile of the spatial distribution give the envelope. Sim, simulated; Obs, observed amplitude. In December and January ET reaches a minimum (close to zero in the dry season) and remains high during spring and summer months until it reaches a minimum again in July. The simulation also reproduces well the lowering of daily ET in the core of the rainy season (August) associated with a higher cloud cover, as there is no water limitation during the period neither LAI lowering. Simulated ET for tree cover shows more spatial variability (shaded blue) in the dry season than herbaceous cover. For trees, ET does not entirely reach zero in the dry season (∼0.5 mm d −1 ), as some trees keep leaves and continue transpiring (Mamadou et al., 2016). For this period, the LAI is still high as trees renew their leaves along a 1-mo period in the dry season even without rain. This variability in the model is then associated with the root distribution and the water availability and WTD.
We additionally compared simulation results for GW WTD and WTA against observations from a woody and an herbaceous hillslope (Figure 4).Hereby, the hypothesis was made that Sy remains similar at the hillslope scale so the observed WTA change with WTD can effectively be compared with the model (with constant Sy). The different colors correspond to the different vegetation type (trees in green, herbaceous in red). In the WTD vs. WTA distribution, an obvious separation of herbaceous and tree cover appears for both observation and simulation outputs. Water table depths are lower below herbaceous pixels (which simulate less transpiration) and higher below pixels with woody savannah. The largest WTA under tree cover corresponds to pixels where WTD is the closest to the surface (2.5-5 m). For WTD between 5 and 10 m, WTA is about 3.5 m with a large variability. Then, simulated WTA drops down below 1 m for WTD between 10 and 12 m. This threshold is related to root depth, which is fixed to 6 m below the surface in all simulations. The low WTA for deep WTD is consistent with observations from piezometers.
Precipitation was 1,105 mm in 2005 and 1,034 mm in 2006. The evaluation of the water budget shows that ET was 93.5%, streamflow was 6.5%, and water storage changes were 3% of precipitation, on average, for the simulated years (2005 and 2006).

Experiment results
The effect of K s variation in H2 and H3/H4 (for each experiment mentioned in Table 2) was tested, and sensitivity on streamflow near the catchment outlet, on the distribution of GW WTA over the catchment area, and on the spatial distribution of the GW WTD was evaluated. Only 2006 results are shown in Figures 4-8 to have a lighter presentation of the analysis, as results and comments for both simulated years are identical.

Streamflow
Streamflow hydrographs for the locations of Beterou for each experiment are shown in Figure 5. For Exp. 1 (Figure 5, Exp. 1a-d), the impact of hydraulic conductivity is low except when hydraulic conductivity was set 10 times higher than the upper bound of the explored range of K s and 40 times higher than the reference case ( Figure 5, Exp. 1d). For this simulation, the streamflow onset happened earlier and streamflow was more reactive at the beginning, peaks slightly higher in comparison to the benchmark simulation ( Figure 5, Exp. 1b). When adding stochastic spatial heterogeneity in hydraulic conductivity in the aquifer level (H3), no sensitivity was observed ( Figure 5, Exp. 2a-d). Considering a two-layer aquifer (H3 + H4) did not change the streamflow pattern ( Figure 5, Exp. 3a-d) either. When K s in the ZFA (H4) is set to 100 times K s in ZA (H3) (Figure 5, Exp. 3d), the beginning of streamflow occurred earlier as for Exp. 1. The overestimation of streamflow with regard to the benchmark simulation was even stronger than in Exp. 1d probably because of the increased thickness and therefore transmissivity of the aquifer in Exp. 3d. However, overestimation only occurred in the beginning of the streamflow period. Experiment 4 for varying constant K s of H2 showed more sensitivity than the preceding experiment ( Figure 5, Exp. 4a-d). Indeed, tested Ks range is much less than in Exp. 1-3 ( Figure 5, Exp. 1-3). The lower the Ks in H2, the higher the streamflow. A strong sensitivity is also observed for different percentages of clay below rivers (Figure 5, Exp. 4a-d) in terms of reactivity (the higher the clayed area the higher the peaks) but also in terms  Table 2. The annual mean WTD for 2006 is shown of streamflow departure (the higher the clayed area, the earlier the departure).
Overall, a comparison of the modeled and simulated streamflow hydrographs at the station of Beterou reveals that there are only minor differences with regard to the hydrogeological configuration chosen. The high K s chosen in 1d and 4d probably correspond to the limit of permeability where aquifer properties start to control streamflow generation. Below this limit, the differences in simulated streamflow hydrographs for Beterou and Cote_238 are likely to be more controlled by H2 than by H3 and H4, as shown in Exp. 4a-d and 5a-d.
We evaluate the model performance for streamflow using the KGE scores with reference to the streamflow observa-tions The best model scores (for the two simulated years) were obtained for Beterou for the model runs 3d (KGE = .90), 4b (KGE = .74), and 5c (KGE = .81). At Cote_238, the best scores were calculated for 3d (KGE = .41), 4a (KGE = .37), and 5c (KGE = .58). The r 2 value shows that the temporal correlation is good in most cases as it is >.7 for all model runs. For a more detailed analysis of the three terms of the KGE we used polar plots to intercompare the scores of all simulation (Schwemmle, 2020). The polar plots in Figure 5 show that the largest differences of the KGE is obtained for Exp. 5 and 4. The variability of KGE is the highest for the different runs in Exp. 4 and 5 and the lowest for Exp. 1 and 2, especially for bias and variability errors. In general, the mean is likely to be F I G U R E 7 Simulated water table (WT) amplitude maps for the northern Ouémé for each experiment listed in Table 2. The annual mean water table amplitude for 2006 was calculated underestimated. The correlation term is in the same range for most of the runs with the largest deviation of the correlation term for 5a and 4d.

Groundwater
Annual average WTD and WTA maps allow a first analysis of the spatial distribution of water quantities and fluctuations in the subsurface of the Ouémé catchment. Figure 6 presents the yearly average of the WTD. Differences in WTD are particularly visible along the river network as the WT is connected to the river network. We can also observe the vegetation pattern on almost every case. For the benchmark simulation (1b, 2a, 3a, 4a, and 5c), average WTDs are deeper below woody savannah (∼8 m) in comparison with herbaceous areas (2-5 m). The WTD maps reflect the pattern of the vegetation cover in almost all model runs in Exp. 1, 2, and 3, especially for lower K s in the aquifer (1a). Only 1d and 3d differ. With increasing K s the WT below river cells rises and equilibrates (1c). With extremely high K s = 0.6 m h −1 , the WT becomes shallower below river cells (1d).
The vegetation pattern disappears as K s increases. Instead, there is a tendency of deeper WT to occur in the western F I G U R E 8 Water table amplitude (WTA) plotted against water table depth (WTD) for all simulations colored according to the respective vegetation class parts of the catchment. For a bimodal aquifer and a 100fold increase of K s in the fissured layer, the WTD is greater than 15 m everywhere at distance of the river network (3d). The WT in cells around and below the river network connects with the surface. Different correlation lengths seem to affect the distribution of WTD along the river network but otherwise do not largely change the WTD distribution (Exp. 2).
For the model runs with K s variations in the near subsurface (Exp. 4 and 5), remarkable differences in the WTD occur, but the reaction of the WT below rivers cells is less marked than in the H3 experiment. With low K s in the near subsurface, there seems to be a disconnection of the deep WT and the rivers (4b). With increasing K s in the subsurface, the WT gets more and more connected to the surface, whereas the width of the connected areas remains relatively thin (4c and 4d). The relation of vegetation cover and WTD is also supported by histograms of the WTD distribution (see Supplemental Figure S1). In all model runs 4a-4d, the WTD reached the surface in areas with herbaceous cover. The only experiment showing fine-scale variations in distributions of the WT along the river network is Exp. 5 with varying clay contents in the near subsurface. Additionally, the effect of vegetation cover is less visible in the maps for the runs 5a-5d. For high amounts of clay (48%) below river cells (Exp. 5c), the WT is significantly deeper at larger distances of the rivers. With less clay content below rivers, variations of WTD along the river network happen in larger sections.
Water table amplitudes are displayed in Figure 7. Again, no large differences are visible in Exp. 1 and 2. In contrast, differences in WT fluctuations are particularly strong amongst the runs in Exp. 3 and 4. Experiment 3 shows, that amplitude changes are less pronounced and smoother with higher K s . In 3d, amplitude changes are more evenly distributed over the whole domain and differences in WTA of cells with different vegetation and cells in proximity of the river network are wiped out. This is also the case for very low conductivities in 4b, where the WT is disconnected from the surface. In Exp. 4, the amplitude fluctuations along the river network become particularly strong between 0.0003 and 0.0004 m h −1 . In all the other experiments, amplitude changes along the river cells remain relatively small. Below herbaceous vegetation cover, amplitude fluctuations are generally a little bit larger around 3 m. In general, where WT is deeper (3 m < WTD < 10 m), fluctuations are also greater (see benchmark simulation in 1b, 2a, 3a, 4a, and 5c). During this stage, WT becomes less and less limited by the surface or the high permeability of H1. Vegetation type shows a distinct pattern of the WTA vs. WTD relationship (similar to the pattern in Figure 4) for most cases, except for 1d and 3c and d, which are the experiments with the highest aquifer transmissivity. For woody-savannah, WTA is stronger for WTD between 3 and 10 m, whereas WTA and WTD are smaller for cells with herbaceous vegetation. For experiments with high aquifer transmissivities, herbaceous and woody savannah WTA vs. WTD relationships mix (3d). The curve flattens, while deeper WTD also occur for grid cells with herbaceous vegetation and below woody savannah deeper WTD (5-10 m) do not co-occur with increasing WTA variability anymore. A slightly higher sensitivity of simulated WTA and WTD on different orders of magnitude of K s in the near subsurface is also visible in 4b (K s = 0.0001 m h −1 ). The reason could be a larger gradient between H2 and H3. For lower K s in H2, less water percolates into the subsurface. Figures 7 and 8 denote that the GW table and the stream network do not connect. A higher amount of water flows to the stream.

3.2.3
Water budget Figure 9 displays the water budget for all experiments. Experiments 4 and 5 have the largest impact on the water budget, especially on streamflow and ET. The highest streamflow in 4b indicates that the gradient between H2 and H3 is also influencing the water budget terms. A similar effect emerges with a significant increase of the clay content in cells spreading from the river network. The least impact on the water budget is for Exp. 2.

DISCUSSION
The observed flows at stations Beterou and Cote_238 are well reproduced, as well as GW and ET behaviors despite the hypothesis of spatially homogeneous properties in our model and without model calibration. This suggests that our catchment scale model is satisfying at the first order and captures the right hydrological processes as described in prior findings (Richard et al., 2013;, with K s values consistent with observations.

Model sensitivity to the subsurface parameterization
Streamflow is extremely sensitive to H2 permeability in a narrow K s range (between 1 × 10 −4 and 4 × 10 −4 m h −1 ) where either all is infiltrated or nothing is infiltrated below H2. Experiments 4 and 5 show the same direction in KGE components changes ( Figure 5, right panel) and, while keeping high r 2 values, mostly exhibit changes in the amount of water flow pulses and on the onset of the flowing period. This suggests that representing a shallow reservoir, at least in the topographic lows (the common feature between Exp. 4 and 5) is key for streamflow simulations.

Model sensitivity to the deeper aquifer parameterization
Simulated streamflow does not seem to vary with K s in the deeper aquifer (H3) for realistic values of K s , whether a fissured layer is taken into account or not ( Figure 5). An overestimation of streamflow occurs only in the beginning of the streamflow period, in simulations where WTD reaches the surface below river cells, indicating that the WT rises to the surface earlier ( Figure 5, Exp. 4a-d). A stochastic distribution of K s in H3 (Exp. 2) had the least impact on streamflow and on the water budget, and this finding is additionally confirmed by the detailed evaluation of WTD and WTA (Figures 6, 7,  and 8). This indicates that kilometer-scale spatial variability of subsurface properties in H3 is not important for hydrological modeling with 1-km scale in the study area, and that an average K s value is enough at the first order for the entire region. This is in contradiction to findings from Atchley and Maxwell (2011), who found that spatial variability matters. However, their simulations were at much finer scales, which suggests that spatial variability has to be considered in front of represented spatial and temporal transfer scales. In our case, the observed K s variability is small compared with water transfers within a catchment size of 10 4 km 2 and seasonal time scales. In effect, subgrid scales variability could have larger impacts on water transfers through H3 to small surface drains.
Water table depth and WTA react differently for different vegetation types (Figure 8). At low K s up to a certain limit (around 10 −5 m h −1 ) in the deeper subsurface, the forests extract water from the near subsurface, which is thus not contributing to infiltration or streamflow. Regarding the relationship of WTD and WTA, a dependency on vegetation was also clearly visible and confirmed by hillslope observations. Where the WT is deep (below 10 or 15 m) the impact of trees on WT fluctuations is decreasing (Figure 8). The decrease of WTA with WTD for deep-rooted woody savannah is probably due to the buffer effect of the vadose zone where trees uptake most of the water before it reaches the WT. Water rapidly infiltrates into the subsurface first and the WT rises again once the aquifer is saturated.

Model behavior when departing from known range of K s values
Our results show that the order of magnitude of K s in the aquifer layer may lead to a wrong process representation in the model. The transition from a topography-controlled WT (where WT is a subdued replica of the surface topography) towards a recharge controlled WT when increasing the aquifer transmissivity is in agreement with Haitjema and Mitchell-Bruker (2005). Experiments 1 and 3 had a strong influence on WTD vs. WTA pattern (Figure 8). For higher K s , topography is visible on the maps (Figure 6).This leads to mixed WTD and WTA relationships for all vegetation classes (Figure 6, Exp. 3d and 1d), which is not reflected by measurement data (Figure 4). Furthermore, the simulated low WTD below river cells (Figures 6, 3d and 1d) are not consistent with geochemical investigation . Still, these model runs achieved good scores for streamflow (KGE) (Figure 5, right panel). This underlines the value of evaluating CZM on different variables of the hydrological cycle.
An increase of K s in the deeper subsurface (between 5and 48-m depth) therefore impedes the regulating effect that vegetation has on the WT (Exp. 3d). Thus, the impact of water uptake through tree roots on the WT decreases for higher transmissivities. The evaluation of WTD and WTA enables to determine a limit in K s where the influence of one-dimensional vertical processes (direct percolation, infiltration, diffusion in the vadose zone, and water uptake by roots of the trees) is likely to diminish in favor of lateral redistribution of water, controlled by local slopes or the distance to the drainage network. However, observations do not confirm a 100-fold transmissivity of the saprock in comparison with the saprolite, nor a flattened WT.
Some uncertainty remains, notably because we assume that Sy is constant everywhere in the aquifer layer H3 for evaluation of WTD and WTA. In reality, Sy in this layer is more heterogeneous and the amplitude of the WT is a function of Sy. Also, a correct validation of the simulated WTD and WTA requires other metrics and further investigation of WTD vs. WTA relationships from observations (e.g., for different vegetation types). To define these metrics, more research on upscaling of models and the effect of coarsening model resolutions on WTD and WTA is necessary. This is a crucial step to be able to reach the next level and to be able to do scenario modeling for the catchment scale in WA and to assess whether climatic changes have similar impact in comparison with catchments in the global north.

Implication for future modeling work
We designed the sensitivity analysis from a critical zone perspective, by combining approaches from different disciplines to design the experiments. Although it is a common practice to use stochastic aquifer properties distributions (Dell'Oca et al., 2020;Guadagnini et al., 2018;Tartakovsky et al., 2017) in hydrogeology, recent paradigm changes in HRA structure (Lachassagne et al., 2009(Lachassagne et al., , 2021 also calls for developing new Vadose Zone Journal conceptual models. This served as a background for Exp. 1-3. On the other hand, we adopted a traditional pedological approach by distributing clay material following the concept of catenas (Khomo et al., 2011). These choices remain highly debatable, and we hope to gather feedbacks from different communities.
Within this framework, we found the subsurface layer (H2) to be much more sensitive than the deeper aquifer (H3) within realistic K s ranges, and across all experiments. This holds true for streamflow, GW behavior, and simulated water budget at the catchment scale. This calls for developing research on the soil organization from a hydrological perspective. Clay formation is strongly connected to alteration and erosion processes. In this study, we assume that the primary control of the extension of clay upstream is the WT, which corresponds to a classical point of view in pedology. A promising path is to focus on surface redistribution processes of material through pedimentation, as recent research evidenced the major regional role . Erosion leads to a lateral redistribution of soils and weathering profiles develop successively above, leading to complex spatial distribution of clays. This could be a starting point for a new approach to clay formation mapping. The lack of subsurface mapping and inland valley distribution maps, together with the still relatively coarse spatial model resolution (1 km 2 ), did not allow for a locally detailed discretization of subsurface clays in this study. The validation and comparison of different clay distributions with global soil profile data (e.g., WoSIS; Batjes et al., 2020) could be an ensuing task.
At kilometric resolutions, vertical variations of subsurface properties within a model's column are often taken into account while lateral heterogeneity is neglected. With coarse model resolution, the use of scaling relationships has been proposed by Foster and Maxwell (2018) for the ParFlow model. With a parametric sensitivity analysis, they found a dependency of effective K s on topography for headwater catchments. In the context of the flat peneplains of the Ouémé catchment, it is likely that such topographic resolution issue is negligible, whereas alteration and erosion processes that are driving the distribution of clay content at scales smaller than the model resolution become dominant. Also, baseflow from perched aquifers and interflow is an important contributor to streamflow in the Ouémé, and the amount of lateral flows in the near subsurface are likely to be high. Future work on resolution issues should therefore focus on scaling parameters related to soil distributions.
Results also showed the important role of vegetation in controlling WT depth and dynamics. On the background of the fast transitions that WA landscapes are currently undergoing, the impact of different land cover maps, vegetation clearing, and the increase in cultivated area is of primary importance when modeling water resource availability in WA catchments. As an example, the subcatchment of Cote_238 has been severely deforested starting from 2000, as the land cover maps for the years 2000 (chosen in this study) and 2013 show (Cotillon & Tappan, 2016). We suspect that this clearing most likely already influenced observed streamflow from the years 2005 and 2006 used as validation data in this study (Figure 2). Simulating the streamflow with a vegetation map with more woody savannah areas could explain the underestimation of simulated streamflow for years 2005 and 2006, considering that less woody savannah leads to less ET and then more streamflow.
We can conclude that the fine-scale heterogeneity of K s in H2, notably the distribution of clay in this layer together with the vegetation type, is the first-order target to represent HP in CZMs for the study area (deep-weathered tropical soil overlying a hard-rock basement and ∼1,000-1,500 mm annual rainfall and P ∼ ET). This is a crucial step before using scenario modeling as a tool to understand the complex interplay among SW, GW bodies, vegetation, and subsurface properties during no flow or streamflow intermittency events (Schilling et al., 2020). Although our results show that large-scale spatial variability in K s in the saprolite (H3 and H4) is not important for modeling water transfers, it may be important on the smaller scale. Future work could therefore also be focused on local CZM, bounded by the present study findings in order to obtain estimations for well yields and GW availability for pumping.
Our model environment is representative for catchments with deep-weathered tropical soils over hard rock basement and subhumid climate, and therefore our results are likely also transferable to other catchments in WA. From a global perspective, our findings are confirmed by other studies that obtained comparable results for similar geological context in a tropical climate Ruiz et al., 2010).

CONCLUSION
Our objective was to understand the relevant degree of knowledge on subsurface properties in hard rock areas of Africa for modeling the critical zone at the catchment scale. We identified this as a research gap due to the lack of experimental data and other datasets for regolith properties, on the one hand, and due to the need for developing CZM at managed hydrosystem scales on the other hand. We conducted a parametric sensitivity analysis with five experiments and 17 model runs to investigate the role of hydrogeologic properties and geometry in modeling of the hydrological processes (GW and surface flows, as well as the influence of vegetation and the total water budget) in the critical zone in a WA catchment (14,000 km 2 ). The comparison of the reference case with field observations showed that the model performs well for simulation of streamflow, ET, and the GW table. The best results for both stream gauges (Beterou and Cote_238) were obtained assuming 24% of clay content spreading from river cells in the subsurface and a single aquifer layer with a constant K s = 0.015 m h −1 (Figures 2-4). The intercomparison of streamflow hydrographs (Figures 2 and 5) using different model scenarios for K s heterogeneity (Table 2) shows that the model is not highly sensitive to the values of hydraulic conductivity characteristic for the deeper subsurface that were explored in this study.
The emphasis in this study was not on how to achieve the best model results, though. We identified a range of K s (between 1 × 10 −4 and 4 × 10 −4 m h −1 ) in the near subsurface (H2) where modeled streamflow was highly sensitive. An evaluation of the water budget also showed that variations of K s in this layer lead to different model results. A comparison of model results with observations from two monitoring sites situated along hillslopes in different vegetation context indicates furthermore that there are important feedbacks between vegetation and the GW table. Simultaneously, there is a topographic effect, which also influences the modeling results (WTA and WTD). Clay distribution was determined depending on topography and had a significant impact on simulated streamflow and on variations of the GW table along the river network. The sensitivity analysis revealed that a stochastic distribution of K s (for K s values confirmed by literature) in the aquifer layer does not have an impact on simulated streamflow and GW table depth. Furthermore, the conceptual model used with the unconsolidated layer (H3) above the fissured aquifer layer (H4) does not significantly change model results.
In order to determine trajectories for GW and SW resources, we found proof that a correct representation of the catchment properties in CZMs requires a correct representation of the upper part of the weathering layer (clay content and vegetation) in addition to average information on hydraulic properties of the aquifer. Other critical zone modelers should consider these findings as guidance in their choice of parameterization of aquifer properties. Moreover, this is also relevant information for experts working in the field of GW prospection because varying vegetation may have an impact on saturated aquifer thickness at the local level, and therefore on the exploitable water quantity as well.

A C K N O W L E D G M E N T S
We would like to acknowledge the valuable contribution of Carolin Winter for her MSc thesis work on defining the geometric mean hydraulic conductivities, which we used as a parameter in the runs done (Winter, 2018). We also acknowledge the value of the data from the AMMA-CATCH observatory (http://www.amma-catch.org) and the data collected and summarized in the ALMIP project. We are grateful to Direction General de l'Eau in Cotonou (Benin) for sharing the discharge datasets. Part of this research was directly supported by a grant of the German Academic Exchange Service for Project Related Personal Exchange (DAAD-PPP) within the project "Development of Subsurface Parameteri-zation focusing on hardrock aquifers" assigned to I. de Graaf, and the French counterpart (MEAE and MESRI through the PHC-PROCOPE program) assigned to B. Hector. For this study, we used the tagged version 3.3.1-IGE of Parflow-CLM (PARFLOW-CLM, 2018). All (or most of) the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the Equip@Meso project (Reference ANR-10-EQPX-29-01) of the programme "Investissements d'Avenir supervised by the Agence Nationale pour la Recherche."

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.