Recession discharge from compartmentalized bedrock hillslopes

. We used numerical modelling to explore the role of vertical compartmentalization of hillslopes on groundwater flow and recession discharge. We found that when hydraulic properties are vertically compartmentalized, streamflow recession behaviour may strongly deviate from what is predicted by groundwater theory that considers the drainage of shallow reservoirs with homogeneous properties. We further identified the hillslope configurations for which the homogeneous theory derived from the Boussinesq solution approximately hold and conversely for which it fails. By comparing the modelled recession 15 discharge and the groundwater table dynamics, we identified the critical hydrogeological conditions responsible for the emergence of strong deviations. The three main controls are i) the contribution of a deep aquifer connected to the stream, ii) the heterogeneity in hydraulic properties, and iii) the slope of the interface between shallow permeable compartment and deep bedrock one with lower hydraulic properties. Our results confirm that a physical interpretation of the recession discharge exponent 𝑏 from the classical equation −𝑑𝑄/𝑑𝑡 = 𝑎𝑄 𝑏 , and its temporal progression, are information that can only be 20 interpreted with knowledge on structural configuration and heterogeneity of the aquifer. homogeneous porosity. We explore contrasts in hydraulic conductivities between both compartments according to the ratio 𝑟 𝐾 as 𝐾 𝑈 = 𝑟 𝐾 𝐾 𝐿(𝑥=0,𝑧=0) . The porosity of the upper compartment and the interface remain constant across simulations with 𝜃 𝑈 = 𝜃 𝐿(𝑥=0,𝑧=0) . as recession advances than the ones predicted by the analytical solution. We expect that for low values of 𝑟 𝐾 , 𝑏 will take values close to 1, similar to the case of a thick aquifer with only marginal transmissivity variations with head. In contrast, for a high ratio 𝑟 𝐾 , the change in transmissivity due to the fast drainage of the 170 upper, more permeable, compartment satisfies the assumption of the Boussinesq solution for a relatively thin unconfined aquifer over a large range in Q . This comparison allows us to identify the influence of the deeper compartment in shaping the recession hydrograph with respect to 𝑟 𝐾 It is urgent to further develop knowledge on the role of such vertical compartmentalization and its representation at the catchment scale to improve prediction in stream baseflow regimes. Our results allow us to identify that both slope interface between deep bedrock and fractured/weathered shallow compartment and their contrast in hydraulic properties are key parameters to account for when modelling storage-discharge functions of hillslopes. Such characterization, although demanding, is possible through geophysical investigation with combined surface and borehole techniques. While our approach assumes continuity of the 355 compartment along the hillslope, another key parameter influencing stream discharge dynamics during baseflow is related to the volume and spatial connectivity of shallow and deep compartments at the catchment scale. While mechanical and bio-geochemical processes may increase the permeability and storage capacity of the bedrock, the spatial representation of this


Introduction
Recent streamflow recession analyses have provided new insight into the spatial and temporal variability of hydrological behaviour across the globe (Fan, 2015;Roques et al., 2021;Schmidt et al., 2020;Tashie et al., 2020). These analyses have confirmed that the relationship between groundwater storage and stream discharge (hereafter called storage-discharge 25 functions) may be highly nonlinear and deviate from what would be predicted from groundwater theory under the assumption of homogeneous subsurface properties. They have also shown that such non-linearities in storage-discharge functions may be variable across different recharge events, suggesting strong dependencies to initial saturation conditions (Jachens et al., 2020;Tashie et al., 2020). In the current context where hydrological systems are increasingly impacted by more frequent and intense 1999; Saar and Manga, 2004), and 2) hydraulic properties may be enhanced in the shallower part (up to a depth of 100 meters) due to supergene weathering processes and the presence of elevated joint densities in response to unloading and topographic stresses (Dewandel et al., 2006;Harman and Cosans, 2019;Moon et al., 2017;Rempe and Dietrich, 2014;Riebe et al., 2017). 65 Besides such evidence of the vertical compartmentalization of hillslopes, efforts aimed at modelling groundwater storage variations and streamflow discharge often simplify this heterogeneity by considering a homogeneous aquifer with effective hydraulic properties. Here we aim at identifying in which context this simplification may hold and conversely in which ones it could be a source of large uncertainty and error when estimating storage-discharge functions of hillslopes.

