Recession discharge from compartmentalized bedrock hillslopes

We used numerical modelling to explore the role of the 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 holds and, conversely, for those for which it does not. By comparing the modelled recession discharge Q and the groundwater table dynamics, we identified the critical hydrogeological conditions controlling 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 a shallow permeable compartment and deep bedrock one with lower hydraulic properties. Our results confirm that a correct physical interpretation of the recession discharge exponent b from the classical equation $-\mathrm{d}Q/\mathrm{d}t=a{Q}^b$, and its temporal progression, requires knowledge of the structural configuration and heterogeneity of the aquifer.

1 Introduction

Recent streamflow recession analyses have provided new insight into the spatial and temporal variability in hydrological behaviour across the globe (Fan, 2015; Tashie et al., 2020; Schmidt et al., 2020; Roques et al., 2021). These analyses have confirmed that the relationship between groundwater storage and stream discharge (hereafter called the storage–discharge function) 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 nonlinearities in storage–discharge functions may be variable across different recharge events, suggesting dependencies on antecedent catchment conditions (Jachens et al., 2020; Tashie et al., 2020). In the current context where hydrological systems are increasingly impacted by more frequent and intense periods of drought, pushing aquifers and streams to new extremes (Gudmundsson et al., 2021; Famiglietti, 2014), it is urgent that we increase our understanding of the main controlling factors responsible for low-flow discharge dynamics (Fan et al., 2019; Blöschl et al., 2019).

Groundwater theory assuming homogeneous subsurface properties shows that the rate of change in discharge $-\mathrm{d}Q/\mathrm{d}t$ of the falling limb of the hydrograph may scale as a power law of the actual discharge Q, as follows:

$$\begin{array}{}\text{(1)}& -{\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}=a{Q}^b.\end{array}$$

The recession constants a and b arise from the geometric and hydraulic properties of the aquifer system. They are also highly sensitive to the predominant diffusive or kinematic flow conditions (Rupp and Selker, 2006). b is specifically defined as the constant of nonlinearity in the storage–discharge function (Kirchner, 2009). Large values of b are characteristic of resilient streams whose low flows decrease slowly with time (Jachens et al., 2020; Berghuijs et al., 2016; Kirchner, 2009). Models derived from solutions to the 1D Boussinesq equation, assuming homogeneous and horizontal aquifers, predict values from 1 to 1.5 under low-flow conditions (Brutsaert and Nieber, 1977). b ranges from b=1 for confined and thick unconfined aquifers (i.e. much thicker than the depth of the penetrating stream channel), whose saturated thickness only marginally changes, to b=1.5 for a relatively thin unconfined aquifer (Brutsaert and Nieber, 1977; Troch et al., 2013; Wittenberg, 1999). When b=1, the aquifer drainage is exponential, with a characteristic timescale $\mathit{\tau }=\mathrm{1}/a$. When b=1.5, the decrease in the transmissivity with head slows down the aquifer drainage. There is currently a dearth in the definition of contexts, where such simplified effective representations are likely to succeed or fail in accurately describing groundwater storage–discharge dynamics.

Indeed, observed values of b for low-flow conditions may exceed the theoretical values previously described and question the use of homogeneous groundwater flow theories (Tashie et al., 2020; Roques et al., 2017; Jachens et al., 2020; Roques et al., 2021). Previous authors have explored the processes that might be responsible. They include the effect of the spatial and vertical distribution of flow paths (Harman et al., 2009; Rupp and Selker, 2005; Clark et al., 2009), the influence of planform shape of elementary area on flow paths (Paniconi et al., 2003), the spatial variability in recharge and evapotranspiration (Tashie et al., 2020), and the storage from the unsaturated zone (Luo et al., 2018). Quantifying the relative contribution of all those factors in controlling the storage–discharge nonlinearity remains a major challenge in hydrology.

The storage–discharge function becomes nonlinear in catchments where geological heterogeneities affect the partitioning of groundwater flow. Heterogeneity may arise from the presence of different lithologies in the catchment with contrasted permeability. Another source of heterogeneity, which has been described in the referenced studies, is the general decrease of permeability with depth (Anderson, 2015; Boutt et al., 2010; St. Clair et al., 2015; Hencher et al., 2011; Welch and Allen, 2014; Guihéneuf et al., 2014). All conceptual models agree that (1) hydraulic conductivity decreases with depth in response to increasing confining stress and illuviation. Permeability data inferred from borehole tests show an exponential trend with depth that has been quantified in several studies around the world (Saar and Manga, 2004; Ingebritsen and Manning, 1999; Achtziger-Zupančič et al., 2017), and (2) hydraulic properties may be enhanced in the shallower part (up to a depth of 100 m), due to supergene weathering processes and the presence of elevated joint densities in response to unloading and topographic stresses (Harman and Cosans, 2019; Riebe et al., 2017; Rempe and Dietrich, 2014; Moon et al., 2017; Dewandel et al., 2006). 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 the storage–discharge functions of hillslopes.

Figure 1 (a) Conceptual model of a compartmentalized hillslope with groundwater flow distributed between fast flow within the shallow regolith/fractured aquifer and slow flow in the underlain bedrock. (b) Theoretical streamflow behaviour considering two linear reservoirs (b=1) in parallel with fast (blue line; τ=10 d) and slow (red line; τ=100 d) drainage timescales. (c) Resulting recession plot $-\mathrm{d}Q/\mathrm{d}t=f\left(Q\right)$. (d) The evolution of b as a function of discharge identifying high values during the transitional flow regimes up to 3 in this case.

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 in transmissivity with the hydraulic head. The decrease in hydraulic conductivity with depth and vertical compartmentalization plays a significant role in the emergence of anomalous b values. From this perspective, Rupp and Selker (2005, 2006) derived analytical solutions for the case of horizontal (Rupp and Selker, 2005) and sloping (Rupp and Selker, 2006) aquifers, with a power law decrease of hydraulic conductivity with depth such that K≅zⁿ, where z is the elevation, and n is the power law exponent (n>0). The authors showed that b is determined by n, according to $b=\frac{\mathrm{2}n+\mathrm{3}}{n+\mathrm{2}}$ for a horizontal aquifer, such that b is an increasing function of n with a lower limit of 1.5 for n 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 uppermost material, implicitly supposing that there is a shallow, impermeable boundary parallel to the topography. Although this may be a reasonable assumption in some cases, we hypothesize in this paper that, in other contexts, deeper aquifers may also significantly contribute to stream discharge (Fig. 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 contrasting hydraulic conductivities. This structure typically represents a permeable near-surface regolith aquifer underlain by a deeper, less permeable bedrock compartment (Fig. 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 near-surface layer with low characteristic timescales (high value of a as shown in blue in Fig. 1b and c) to the drainage of the deeper less permeable compartment with consequently longer drainage timescale (lower value of a as shown in red in Fig. 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 linear reservoirs (Fig. 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). We will use the two parallel linear reservoirs theory as a reference to describe our results and explore in which context such a simplification might hold to describe storage–discharge functions of vertically compartmented hillslopes.

2 Methodology

We investigate the time dependency of streamflow recession behaviours during the drainage of hillslopes where both shallow and deep compartments are involved (Fig. 1a). We specifically provide a thorough analysis of the respective controls of both geometrical and hydraulic properties in controlling storage–discharge nonlinearity in our two-compartment representation of the subsurface.

2.1 Model geometry and structure

We consider a 2-dimensional (2D) representation of the subsurface (Fig. 2a), which consists of a lower compartment with a horizontal impermeable base and a length of L=100 m from the stream location to the groundwater divide. The elevation z=0 is set at the river elevation located on the left boundary. The bottom of the model is at −ζ at x=0. We define an upper compartment of constant thickness D, leading to a total thickness (D+ζ) at x=0 that is equal to L and a slope angle $\mathit{\beta }=e/L$, where e is the elevation of the interface at x=L.

Figure 2 (a) A 2D model cross section with the main geometrical parameters identified along with model boundaries and groundwater flow lines. Panels (b) and (c) are vertical profiles taken at x=0 of K and θ, respectively, for the different α values investigated.

In the lower compartment, the saturated hydraulic conductivity K_L and the porosity θ_L decrease exponentially with depth (Fig. 2b and c), following the classically used relationship (e.g. in Cardenas and Jiang, 2010):

$$\begin{array}{}\text{(2)}& K_{\mathrm{L}}(x,z)=K_{\mathrm{L}(x=\mathrm{0},z=\mathrm{0})}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}\left(z_{\mathrm{s}}\left(x\right)-z\right)\right),\end{array}$$

and

$$\begin{array}{}\text{(3)}& {\mathit{\theta }}_{\mathrm{L}}(x,z)={\mathit{\theta }}_{\mathrm{L}(x=\mathrm{0},z=\mathrm{0})}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }\mathit{\eta }}}\left(z_{\mathrm{s}}\left(x\right)-z\right)\right),\end{array}$$