70
Storage-discharge functions with values of b larger than 1.5 during baseflow suggest that the resistance to drainage can be even stronger than the one induced by the decrease of transmissivity with the hydraulic head. The decrease of hydraulic conductivity with depth and vertical compartmentalization plays a significant role in the emergence of anomalous values. In this perspective, Rupp and Selker derived analytical solutions for the case of horizontal (2005) and sloping (2006) aquifers with a power law decrease of hydraulic conductivity with depth such that ≅ z n , where z is the elevation and n is the power-75 law exponent ( > 0). The authors showed that b is determined by according to = 2 +3 +2 for a horizontal aquifer, such that is an increasing function of n with a lower limit of 1.5 for near zero to an upper limit of 2 for large values of n. As expected, the stratification of the hydraulic conductivity increases b and slows down the drainage. In this approach, groundwater flows are constrained to the upper-most material, implicitly supposing that there is a shallow impermeable boundary parallel to topography. Although this may be a reasonable assumption in some cases, we hypothesize in this paper that, in other contexts, 80 deeper aquifers may also significantly contribute to stream discharge (Figure 1a.). To test our hypothesis, we determine the recession discharge behaviour arising from the presence within the hillslope of a shallower compartment overlaying on a deeper one with contrasted hydraulic conductivities. This structure typically represents a permeable near-surface regolith aquifer underlain by a deeper less permeable bedrock compartment ( Figure 1a). Such structures are expected to induce the two following effects. First, the aquifer discharge is expected to slow down progressively from the fast drainage of the permeable 85 near-surface layer with low characteristic timescales (high value of a as pictured in blue in Figure 1b and c) to the drainage of the deeper less permeable compartment with consequently longer drainage timescale (lower value of a as pictured in red in Figure 1b and c). Second, the sloping interface between the two compartments conditions the shape of the water table profile and progressively modifies the partitioning of flows on either side of the interface. The transition between the trends given by each of the compartments is likely to cover a wide range of times and behaviours. For a simpler system of two parallel 90 reservoirs (Figure 1b-d), Gao et al (2017) have shown a broad variety of transitioning flows and found that b may be written as a function of the ratio of the fast flow to the total flow and of the hydraulic conductivity ratio between the two compartments (Eq. 19 of their paper and in Appendix A.1). We will use the two-linear reservoir theory as a reference to describe our results and explore in which context such simplification might hold to describe storage-discharge functions of vertically compartmented hillslopes. https://doi.org/10.5194/hess-2022-7 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.

Methodology
We investigate the time dependency of streamflow recession behaviours during the drainage of hillslopes where both shallow and deep compartments are involved (Figure 1a.). We specifically provide an exhaustive analysis of the respective controls of 105 both geometrical and hydraulic properties in controlling storage-discharge non-linearity in our two-compartment representation of the subsurface.

Model geometry and structure
We consider a 2-dimensional representation of the subsurface ( Figure 2) which consists of a lower compartment with a horizontal impermeable base and a length = 100 m from the stream location to the groundwater divide. The elevation = 0 110 is set at the river elevation located on the left boundary. The bottom of the model is at -at = 0. We define an upper compartment of constant thickness , leading to a total thickness ( + ) at = 0 that is equal to , and a slope angle = ⁄ , where is the elevation of the interface at = .
In the lower compartment, the saturated hydraulic conductivity and the porosity decrease exponentially with depth 115 following the classically used relationship (e.g. in Cardenas and Jiang (2010)): where = is the elevation of the interface between the two compartments and is the characteristic depth over which hydraulic properties decrease. Initial values of ( =0, =0) and ( =0, =0) were taken arbitrarily at 5 10 −6 / and 0.3 respectively. Note that we found that changing the parameter does not significantly impact our results (see Appendix A. 2). 120 is a coefficient related to the medium structure that we choose equal to 2 as commonly reported in the literature (Bernabé et al., 2003;Cardenas and Jiang, 2010). The upper layer has a homogeneous isotropic saturated hydraulic conductivity and a homogeneous porosity. We explore contrasts in hydraulic conductivities between both compartments according to the ratio as = ( =0, =0) . The porosity of the upper compartment and the interface remain constant across simulations with =

Groundwater flow simulation
We performed Computational Fluid Dynamics (CFD) simulations in the 2-dimensional (2D) hillslope. The flow for variably 135 saturated porous media with incompressible fluid is modelled by solving the Richards' equations with a finite element approach. We set the vertical uphill boundary to no-flow. The vertical downhill boundary is divided into a no-flow boundary below the river elevation and a fixed-head boundary on the top layer. The fixed-head boundary allows groundwater to discharge into the stream at a rate that is driven by the head difference between the aquifer and the river and a conductance term. For simplicity, we set the head of the river to 0 m, meaning that no river stage variation is allowed in our model. We set the 140 conductance between the aquifer and stream equal to the hydraulic conductivity of the upper most aquifer to avoid any additional resistance at the boundary. We assume that the rock compressibility is negligible -by setting it equal to the water compressibility -implying that the storage coefficient only depends on the effective porosity.
The mesh is composed of tetrahedral elements. We imposed a finer mesh in the upper most part of the model where the fastest 145 changes in water table elevation and saturation occur. The initial conditions are of a fully saturated aquifer which then drains towards an asymptotically flat-water table. We excluded the initial times of the simulations, until the water table progressed https://doi.org/10.5194/hess-2022-7 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License. toward a free water table profile from the river to the uphill no-flow boundary. The case of a horizontal homogeneous shallow aquifer, i.e. without the deeper compartment, was checked against the analytical solution from Boussinesq (Boussinesq, 1904) to ensure accurate convergence, i.e. considering Eq. 1 for discharge per unit of river length , with and defined as: 150 a = 4.804 1 2 ⁄ (2 ) 3 2 ⁄ ; = 1.5 Eq. 4

Effect of the deep compartment in the case of a flat aquifer ( = )
First, we review the numerical results for a horizontal homogeneous shallow aquifer overlying a deeper, less permeable, compartment, i.e. the case where β = / = 0. In order to identify the relative influence of the deep aquifer on the recession behaviour for the different values of , we compare in Figure 3 our results with predictions obtained from the classical 155 analytical solution Eq. 1 with recession constants from Eq. 4 which consider only the drainage of a horizontal shallow aquifer with homogeneous hydraulic properties (Boussinesq, 1904). For a review of the different analytical solutions available, the reader is referred to Rupp and Selker (2006) and Dewandel et al. (2003).
The − / vs recession curves shift upward as increases, i.e. increasing rate of change in discharge for the same 160 discharge value as increases, which corresponds to an increase in the value of the recession constant in Eq. 1. This is expected from Eq. 4, with being directly proportional to √ . We found that the analytical solution provides accurate results for high values of (> 16), while deviations are seen when is smaller than 16. At the late times, the analytical solution significantly underestimates − / with tending to lower values as recession progresses. For higher values of , the source of discharge is mainly the upper compartment which satisfies the assumption of only considering the shallower flat 165 aquifer. However, for lower values of , flow tends to be distributed over the full depth of the aquifer system, leading to higher − / as recession advances than the ones predicted by the analytical solution.
We expect that for low values of , will take values close to 1, similar to the case of a thick aquifer with only marginal transmissivity variations with head. In contrast, for a high ratio , the change in transmissivity due to the fast drainage of the 170 upper, more permeable, compartment satisfies the assumption of the Boussinesq solution for a relatively thin unconfined aquifer over a large range in Q. This comparison allows us to identify the influence of the deeper compartment in shaping the recession hydrograph with respect to .