where z_s=βx is the elevation of the interface between the two compartments, and α is the characteristic depth over which hydraulic properties decrease. Initial values of $K_{\mathrm{L}(x=\mathrm{0},z=\mathrm{0})}$ and ${\mathit{\theta }}_{\mathrm{L}(x=\mathrm{0},z=\mathrm{0})}$ were taken arbitrarily at $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ and 0.3, respectively. η is a coefficient related to the medium structure that we chose to be equal to 2, as commonly reported in the literature (Cardenas and Jiang, 2010; Bernabé et al., 2003). The upper layer has a homogeneous isotropic saturated hydraulic conductivity, K_U, and a homogeneous porosity, θ_U. We explore contrasts in hydraulic conductivities between both compartments according to the ratio r_K, as $K_{\mathrm{U}}=r_K{K}_{\mathrm{L}(x=\mathrm{0},z=\mathrm{0})}$. The porosity of the upper compartment and the interface remain constant across simulations with ${\mathit{\theta }}_{\mathrm{U}}={\mathit{\theta }}_{\mathrm{L}(x=\mathrm{0},z=\mathrm{0})}$.

The geometrical parameters β and D are varied in order to explore different hillslope configurations, with $D/L=[\mathrm{0.1},\mathrm{0.2},\mathrm{0.5}]$ and $\mathit{\beta }=e/L$ = [0.0, 0.1, 0.2, 0.5]. The exponential decrease in hydraulic properties of the lower compartment is changed by varying α = [20, 100, 200] m in Eqs. (2) and (3). Finally, the contrast in hydraulic conductivities is explored over 2 orders of magnitude with r_K = [1, 4, 16, 64, 256].

2.2 Groundwater flow simulation

We performed computational fluid dynamics (CFDs) simulations in the 2D hillslope. The flow for variably saturated porous media with incompressible fluid is modelled by solving the Richards' equations with a finite element approach implemented in COMSOL Multiphysics^®. 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 the 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 conductance between the aquifer and stream equal to the hydraulic conductivity of the uppermost aquifer K_U 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 uppermost part of the model where the fastest 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 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 (1904) to ensure accurate convergence, i.e. considering Eq. (1) for discharge per unit of river length Q, with a and b defined as follows:

$$\begin{array}{}\text{(4)}& a={\displaystyle \frac{\mathrm{4.804}{K}^{\mathrm{1}/\mathrm{2}}}{\mathit{\theta }(\mathrm{2}L{)}^{\mathrm{3}/\mathrm{2}}}};b=\mathrm{1.5}.\end{array}$$

3 Results

3.1 Effect of the deep compartment in the case of a flat aquifer (β=0)

First, we review the numerical results for a horizontal homogeneous shallow aquifer overlying a deeper, less permeable, compartment, i.e. the case where $\mathit{\beta }=e/L=\mathrm{0}$. In order to identify the relative influence of the deep aquifer on the recession behaviour for the different values of r_K, we compare in Fig. 3 our results with the predictions obtained from the classical 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).

Figure 3 Recession plots displaying $-\mathrm{d}Q/\mathrm{d}t$ vs. Q in a log–log space obtained for the horizontal interface, i.e. β=0, and different hydraulic conductivity ratios r_K. D is fixed at 20 m and α at 100 m. The dashed lines show predictions for a flat and homogeneous aquifer, without the deep compartment combining Eqs. (1) and (4) (b=15; Boussinesq, 1904). The slope for b=1, characteristic of the drainage of a linear reservoir, has been drawn for reference.

The $-\mathrm{d}Q/\mathrm{d}t$ vs. Q recession curves shift upward as r_K increases, i.e. increasing the rate of change in the discharge for the same discharge value as r_K increases, which corresponds to an increase in the value of the recession constant a in Eq. (1). This is expected from Eq. (4), with a being directly proportional to $\sqrt{K}$. We found that the analytical solution provides accurate results for high values of r_K(>16) for long portions of the recession, while deviations are seen when r_K is smaller than 16. At late times, the analytical solution significantly underestimates $-\mathrm{d}Q/\mathrm{d}t$, with b tending to lower values as recession progresses. For higher values of r_K, the source of discharge is mainly the upper compartment, which satisfies the assumption of only considering the shallower flat aquifer. However, for lower values of r_K, flow tends to be distributed over the full depth of the aquifer system, leading to higher $-\mathrm{d}Q/\mathrm{d}t$, as the recession advances, than the ones predicted by the analytical solution.

We expect that, for low values of r_K, b 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 r_K, the change in transmissivity due to the fast drainage of the 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 r_K.

3.2 Effect of an inclined interface (β>0)

The variety of recession behaviours obtained by changing both r_K and β is displayed in Fig. 4. Considering cases where β>0, and by increasing r_K, values of $-\mathrm{d}Q/\mathrm{d}t$ increase near the onset of the recessions. With time, the water table progressively lowers below the interface, inducing a sharp reduction in the effective hydraulic conductivity that translates to larger values of b. Finally, the recession progressively converges to a value of a that is characteristic of the geometrical and hydraulic properties of the deeper compartment, i.e. the blue lines at late times for r_K=1. The changes in local b values as a function of Q for each hillslope configuration b_i are shown in Fig. 5. It is computed as follows:

$$\begin{array}{}\text{(5)}& b_i=\log \left({\left(-{\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}\right)}_{i+\mathrm{1}}-{\left(-{\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}\right)}_i\right)/\log \left(Q_{i+\mathrm{1}}-Q_i\right)\end{array}$$

where i and i+1 are indexes of the discharge range considered [Q_min:Q_max].

Figure 4 Recession plots displaying $-\mathrm{d}Q/\mathrm{d}t$ vs. Q in a log–log space obtained for different hydraulic conductivity ratios r_K and different slope angles of the interface between shallow and deep compartments β. D is fixed at 20 m and α at 100 m. Lines of b=3 and b=15 are drawn for reference.

Figure 5 Evolution of local b values as a function of discharge Q for the range of hydraulic conductivity ratios r_K and different slope angles of the interface between shallow and deep compartments β. D is fixed at 20 m and α at 100 m.

Overall, sharp differences in recession behaviours are evident across all hillslope configurations. At late times, i.e. low discharge Q, all recessions show increasing b towards, and ultimately above, 1.5 when the flow is mostly controlled by the lower compartment. Theoretically, the drainage of a deep homogenous aquifer would lead to b=1, i.e. a negligible change in transmissivity during recession, whereas, in our simulations, the exponential decrease in hydraulic conductivity with depth results in b 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=2. However, our configuration with a discontinuity of hydraulic conductivity also shows how b may not be constant during late times.

As expected from the solution considering the combined discharge from two parallel linear reservoirs (Fig. 1; Gao et al., 2017), the maximum values of b is strongly influenced by the contrast in hydraulic properties r_K between the two overlapping compartments. Here we investigate the impact of hillslope structural configurations on recession discharge by comparing our simulated b with the ones predicted from the two parallel linear reservoirs. We derived the equation for the maximum value of b from the equation provided by Gao et al. (2017) used to predict the evolution of b as a function of discharge (see Appendix A). We found max(b) to only be a function of r_K between the two parallel linear reservoirs with the following:

$$\begin{array}{}\text{(6)}& \max \left(b\right)={\displaystyle \frac{{\left(r_K+\mathrm{1}\right)}^{\mathrm{2}}}{\mathrm{4}{r}_K}}.\end{array}$$

We compare, in Fig. 6, the progression of max(b) as a function of r_K and, for different β, our simulations against the prediction from Eq. (6). Our two-compartment model shows a nonmonotonic behaviour. For r_K≤4, max(b) remains close to the case of the horizontal interface, i.e. max(b) ranging from 1 to 2. However, increasing r_K above 4 leads to a progressive increase in max(b) before decreasing back to lower values at high r_K. For example, in the case of β=0.5, max(b) takes value of ∼5 at r_K=16 before decreasing to approximately 4 at r_K=256. This behaviour results from the different flow regimes involved.

Figure 6 Maximum values of b as a function of r_K. The different slope angles of the interface β are pictured by different symbols. The grey dashed line represents the Gao et al. (2017) model, predicting the evolution of maximum values of b as a function of r_K, considering the drainage of two linear reservoirs in parallel.