Effect of an inclined interface ( >0) 180
The variety of recession behaviours obtained by changing both and is displayed in  Overall, sharp differences in recession behaviours are evident across all hillslope configurations. At late times, i.e. low 195 discharge , all recessions show increasing towards, and ultimately above, 1.5 when flow is mostly controlled by the lower compartment. Theoretically, the drainage of a deep homogenous aquifer would lead to = 1, i.e. negligible change in transmissivity during recession, whereas in our simulations, the exponential decrease in hydraulic conductivity with depth results in values > 1 at low Q. This latter result agrees qualitatively with the predictions made by Rupp and Selker (2005) showing how a permeability that decreases smoothly with depth increases b above the homogeneous case, up to a value of b 200 = 2. However, our configuration with a discontinuity of hydraulic conductivity, also shows how may not be constant during late times.
As expected from the solution considering the combined discharge from two linear reservoirs (Figure 1, Gao et al. 2017), the maximum values of is strongly influenced by the contrast in hydraulic properties between the two overlapping 205 compartments. We explore the impact of hillslope structural configurations on recession discharge by comparing our simulated https://doi.org/10.5194/hess-2022-7 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.
with the ones predicted from the two-linear reservoir theory. We derived the equation for the maximum value of from the equation provided by Gao et al. (2017) used to predict the evolution of b as a function of discharge (see Appendix A1). We found max( ) to be only a function of between the two linear reservoirs: Eq. 6 We compare in Figure 6 the progression of max( ) as a function of and for different from our simulations against the 210 prediction from Eq. 6. Our two-compartment model shows a non-monotonic behaviour. For ≤ 4, max( ) remains close to the case of the horizontal interface. However, for > 4, max( ) significantly increases before decreasing back to lower values at high . For example, in the case of = 0.5, max( ) increases toward values close to 5 at = 16 before decreasing back to values close to 4 at = 256. This behaviour results from the different flow regimes involved.