We found that Eq. (6) (dashed line in Fig. 6) fails to predict the variety of behaviours arising from the different subsurface configurations investigated here. While it might be useful in the case of lateral heterogeneity (Harman et al., 2009; Gao et al., 2017), the parallel reservoir paradigm is not able to predict the recession dynamics involved in the case of a vertically compartmentalized hillslope.

3.3 Effect of changing thickness of the upper compartment, D, and characteristic length of permeability and porosity decrease with depth α

We display, in Fig. 7, the variability in the recession behaviours obtained by jointly changing r_K, the upper compartment thickness D (Fig. 7a–c), and the characteristic length of permeability and porosity decrease with depth α (Fig. 7d–f). In this exploration, we fixed the value of β at 0.2. The results revealed that r_K controls most of the variability in the recession behaviours obtained across the different hillslope configurations, while changing D and α only have little impact. This is confirmed by analysing the progression of max(b) as a function of r_K for different D and α (Appendix B; Fig. B1), which generally shows low sensitivity to these parameters. However, one may note that the impact of both D and α appear not to be negligible when r_K≥64. In this context, increasing the values of both parameters tends to lower the average and maximum values of b. This decrease in b is a consequence of the homogenization of the hydraulic properties when increasing D and α across the hillslope, which tends to simplify groundwater flow and diminish the nonlinearity in its storage–discharge functions.

Figure 7 Evolution of local b values as a function of discharge Q for the range of hydraulic conductivity ratios r_K and different upper compartment thickness D (a–c) and characteristic length of permeability and porosity decrease with depth α (d–f). β is fixed at a value of 0.2.

3.4 Impact of transient flow regimes on b values with steep interfaces (β≫0): transition from gravity- to diffusion-controlled regimes

To illustrate the link between the transient pattern of the water table – and, consequently, the flow field – and its impact on the evolution of b, we show, in Fig. 8, the progression of the water table for two hillslopes under β=0.5, with r_K=1 (Fig. 8a) and r_K=64 (Fig. 8b). We also plot, in Fig. 8c, the progression of the average head gradient ∇H computed across the hillslope, normalized by the slope of the interface β, and, in Fig. 8d, the change in the saturation area in the upper compartment, ν_u, with respect to the total saturation area, ν, as a function of normalized discharge rates $Q/Q_{\mathrm{0}}$ for the same hillslope configurations.

Figure 8 The uppermost plots in panels (a) and (b) represent the groundwater table profiles taken at different $Q/Q_{\mathrm{0}}$ during the recession, considering (a) β=0.5 and r_K=1 and (b) β=0.5 and r_K=64. The actual values of $Q/Q_{\mathrm{0}}$ correspond to the ones used in the lowermost plots in panels (c) and (d). The upper–lower compartment interface is represented by the dashed black line in the background. Water tables are shown as a continuous line, with colours scaling with the corresponding b values. Panel (c) indicates the progression of the average head gradient ∇H, normalized by the slope of the interface β. In panel (d), the change in the saturation area of the upper compartment ν_u, with respect to the total saturation area ν, as a function of $Q/Q_{\mathrm{0}}$ for the same hillslope configurations as in panels (a) and (b), is shown. The triangular symbols are for the case where r_K=1, while diamond ones are for r_K=64. The colours of the symbols show the magnitude of b during the recession.

The water table profile for a hillslope with r_K=1 maintains a shape typical of diffusion-controlled drainage. This is confirmed by the linear decrease in the head gradient relative to discharge for r_K=1, as expected from the Boussinesq equation (Fig. 8c is provided in linear scales in Fig. C1) and the progressive desaturation of the upper compartment (Fig. 8d). In the case of r_K=64, the early stage of the recession follows almost constant $-\mathrm{d}Q/\mathrm{d}t,$ with b tending toward values close to 0 (and even negative at the beginning of the recession) before increasing during a transition period (Fig. 5d). This first regime, where b 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. $\mathrm{\nabla }H/\mathit{\beta }\approx \mathrm{1}$ (Fig. 8c), and a constant or increasing rate of change in discharge is involved. This behaviour is visible in Fig. 8b, where the water table for r_K=64 is mostly parallel to, and above, the interface near the start of the recession ($Q/Q_{\mathrm{0}}\approx \mathrm{1}$). As the recession advances, the water table drops rapidly below the interface (Fig. 8b), with a low saturation of the upper compartment ν_u (Fig. 8d), whereby $\mathrm{\nabla }H/\mathit{\beta }$ starts to decrease and discharge rates transition toward a diffusion-controlled regime. During the transition period, the hydraulic gradient remains steep, and the peak in b values occurs. Finally, at the late times, the water table progressively disconnects from the interface, and b values decrease back toward 1.5, as expected from the diffusion-controlled regimes under homogeneous conditions.

4 Discussion

4.1 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. For example, St. Clair et al. (2015) presented a high-resolution geophysical imaging of bedrock in three U.S. Critical Zone Observatories. They revealed a complex internal structure of hillslopes, with a clear differentiation between shallow permeable compartments overlying a deep and less permeable bedrock. There are multiple factors that may explain the presence of such a compartmentalization. Several authors have revisited the mechanical processes that may be responsible for such a compartmentalization (Moon et al., 2017; Martel, 2006; St. Clair et al., 2015) and argued that near-surface tectonic and topographic stresses could be the main drivers for the development of permeable fractures organized subparallel to the topographic surface, the so-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 the cooling of the lava flow. Manga (1996), followed by several other works (Tague et al., 2008; Jefferson et al., 2006), has 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 (Wolter et al., 2020; Vallet et al., 2015). Consequently, increased permeability and storage capacities in the shallow subsurface are common characteristics 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 (Riebe et al., 2017; Rempe and Dietrich, 2014; Dewandel et al., 2006; Guihéneuf et al., 2014; Worthington et al., 2016; Litwin et al., 2022).

While geophysics reveals the presence of such complexity with unprecedented details (Flinchum et al., 2018; St. Clair et al., 2015; Leone et al., 2020; Callahan 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; Maréchal et al., 2004; Guihéneuf et al., 2014) and the Armorican Massif (Roques et al., 2014; Ayraud et al., 2008; Le Borgne et al., 2006) that have revealed contrasting hydraulic conductivities distributed over 2–3 orders of magnitude between shallow fractured/weathered aquifer and deeper, fresh bedrock. Yet the question remains in the representativity of these local estimates in the description of a hillslope to catchment-scale flow and transport processes.

4.2 To what extent does vertical heterogeneity impact the discharge behaviour of hillslopes?