215
Eq. 6 fails in predicting the variety of behaviours arising from the different subsurface configurations investigated here. While it might be useful in the case of lateral heterogeneity (Gao et al., 2017;Harman et al., 2009), the parallel reservoir paradigm is not able to predict the recession dynamics involved in the case of a vertically compartmentalized hillslope.  hillslope, normalized by the slope of the interface , and, in Figure 7d, the change in saturation area in the upper compartment, ν u , with respect to the total saturation area, ν, as a function of normalized discharge rates / 0 for the same hillslope configurations. 230 The water table profile for hillslope with = 1 maintains a shape typical of diffusion-controlled drainage. This is confirmed by the linear decrease of the head gradient relative to discharge for = 1 as expected from the Boussinesq equation (Figure   7c provided in linear scales in Appendix A3) and the progressive desaturation of the upper compartment (Figure 7d). In the case of = 64, the early stage of the recession follows almost constant − / , with tending toward values close to 0 235 (even negative at the beginning of the recession) before increasing during a transition period (Figure 5d). This first regime where is close to 0 or negative suggests that drainage is controlled by a kinematic wave regime (Beven, 1982;Rupp and Selker, 2006) where the hydraulic head gradient mostly follows the slope interface, i.e., ⁄ ≈ 1 ( Figure 7c) and a constant or increasing rate of change in discharge is involved. This behaviour is visible in Figure 7b where the water table for = 64 is mostly parallel to, and above, the interface near the start of the recession (Q/Q 0 ≈ 1). As the recession advances, the water 240 table drops rapidly below the interface (Figure 7b) with low saturation of the upper compartment ν u (Figure 7d

Discussion 255
In which geomorphological context might vertical compartmentalization impact recession discharge? Heterogeneity is ubiquitous in the subsurface, inherited from the interplay between tectonic, erosion and weathering processes. St-Clair et al. (2015) presented a high-resolution geophysical imaging of bedrock in three US Critical Zone Observatories. They revealed a complex internal structure of hillslopes with clear differentiation between shallow permeable compartments overlying a deep and less permeable bedrock (Clair et al., 2015). There are multiple factors that may explain the presence of such 260 compartmentalization. Several authors have revisited the mechanical processes that may be responsible for such https://doi.org/10.5194/hess-2022-7 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.
compartmentalization (Clair et al., 2015;Martel, 2006;Moon et al., 2017) and argued that near-surface tectonic and topographic stresses could be the main drivers for the development of permeable fractures organized sub-parallel to the topographic surface, also called exfoliation joints. Similarly, in volcanic regions, basalt formations have been described as highly compartmentalized due to the presence of high densities of fractures formed during cooling of the lava flow. Manga 265 (1996), followed by several other works (Jefferson et al., 2006;Tague et al., 2008), have discussed the role of erosion and weathering in setting shorter drainage characteristic timescales in the Oregonian volcanic Cascades range. It has also been shown that, when hillslopes are under critical stress conditions, the permeability of the shallower compartment may be enhanced by gravitational movements which favour the development of rock damage and fracture density (Vallet et al., 2015;Wolter et al., 2020). Consequently, increased permeability and storage capacities in the shallow subsurface is a common 270 characteristic in bedrock regions. Characterizing the depth and connectivity of such shallow compartments, set during the evolution of the landscape and overprinted by current stress and hydrogeological conditions, remains a main challenge (Dewandel et al., 2006;Guihéneuf et al., 2014;Litwin et al., 2021;Rempe and Dietrich, 2014;Riebe et al., 2017;Worthington et al., 2016).

275
While geophysics reveals the presence of such complexity with unprecedented details (Callahan et al., 2020;Clair et al., 2015;Flinchum et al., 2018;Leone et al., 2020), a challenge resides in attributing effective hydraulic properties to each compartment. This is classically achieved through borehole testing which provides precious local estimates. As examples, we could cite the studies performed in India (Dewandel et al., 2006;Guihéneuf et al., 2014;Maréchal et al., 2004), and the Armoricain massif (Ayraud et al., 2008;Le Borgne et al., 2006;Roques et al., 2014) that have revealed contrast in hydraulic conductivities 280 distributed over 2-3 orders of magnitudes between shallow fractured/weathered aquifer and deeper "fresh" bedrock. Yet the question remains in the representativity of these local estimates in the description of hillslope to catchment-scale flow and transport processes.
To what extent does vertical heterogeneity impact the discharge behaviour of hillslopes? Despite our simplification of 285 the actual heterogeneity of hillslopes, we reveal a strong diversity of impacts on recession behaviour across the different hillslope configurations investigated here. We found that, when compartmentalization is accounted for in the spatial distribution of hydraulic properties, recession behaviour strongly deviates from what is predicted by groundwater theory that considers the drainage of shallow reservoirs with homogeneous properties. Results specifically demonstrate that the critical parameters in the emergence of anomalous behaviour are i) the presence of a deep bedrock aquifer connected to the stream 290 network, ii) the contrast in hydraulic properties between shallow and deep compartments, and iii) the slope of the interface ( Figure 1). Other parameters such as the thickness of the shallow aquifer ( ) and the decay exponent in hydraulic conductivity of the deeper compartment (parameter in Eqs. 2 and 3) have less influence on the variability of recession exponent obtained in our simulations (Appendix A. 2). Noteworthy, our exploration of the decay exponent in hydraulic conductivity remains https://doi.org/10.5194/hess-2022-7 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.
restricted to a few cases of an exponential function. Other functions may deliver different results which would deserve a 295 dedicated study.
In the case of hillslopes with a horizontal interface between shallow and deep compartments, we reveal that there exists a threshold in hydraulic conductivity ratio above which the assumption of only considering the drainage of a homogeneous shallow aquifer derived from the Boussinesq solution seems to hold. In our simulation scheme, the assumption is satisfied for 300 Knowledge regarding the structural configuration and heterogeneity of the aquifer is crucial when interpreting discharge recession behaviours. Our results confirm that the exponent b and its temporal progression are information that can only be interpreted with some amount of local geological knowledge. One of the main challenges in hydrology is to identify 330 the relative contribution among climatic, geomorphological, geological and land use factors in controlling spatial and temporal variabilities of discharge recession behaviour (Blöschl et al., 2019;Gnann et al., 2021). Our results suggest that different recession behaviour is expected depending on the structural configurations, the initial saturation condition and location of the water table prior recession. These are likely to be major factors in explaining the high variabilities in recession constants that are observed in some catchments by recent works (Roques et al., 2017(Roques et al., , 2021Tashie et al., 2020) as it is for transit times and 335 solute transport processes (Ameli et al., 2016;Berghuijs and Kirchner, 2017). We reveal that the temporal progression of head gradient and the relative proportion of saturation volume in the upper compartment with respect to the deeper one are key indicators to predict recession discharge behaviour. These indicators can be assessed through borehole logging and monitoring.
The vertical compartmentalization is classically described through detailed core analysis and geophysical investigation (in hole and surface). The interface can also be approximated through the maximum depth of the boreholes in a region. Indeed, it 340 is common that the drilling of a borehole is stopped when a less productive zone is reached, i.e. when airlift discharges stop increases with depth. While such methods will inform about the geometry and heterogeneity of the hillslope, the pressure monitoring of a network of boreholes distributed along the hillslope will help identify the different flow regimes involved during recession. It can specifically provide estimates in both the progression of head gradients and saturation distribution between upper and deeper compartments that, as we have shown here, are critical parameters to consider when interpreting 345 streamflow recession dynamics.

Conclusion
While it is well accepted from structural and geophysical investigations that hillslopes settled in bedrocks are compartmentalized, its consideration is often overlooked in hydrological models used to predict the evolution of groundwater storage and stream discharge dynamics (Blöschl et al., 2019;Clark et al., 2017;Gnann et al., 2021;Stephens et al., 2021). It 350 is urgent to further develop knowledge on the role of such vertical compartmentalization and its representation at the catchment scale to improve prediction in stream baseflow regimes. Our results allow us to identify that both slope interface between deep bedrock and fractured/weathered shallow compartment and their contrast in hydraulic properties are key parameters to account for when modelling storage-discharge functions of hillslopes. Such characterization, although demanding, is possible through geophysical investigation with combined surface and borehole techniques. While our approach assumes continuity of the 355 compartment along the hillslope, another key parameter influencing stream discharge dynamics during baseflow is related to the volume and spatial connectivity of shallow and deep compartments at the catchment scale. While mechanical and biogeochemical processes may increase the permeability and storage capacity of the bedrock, the spatial representation of this https://doi.org/10.5194/hess-2022-7 Preprint. Discussion started: 31 January 2022 c Author(s) 2022. CC BY 4.0 License.
shallow permeable compartment might be highly heterogeneous, with poor connectivity at the catchment scale (Guihéneuf et al., 2014) and, consequently, displaying higher sensitivity to droughts. 360

Code and data availability
Model results are available at Roques, C. (2022). Model results presented in HESS-2022-7, HydroShare, http://www.hydroshare.org/resource/882cbdcddfdf49ec9ba90c9e43c19cb0. Matlab codes used to process model results are available on request from the authors.

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

Author contribution
CR, DER, GG, ERJ and JSS defined the initial research question. All co-authors were involved in conceptualization. CR