Despite our simplification of 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 vertical 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 network, (ii) the contrast in hydraulic properties between shallow and deep compartments, and (iii) the slope of the interface (Fig. 1). Other parameters, such as the thickness of the shallow aquifer (D) and the decay exponent in hydraulic conductivity of the deeper compartment (parameter α in Eqs. 2 and 3), have less of an influence on the variability in recession exponent b obtained in our simulations (Fig. 7 and Appendix B). It is noteworthy that our exploration of the decay exponent in hydraulic conductivity and porosity remains restricted to a few cases of an exponential function. Other functions may deliver different results which would deserve a 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 r_K above which the assumption of only considering the drainage of a homogeneous shallow aquifer derived from the Boussinesq solution seems to hold. The assumption is satisfied when the lateral flow in the shallower compartment prevails, i.e. with only limited contribution from the deep bedrock to the discharge dynamics (when r_K≥16 in our hillslope configuration). However, for lower values of r_K, the continuity in hydraulic properties between shallow and deep compartments favours deeper flow partitioning. In this case, the recession behaviour deviates from the prediction of the Boussinesq equation. The values of b tend toward values close to 1, typical of the drainage of thick aquifers where the change in transmissivity involved by the desaturation during recession is negligible. At later times, b values were found to increase, which are again controlled by the decay in hydraulic conductivity, as shown in Fig. 5a, for $Q<{\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹ and r_K=1 and 4. This shows the intrinsic variability in recession behaviour, even in a relatively simple configuration, which challenges the use of homogeneous groundwater theories.

Considering a steep slope interface, β=0.5, and a strong contrast in hydraulic conductivities, r_K=64, sharp changes in flow regimes are involved when the water table transitions from the shallower compartment toward the deeper one with three separated phases (Fig. 8c and d). In the first phase, during high-flow regimes, when $Q/Q_{\mathrm{0}}>\mathrm{0.8}$ (i.e. the two first points on the left of Fig. 8c and d), the saturated volume drops by a factor of 4 with a minimal change in discharge rates (Fig. 7d). The upper zone drains quickly because of the combined effect of its higher conductivity and of the inclined interface. The slope has a strong control on the head gradient, causing high discharge rates and even achieving a kinematic flow regime. Flows occur predominantly in the upper compartment. In the second phase, discharge rates decrease sharply while transitioning to the lower compartment at an almost constant head gradient. The sharp increase in the value of b can be explained in the following way. As $b=\mathrm{d}\left(\log \right(-\mathrm{d}Q/\mathrm{d}t\left)\right)/\mathrm{d}\left(\log Q\right)$ and Q∼HK∇H, for b to be large, K must be decreasing rapidly with the otherwise steady H∇H values. In the third phase, the head gradient decreases again with a flow mostly occurring in the lower compartment. The value of b decreases quickly with Q. The water table becomes parabolic uphill and remains influenced by the proximity of the interface close to the river (Fig. 8b), thus explaining the still significantly larger values of b than in the horizontal interface case (Fig. 8a). The presence and geometry of the interface between high and low permeable compartments and the decreasing trend of permeability with depth are keys in conditioning the emergence of complex groundwater flow processes and recession discharge behaviours. It is specifically true in this hillslope configuration for which the river channel fully penetrates the upper compartment. This can be considered as being a common case for landscapes under sufficient erosion rates to carve the weaker compartment, i.e. concerned by higher fracture density or higher porosity, until a more resistant fresh bedrock is reached. The case where high erosion rates may lead to the formation of gorges and canyons in the fresh bedrock may result in different behaviours. The impact of such specific configurations on recession behaviour remains to be explored.

With K_L chosen at $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹, simulated recession behaviours are typical of fractured igneous and metamorphic bedrock (Dewandel et al., 2006). Initial values of K_L, and the porosity θ_L, may impact results and interpretation if different flow regimes are involved when increasing r_K. Exploring different values of K_L and θ_L typical of other lithologies is an important perspective of this work. Similarly, our experiment considers a constant value of porosity for the upper compartment across simulations. While this assumption neglects the coevolution of θ_U and K_U during geomechanical and geochemically driven processes shaping the upper compartment of the aquifer, it allowed us to focus our interpretation on the impact of hydraulic conductivity enhancement. Investigating the impact of porosity contrast remains also to be explored.

4.3 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 provide information that can only be interpreted with some amount of local geological knowledge. One of the main challenges in hydrology is to identify the relative contribution among climatic, geomorphological, geological, and land use factors in controlling spatial and temporal variabilities in discharge recession behaviour (Blöschl et al., 2019; Gnann et al., 2021; Hopp and McDonnell, 2009). 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 to recession. These are likely to be major factors in explaining the high variabilities in recession constants that have been observed in some catchments by recent works (Tashie et al., 2020; Roques et al., 2017, 2021) and for transit times and solute transport processes (Berghuijs and Kirchner, 2017; Ameli et al., 2016, 2021). We reveal that the temporal progression of the head gradient and the relative proportion of saturation volume in the upper compartment with respect to the deeper one are key indicators for predicting 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 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 us about the geometry and heterogeneity of the hillslope, the pressure monitoring of a network of boreholes distributed along the hillslope will help us identify the different flow regimes involved during recession. It can specifically provide estimates of 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 streamflow recession dynamics.

5 Conclusion

While it is well accepted from structural and geophysical investigations that hillslopes settled in bedrock are compartmentalized, this 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; Stephens et al., 2021; Gnann et al., 2021). It is urgent to further develop knowledge on the role of such a vertical compartmentalization and its representation at the catchment scale to improve the prediction in stream baseflow regimes. Our results allow us to identify that both the slope interface between deep bedrock and a fractured/weathered shallow compartment and their contrast in hydraulic properties are key parameters to account for when modelling storage–discharge functions of hillslopes. Such a characterization, although demanding, is possible through geophysical investigation with combined surface and borehole techniques. While our approach assumes the continuity of the 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 shallow permeable compartment might be highly heterogeneous, with poor connectivity at the catchment scale (Guihéneuf et al., 2014) and, consequently, display higher sensitivity to droughts.

Appendix A: The two parallel linear reservoirs theory

We can define the total outflow from the drainage of two linear reservoirs as follows:

$$\begin{array}{}\text{(A1)}& Q\left(t\right)=Q_{\mathrm{f}}\left(t\right)+Q_{\mathrm{s}}\left(t\right),\end{array}$$

with subscripts f representing fast drainage and s slow drainage, respectively.

Linear reservoir theory defines the evolution of the discharge as a function of time, as follows:

$$\begin{array}{}\text{(A2)}& Q_{\mathrm{f}}\left(t\right)=Q_{\mathrm{f}}\left(\mathrm{0}\right)\exp \left(-{\displaystyle \frac{t}{k_{\mathrm{f}}}}\right),\end{array}$$

and

$$\begin{array}{}\text{(A3)}& Q_{\mathrm{s}}\left(t\right)=Q_{\mathrm{s}}\left(\mathrm{0}\right)\exp \left(-{\displaystyle \frac{t}{k_{\mathrm{s}}}}\right),\end{array}$$

with k_f and k_s as the recession constants of the fast and slow reservoirs, respectively, and Q_f(0) and Q_s(0) as the initial flow rates prior recession from the fast and slow reservoirs, respectively.

k is related to the geometrical and hydraulic properties of the reservoir by the following (Boussinesq, 1904):

$$\begin{array}{}\text{(A4)}& k={\displaystyle \frac{\mathrm{4}\mathit{\varphi }{L}^{\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{2}}}}KH.\end{array}$$

Considering r_f the ratio of fast flow to the total flow, we can write the following:

$$\begin{array}{}\text{(A5)}& Q\left(t\right)=Q\left(\mathrm{0}\right)\left[\left(\mathrm{1}-r_{\mathrm{f}}\left(\mathrm{0}\right)\right)\exp \left(-{\displaystyle \frac{t}{k_{\mathrm{s}}}}\right)+r_{\mathrm{f}}\left(\mathrm{0}\right)\exp \left(-{\displaystyle \frac{t}{k_{\mathrm{f}}}}\right)\right].\end{array}$$

Rearranging to solve for r_f(t) yields the following:

$$\begin{array}{}\text{(A6)}& r_{\mathrm{f}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\left(\frac{\mathrm{1}}{r_{\mathrm{f}}\left(\mathrm{0}\right)}-\mathrm{1}\right)\exp \left(\left(\frac{\mathrm{1}}{k_{\mathrm{f}}}-\frac{\mathrm{1}}{k_{\mathrm{s}}}\right)t\right)}}.\end{array}$$

In Gao et al. (2017), the authors derived the solution for solving b as a function of r_f(t), as follows:

$$\begin{array}{}\text{(A7)}& b={\displaystyle \frac{\mathrm{1}+\left({\left(\frac{k_{\mathrm{s}}}{k_{\mathrm{f}}}\right)}^{\mathrm{2}}-\mathrm{1}\right)r_{\mathrm{f}}\left(t\right)}{{\left(\mathrm{1}+\left(\frac{k_{\mathrm{s}}}{k_{\mathrm{f}}}-\mathrm{1}\right)r_{\mathrm{f}}\left(t\right)\right)}^{\mathrm{2}}}}.\end{array}$$

If we replace r_f(t) with Eq. (A6) and derive the solution for the maximum value of b, we find the following:

$$\begin{array}{}\text{(A8)}& \max \left(b\right)={\displaystyle \frac{{\left(\frac{k_{\mathrm{s}}}{k_{\mathrm{f}}}+\mathrm{1}\right)}^{\mathrm{2}}}{\mathrm{4}\frac{k_{\mathrm{s}}}{k_{\mathrm{f}}}}}.\end{array}$$

Considering that k_f and k_s are only differentiated by their hydraulic conductivities, we can write Eq. (A8) as a function of r_K.

$$\begin{array}{}\text{(A9)}& \max \left(b\right)={\displaystyle \frac{{\left(r_K+\mathrm{1}\right)}^{\mathrm{2}}}{\mathrm{4}{r}_K}}.\end{array}$$

Appendix B: Sensitivity to parameter D and α

Figure 8c and d are displayed with a linear scale for $Q/Q_{\mathrm{0}}$.

Figure B1 Maximum values of b as a function of r_K for β=0.2 and (a) variations in the thickness of the upper compartment, D, and (b) variations in the characteristic depth over which hydraulic properties decrease, α.

Appendix C: Figure 8c and d displayed with a linear scale for $Q/Q_{\mathrm{0}}$

Figure C1 (a) The progression of the average head gradient ∇H, normalized by the slope of the interface β. (b) The change in saturation area of the upper compartment ν_u with respect to the total saturation area ν, as a function of $Q/Q_{\mathrm{0}}$. Triangular symbols are for the case where r_K=1, while diamond ones are for r_K=64. The colours of the symbols show the magnitude of b during recession.