# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2024GL113428

**Key Points:**

* This study examines the combined effects of topographic and tectonic stress on hydraulic property distribution at the hillslope scale

* Steep hillslopes under stress show heterogeneity in hydraulic properties, diverging from traditional depth-dependent conceptual models

* Stress-dependent heterogeneity causes groundwater flow paths to channelize near the surface, deviating from homogeneous model predictions

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

R. Figueroa,
ronny.figueroa@unine.ch

**Citation:**

Figueroa, R., Abhervé, R., Halloran, L. J. S., Valley, B., & Roques, C. (2025). Impacts of stress-dependent hydraulic properties on hillslope-scale groundwater flow. *Geophysical Research Letters*, *52*, e2024GL113428. https://doi.org/10.1029/2024GL113428

Received 31 OCT 2024

Accepted 10 JAN 2025

© 2025. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

<figure><img src="image_1.png" /></figure>

## Impacts of Stress-dependent Hydraulic Properties on Hillslope-Scale Groundwater Flow

Ronny Figueroa¹ ID, Ronan Abhervé¹ ID, Landon J. S. Halloran¹ ID, Benoît Valley¹ ID, and Clément Roques¹ ID

¹Center for Hydrogeology and Geothermics (CHYN), University of Neuchâtel, Neuchâtel, Switzerland

**Abstract** Groundwater flow at the hillslope scale is controlled by hydraulic properties and topography. Geological media are inherently heterogeneous, with properties depending on both lithology and external factors such as weathering, erosion, and tectonic activity. The effects of stress-dependent variations in porosity and hydraulic conductivity on groundwater flow and residence times remain unexplored. Using process-based hillslope models, we demonstrate how external lateral stresses cause significant variations in stress and strain, particularly on slopes steeper than 30%. These stresses reorganize porosity and hydraulic conductivity, with compressive regimes reducing these properties at valley bottoms and tensional regimes elevating values at crests. This redistribution of porosity and hydraulic conductivity leads to a channeling of groundwater flow near the surface. Channeling increases discharge rates in the upper seepage zone and the proportion of rapid flows. Our results highlight geomechanical-hydraulic property relationships as critical controls on critical zone processes in mountain settings.

**Plain Language Summary** We investigate how groundwater flows through rock under the influence of tectonic and topographic forces. The weight of mountains and movement of tectonic plates can locally alter the pore spaces through which groundwater flows. We determine that, in specific tectonic and topographic settings, traditional hydrological models may inaccurately predict water availability and quality. Groundwater would flow quicker near the surface than previously thought, increasing its vulnerability to climate change. However, deeper, slower flow paths provide steady discharge to valleys, supporting long-term water supply.

## 1. Introduction

Groundwater flow and transport at hillslope to catchment scales are significantly influenced by geomorphology and geology (Ameli et al., 2016; Bresciani et al., 2014; Cardenas & Jiang, 2010; Fan et al., 2019; Gauvain et al., 2021; Goderniaux et al., 2013; Leray et al., 2016, 2019). Geological environments are inherently heterogeneous in terms of hydraulic properties due to variations in lithology, which are further modified by tectonic processes, weathering, and erosion. Although numerous hydraulic tests have been conducted in diverse geological settings worldwide and compiled in public databases (Achtziger-Zupančič et al., 2017; Ingebritsen & Glee-son, 2015; Ingebritsen & Manning, 2010; Manning & Ingebritsen, 1999; Maréchal et al., 2004; Masset & Loew, 2010; Stober & Bucher, 2015), parameterizing hydrological models with hydraulic conductivity and porosity that properly capture groundwater dynamics remains a challenge. In particular, the variability of hydraulic properties at depth is classically implemented in hydrogeological models via a power law or exponential decay, following the trends determined from hydraulic test data (Achtziger-Zupančič et al., 2017; Manning & Ingebritsen, 1999), overlooking the influence of topographic and tectonic stresses.

Steep topography, combined with tectonic strain, affects the stress distribution in the hillslope. This results in compressional regimes at the valley floor and tensional regimes at the slope crest (L. Liu & Zoback, 1992; Martel, 2000; Martel & Muller, 2000; Miller & Dunne, 1996; Moon et al., 2017; Savage et al., 1985; Slim et al., 2015). Since hydraulic properties are stress-dependent (Cappa & Rutqvist, 2011; H.-H. Liu et al., 2013; Oda et al., 2002, 2023; Oestreicher et al., 2023), it is expected that valleys exhibit lower values of hydraulic conductivity and porosity, while the crests should display higher values. This distribution has been recently imaged in geophysical studies (Chen et al., 2023; Moon et al., 2017; Nielson et al., 2021; Riebe et al., 2017; Slim et al., 2015; St. Clair et al., 2015) which reveal that bedrock near the surface at valley floors would have a lower fracture density compared to the crests. This spatial distribution of aquifer properties has been referred to as a surface mirror type, due to the geometric reflection that occurs with respect to topography (Moon et al., 2017; Riebe
et al., 2017; St. Clair et al., 2015). This conceptualization of the structure of the critical zone has significant implications for groundwater recharge, dynamic storage, streamflow generation, and solute transport (Ameli et al., 2016; Cardenas & Jiang, 2010; Chen et al., 2024; Roques et al., 2022; Rupp & Selker, 2006). However, it remains unclear under what specific topographical and stress conditions this heterogeneity must be accounted for in hydrogeological models to accurately predict flow and transport processes.

In this study, we investigate the combined effects of topographic and tectonic stresses on the distribution of porosity and hydraulic conductivity at the hillslope scale and their impact on groundwater flow partitioning. Using a 2D hillslope geomechanical model, we simulate stress and strain distribution across synthetic hillslopes under different topographic and tectonic stress conditions. The resulting strain fields are then converted into porosity and hydraulic conductivity fields, which serve as inputs for a steady state groundwater flow model. The analysis focuses on quantifying the impacts of hillslope scenarios on flow channeling, spatial variations of seepage rates, and residence time distributions as key hydrological indicators.

## 2. Modeling Approach

### 2.1. Geomechanical Model Setup

Two processes were considered to simulate the distribution of stresses and strains in the hillslope: gravitational loading resulting from own weight of the material and Lateral Tectonic Strain (LTS). Both processes interact with the topography and generate heterogeneous stress and strain distribution in the hillslope.

We simulate these processes in a 2D synthetic hillslope model with horizontal dimensions of $L = 1000$ m between the valley bottom and the ridge, and vertical dimension $2L$ between the valley bottom and the base of the model (Figure 1a, red rectangle). To minimize boundary effects, the model includes three valleys, with the central valley representing the hillslope of interest for the groundwater model (Figure 1a). This configuration was determined through a sensitivity analysis on how the number of valleys affects the stability of results in the central area. The stress and strain data from one of the central slopes are then used to compute the hydrogeological properties. We parameterized the model with a purely elastic, homogeneous and isotropic material with the following mechanical properties: Young's modulus, $E = 20$ GPa, Poisson's ratio, $\nu = 0.3$, and specific weight, $W = \rho g = 27$ kN/m$^3$. Using the finite element geomechanical simulator Rocscience-RS2 (Rocscience Inc, 2001), we computed stress and strain fields by superposing the results from gravity loading with lateral constraints (rollers on model sides and bottom and free boundary on model top) and LTS simulated by a lateral (i.e., purely horizontal) displacement boundary condition imposed on the right side of the model (red arrows in Figure 1a).

To test the influence of various slopes and strain conditions, we defined nine scenarios, combining three slope angles with three levels of LTS (Table S1 in Supporting Information S1). The slope gradients considered are 10%, 30%, and 50%. The tectonic strain $\epsilon_{xx}$, defined as the ratio between the imposed boundary displacement and the model width, ranges from 0% (No-LTS, i.e., only gravitational loading), to 0.2% (Low-LTS), and 0.4% (High-LTS). In this study, the scenarios are referred to using the nomenclature $C_{\text{slope-LTS}}$, where the first subscript indicates the slope and the second indicates the strain level, separated by a hyphen. For example, $C_{10-high}$ refers to a scenario with a 10% slope and High-LTS. Volumetric strain is exported from these simulations as an input for hydraulic properties computations.

### 2.2. Groundwater Flow and Transport Simulations

2D steady-state groundwater flow simulations are performed considering the porosity ($\theta$) and hydraulic conductivity ($K$) fields calculated from the volumetric strain fields simulated by the geomechanical models (Figure 1b) using Equations 1 and 2, proposed by Chin et al. (2000):

$$ \theta = 1 - (1 - \theta_0) e^{\epsilon_v} \quad (1) $$

$$ K = K_0 \left( \frac{\theta}{\theta_0} \right)^n \quad (2) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Numerical experiments and hydromechanical results. (a) Diagram of the geomechanical model setup and boundary conditions. The red area was considered for groundwater flow modeling. (b) Volumetric strain from the geomechanical model, (c) Normalized hydraulic conductivity, K/K₀, calculated from the volumetric strain. K₀ is the hydraulic conductivity at ε_v = 0 and Equation 2. (d) K/K₀ as a function of normalized depth below surface for the three Lateral Tectonic Strain (LTS) cases investigated.</figcaption></figure>

Here, $\theta$ [-] is the porosity at a given volumetric strain, $\varepsilon_v$ [-] is the volumetric strain, $\theta_0$ [-] is the porosity at $\varepsilon_v$ = 0 at the surface, K [L/T] is the hydraulic conductivity at a given volumetric strain, K₀ [L/T] is the hydraulic conductivity at $\varepsilon_v$ = 0 at the surface, n [-] is a power-law exponent that describes the sensitivity of K to changes in $\theta$. A larger n indicates that small changes in $\theta$ will result in more significant changes in K. Conversely, a smaller n (e.g., close to 1) indicates a weaker dependence. n is known to vary from 5 to 25 depending on the type of rock and the geomechanical conditions (Akande et al., 2021; Cappa & Rutqvist, 2011; Laouafa et al., 2021; J. Liu et al., 2019; Zhang et al., 2020). Here, we use n = 15 as it provides an optimal representation of the decreasing trend of K over 3 orders of magnitude at 2 km depth (Figures 1c and 1d), similar to those shown in permeability databases (Achtziger-Zupančič et al., 2017; Manning & Ingebritsen, 1999). Here, we fix $\theta_0$ = 0.01 and K₀ = 0.1 m/d, as representative values for crystalline rocks at the surface (Achtziger-Zupančič et al., 2017). Spatially-variable $\theta$ and K values from the geomechanical model are transferred to the hydrogeological model, which has different discretization characteristics.
Steady-state groundwater flow simulations are performed with a 2D finite difference approach implemented in MODFLOW (Harbaugh, 2005; Niswonger et al., 2011) and handled via the python package FloPy (Bakker et al., 2016). The hillslope geometry is discretized in cells of 5 m along its horizontal length (L in Figure 1a) while the vertical discretization comprises 200 layers, resulting in 200 by 200 mesh cells. Furthermore, to enhance the precision of the processes in proximity to the surface, a vertical discretization was implemented through an exponential increase of the mesh size with depth, starting with a first layer of 10 cm thickness. A drainage boundary condition is activated at the surface whenever the groundwater table intersects the surface, using the drain package (DRN). Furthermore, a recharge boundary condition is implemented (RCH) where a constant value was applied across the spatial domain of the model. The magnitude of the steady-state recharge is adjusted for each hillslope scenario in order to obtain a seepage length of 40% of the total length of the slope, L (Table S1 in Supporting Information S1). This is necessary to be able to compare the scenarios and eliminate potential effects of varying seepage on groundwater flow partitioning and residence time distribution (Bresciani et al., 2014). A homogeneous model, characterized by uniform, isotropic $K_0$ and $\theta_0$, is also executed and serves as reference for the comparisons. We consider the distributed discharge rates across the seepage zone, Q, that we normalize by the total recharge, $R_T$, as the sum of the recharge. This allows the cross-comparison of the different cases where $R_T$ might differ. In addition, we consider the specific Darcy flux, q, that we normalize by the recharge R. This allows us to distinguish domains with high ($q/R > 1$) and low ($q/R < 1$) fluxes.

For transport, particle-tracking simulation is performed using MODPATH (Pollock, 2016). Twenty particles were injected regularly spaced for each cell over the top boundary of the water table. We performed an exhaustive convergence study that revealed no influence on the model results when increasing particle densities. Residence times, t, are calculated at the ending point in seepage zones, and they are normalized by the mean $\tau$ to distinguish between fast ($t/\tau < 1$) and slow ($t/\tau > 1$) flow. To compute the probability density function (PDF), the particles were weighted relative to the local effective recharge (Pollock, 2016). Both flow and transport simulations were conducted under steady-state conditions to focus on dominant flow patterns and to avoid the potential impacts of temporal water table variations, which could impact analysis of stress-dependent hydraulic conductivity on our hydrological variables of interest.

## 3. Results

### 3.1. Hydraulic Conductivity and Porosity Fields Under Different Strain Conditions

$\theta$ and K vary significantly across the hillslope domain in consequence to the LTS (Figure 1d). As expected, lower $\varepsilon_v$ values are observed near the crest, while high values are present at the valley bottom (Figures 1b and S1 in Supporting Information S1). For non-zero LTS, $\varepsilon_v$ and K fields are heterogeneous (Figures 2a and 2b, as well as Figures S2, S3, and S4 in Supporting Information S1). The impact is greater on steeper hillslopes (Figure 2c - cases $C_{50-low}$ and $C_{50-high}$, see Table S1 in Supporting Information S1), resulting in a non-homogeneous distribution of K and $\theta$ along the slope. Near the crest, $\theta$ and K are higher and approach the values of $\theta_0$ and $K_0$ respectively, due to the near-zero volumetric strain.

Near the surface, hydraulic properties are more sensitive to LTS (Figures 1d and S4 in Supporting Information S1). As LTS and slope increase ($C_{50-low}$ and $C_{50-high}$), the coefficient of variation of $\theta$ and K (CV in Table S1 in Supporting Information S1) also increases with sharper differences in values between the crest and valley domains. The variability can span up to an order of magnitude for $\theta$ (Figures S2 in Supporting Information S1) and up to four orders of magnitude for K in case $C_{50-high}$ (Figures 1d and 2c as well as Figures S3 and S4 in Supporting Information S1). Figures 1d and S4 in Supporting Information S1 shows that, on average, K is highest at the surface, where the least volumetric strain is involved, and decreases logarithmically with depth, except for the valley bottom where K may be lower. However, $\theta$ and K found at the surface can vary significantly from valley to crest. For case $C_{30-high}$, $K/K_0$ takes a value of ~1 at the crest while and ~10⁻⁴ in the valley bottom. Along the hillslope (i.e., along $x/L$) two patterns are observed, one at the shallow (i.e., $z/L > -1.0$) and at the deeper (i.e., $z/L < -1.0$) part. At shallow depths, we find that $K/K_0$ increases monotonically from the valley to the crest, for example, $C_{50-high}$ exhibits at the valley $K/K_0(z = 0; x = 0) = 10^{-5}$ and at the crest $K/K_0(z = 0.5; x = L) = 1$. For the deeper part, the trend of $K/K_0$ tends to homogenize along horizontal axis of the hillslope ($x/L$), as reflected by the isolines in $K/K_0$ remaining horizontal at depth $z/L < -1.0$ independently from the slope and LTS configurations.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Subsurface flow characteristics under no-LTS and high-LTS (lateral tectonic strain) scenarios ($C_{50-no}$ and $C_{50-high}$). (a, b) Normalized hydraulic conductivity ($K/K_0$). (c) $K/K_0$ distribution in boxplots for all cases. (d, e) Normalized groundwater fluxes ($q/R$). Purple arrows indicate flow vectors and dashed black lines indicate the division of high ($q/R > 1$) and low ($q/R < 1$) flux domains. (f) Flux channeling indicator ($d_F$ – Equation 3) under different LTS scenarios. (g, h) Normalized residence time ($t/\tau$) of particle flow paths. Dashed black lines separate fast ($t/\tau < 1$) and slow ($t/\tau > 1$) residence times. (i) Flowpath channeling indicator ($d_{RT}$ – Equation 4) under different LTS scenarios. (See Figures S3, S5 and S6 in Supporting Information S1 for all combinations of strain scenarios and hillslope gradients).</figcaption></figure>

## 3.2. Discharge Rate Through the Seepage Extent and Channeling Effect

LTS has a pronounced effect on the reorganization of flow across different hillslope scenarios. This reorganization arises from the varying spatial distributions of $\theta$ and $K$ induced by the different combinations of imposed topographic stress and LTS. These distributions modify the pressure fields and, consequently, the distribution of $q/R$ within the hillslope.

In the no-LTS cases, the domain concerned by high flux velocities, that is, where $q/R > 1$ (red area in Figures 2d and 2e), extends toward the lower part of the slope, with the highest velocities concentrated downslope (Figure 2d and $C_{50-no}$ in Table S1 in Supporting Information S1). In the high-LTS case (Figure 2e, $C_{50-high}$), the region where $q/R > 1$ reduces and its lower extent remains limited to $x/L = 0.18$ (Figure 2e). In this case, the valley domain displays low $\theta/\theta_0$ and $K/K_0$ (around 0.55 [−] and $10^{-4}$ [−], respectively) compared to the upper part of the slope, which prevent the flux to develop in this area. To quantify this migration effect of the high flux domain, we
computed the flux channeling indicator, $d_F$, as the ratio between the low-flux area, ($A_{q/R<1}$) and the high-flux area ($A_{q/R>1}$) (values are shown in Table S1 in Supporting Information S1).

$$d_F = \frac{A_{q/R<1}}{A_{q/R>1}} \qquad (3)$$

As LTS increases, $d_F$ tends to increase across all three slope scenarios (Figure 2f), confirming that the high flux domain shrinks as the contrast in K between shallow and deeper parts of the hillslope intensifies in response to stress. We found a general linear trend between $d_F$ and LTS for all three slope cases in the range of LTS investigated here. Furthermore, we found that the relative change in flux channeling for a given change in LTS is greater as the hillslope angles increases (i.e., the slope of the linear relationships in Figure 2f). For each case, the gradient of the linear trend is 0.4, 14.8 and 24.6 for C₁₀, C₃₀ and C₅₀, respectively (Figure 2f). This latter means that, for example, in the C₅₀ cases, a 0.1% increase in the LTS would increase the ratio between low-flux and high-flux by 2.46%.

LTS has a significant impact on the distribution of Q. Considering the homogeneous scenarios in Figures 3a and 3c (Table S1 in Supporting Information S1), Q/R_T decreases exponentially along seepage extent (x/L). A similar distribution is observed in the no-LTS cases (C₁₀-no and C₅₀-no), where Q/R_T range from 5.8% to 6.8% for C₁₀-no and C₅₀-no, respectively, to 0% at the end of seepage extent. In those cases, θ decreases by 20% from the surface to the bottom of the model, while K decreases by ~1.5 orders of magnitude. When LTS is applied, high Q/R_T progressively migrate toward the upper part of the seepage zone as a consequence of the increasing channeling effect (Figures 3a and 3c). Most significant changes occur at higher slopes due to enhanced flux channeling, where the Q/R_T increases in the section between x/L = 0.1 – 0.4. Considering C₅₀-high, Q/R_T at the bottom of the slope is reduced by 85% while it increases by 250% in the upper part of the seepage extent, that is, at x/L = 0.3 (Figure 3c).

### 3.3. Flowpath and Residence Time Distributions

Increasing LTS causes the high-flux domain to become more channelized, confining the flowpaths closer to the surface and enhancing the discharge in the upper part of the slope. This is due to the intensifying variation of θ and K with depth as LTS increases, resulting in lower average hydraulic properties (Table S1 in Supporting Information S1 shows average θ, K for all cases). Through the analysis of the particle tracking simulations, we found that the mean residence times τ also increases, driven primarily by the changes in K spanning order of magnitudes, while θ changes over only a few factors (Table S1 in Supporting Information S1).

As LTS reorganizes the redistribution of θ and K in the hillslope, groundwater flow paths become concentrated in a shallower part of the hillslope and migrate upslope, where θ and K are higher (Figures 2g and 2h). As well as $d_F$, we computed a flowpath channeling indicator $d_{RT}$ (values are shown in Table S1 in Supporting Information S1). It is computed as:

$$d_{RT} = \frac{A_{t/\tau>1}}{A_{t/\tau<1}} \qquad (4)$$

As expected, $d_{RT}$ increase as LTS increases (Figure 2i), indicating that the fast flow path domain (shown in blue in Figures 2g and 2h) progressively shrinks as increasing LTS is involved. This trend is observed across all slopes and LTS configurations (see Figure S6 in Supporting Information S1 for all cases).

The computed residence time distributions (Figures 3b and 3d) show strong deviations from the widely used exponential model in all cases. This is expected, as the model geometry favors flow path partitioning at depth and heterogeneous effective recharge along the surface, which diverges the assumptions of the exponential distribution model (Haitjema, 1995). Instead, the distribution follows a more complex distribution, with a gamma like shape for fast flows (t/τ < 0.2) and transitioning to a power law distribution at slower flows (t/τ > 0.2). In the homogeneous case, the power law has a value of -1.85.

We observed a general increase in the likelihood of fast residence times (Figures 3b and 3d) as LTS increases, particularly pronounced for steeper slopes. In the low slope scenario, the maximum probability shifts from
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Seepage and residence time distributions for the low and high slope scenarios ($C_{10}$ and $C_{50}$). (a, c) Normalized discharge ($Q/R_T$) along the seepage face. (b, d) Probability density function (PDF) of the normalized residence time, semi-log and log-log (inner plot) scale. The black solid line indicates the pdf of the traditional exponential model (Haitjema, 1995) which serves as reference. The lines indicate the PDFs of a homogeneous model (dashed red line), and the No-LTS (blue line), Low-LTS (orange line) and High-LTS (green line) scenarios.</figcaption></figure>

$t/τ = 4 \cdot 10^{-2}$ for $C_{10-no}$ to $t/τ = 2 \cdot 10^{-2}$ for $C_{10-high}$, increasing the maximum probability 0.5. In the high slope case, this shift is even more pronounced, from $t/τ = 4 \cdot 10^{-2}$ for $C_{50-no}$ to $t/τ = 4 \cdot 10^{-3}$ for $C_{50-high}$ and increasing the maximum probability around 8%. This shift is a result of the shrinking of the fast flow domains and the increased flow velocities, which tends to increase the occurrence of fast flow paths in the upper part of the hillslope. Conversely, the distribution of slow paths shows only minor changes. The power law exponent changes slightly from -1.85 for $C_{50-no}$ to -1.62 for $C_{50-high}$ (log-log plot in Figure 3d).

**4. Discussions and Conclusions**

Our results demonstrate that topography and lateral tectonic stress (LTS) play a critical role in shaping the distribution of $\theta$ and $K$ within hillslopes. These variations arise from differences in volumetric strain between the upper and lower parts of the slope. Similar distributions have been previously identified through principal stress calculations using both analytical solutions (L. Liu & Zoback, 1992; Savage et al., 1985) and numerical methods (Martel & Muller, 2000; Miller & Dunne, 1996; Slim et al., 2015). Our findings show that, as slope gradients increase—particularly under high LTS—there is greater variability in $\theta$ and $K$ from the top to the bottom of the
slope. In high-LTS scenarios, the distribution of θ and K agrees with the observations from critical zone studies using geophysical techniques, such as those conducted in Colorado, Carolina, Maryland (St. Clair et al., 2015) and Idaho, US (Chen et al., 2023; Nielson et al., 2021). This geometry, described as “surface-mirroring” (St. Clair et al., 2015), results in not only a reduced average flux and increased groundwater residence times, but also in a compartmentalization of the hillslope with higher θ and K domains located upslope compared to the valley bottom. It has been postulated that such spatial distribution of θ and K, influenced by changes in geomorphic and mechanical factors, might affect hydrogeological processes at the hillslope scale (Anderson, 2015; St. Clair et al., 2015) but without quantification. This study provides a first evaluation, focusing on groundwater seepage discharge and residence time distributions.

Traditional conceptual models of groundwater flow system suggest that groundwater discharge mostly occurs in the lower part of a landscape, where streams are typically located (Freeze, 1972; Freeze & Witherspoon, 1967; Jiang et al., 2009; Tóth, 1962). However, our findings show that in stressed hillslopes, a significant portion of groundwater discharge may occur upslope, avoiding lower K domains at the valley bottom. Due to the challenges of in-situ distributed measurements, the distribution of discharge rates along stream networks has been poorly evaluated in previous field-based studies. Our results suggest that the main springs network might occur at a higher elevation in the hillslope than was previously considered under traditional conceptual models, especially in steep topography settings.

With regard to the residence time distribution, the results significantly differ from the classical exponential model (Haitjema, 1995) which assumes homogeneous aquifer and groundwater flow paths parallel to topography. Instead, the distributions exhibit a gamma-shaped pattern (Leray et al., 2019) for lower values of $t/τ$ and a power law for higher values. Our findings indicate a higher occurrence of fast flows than predicted by the exponential model, due to channeling effects. Previous studies have observed this increase in probability in fast flows (Ameli et al., 2016; Leray et al., 2019), typically linking it to specific geometric characteristics of the hillslope, such as the distance between the discharge zone and the impermeable base. The combined effects of flow partitioning within the hillslope and stress-dependent heterogeneities in θ and K might also cause this distribution. For the high residence times, we found that the tailing of the distributions follows a power law with exponents ranging from -1.62 to -1.85. Notably, these exponent values are similar to those reported by Cardenas & Jiang, 2010, who also shows that, as heterogeneity in hydraulic properties increases, the power of the tailing distribution decreases due to the increasing contrast between fast (shallow) and slow (deep) flow systems.

It is important to acknowledge that real hillslopes likely exhibit more complex flow patterns than those simulated in this conceptual modeling framework. We have performed an analysis on a real 2D topography (Text S1, Figures S7 and S8 in Supporting Information S1). While the distributions of θ and K show some spatial variability due to the smaller scale roughness of the topography, the general patterns are similar to those of our synthetic cases. Both exhibit high values in the crest and lower values in the valley bottom. Similar to the synthetic cases, the main discharge zone occurs upslope and channeling is enhanced near the surface as LTS increases.

Our 2D synthetic analysis allows us to perform a sensitivity analysis on the role of slope and lateral strain on groundwater flow and residence times. Such synthetic models are particularly valuable for isolating and understanding the fundamental processes at play, especially in complex systems where multiple factors interact. Here, we provide some key guidance on the topographic and tectonic contexts at which such stress-dependent heterogeneity might significantly influence flow and transport processes, representing a potential source of uncertainty for classical hydrological models assuming homogeneous aquifers. Future investigations are needed to explore the potential impacts of other geomechanical processes such as fracturing. Indeed, in our analyses, variations in θ and K are solely derived from volumetric strain of a homogeneous and isotropic geomechanical model, assuming purely elastic deformation. Mechanical properties were adjusted to maintain the medium within the elastic range, excluding more complex, non-elastic behavior, such as fracturing and faulting. Additionally, varying lithologies may also introduce variability in local mechanical properties that may significantly impact the θ and K fields. In addition, our 2D modeling framework restricts the consideration of geomechanical and flow effects that may manifest differently in three-dimension. For instance, the stress concentration observed in the valley bottom of a steep topography may redirect groundwater flow to alternative directions, avoiding the low hydraulic conductivity zones.

These findings build-on our conceptualization of critical zone processes involved at the hillslope scale, particularly relevant for mountainous areas where steep slopes and high tectonic stresses are prominent. Such
partitioning of groundwater flow and residence times will have strong controls on erosion and weathering processes, with direct consequences on landscape evolution and deep critical zone architecture (Harman & Cosans, 2019; Litwin et al., 2022; Riebe et al., 2017). They may be considered when describing the occurrence and dynamics of springs and stream networks in steep landscapes, complementing information on controls due to topography, lithology or geological discontinuities (Costigan et al., 2016; Florianic et al., 2024; Keegan-Treloar et al., 2022; Kirchner et al., 2023; Prancevic & Kirchner, 2019; Shanafield et al., 2021). In such heterogeneous hillslopes, groundwater is expected to exhibit a more complex storage-discharge relationship. While shallow, high-K zones promote faster response timescales, deeper flow paths enable longer response timescales (Roques et al., 2022). This implies a high spatio-temporal dynamic of water availability and solute export for groundwater-dependent aquatic ecosystems (Winter, 1999). In this context, shallow groundwater-dominated ecosystems may show high sensitivity to climate variability (Carroll et al., 2024; Kløve et al., 2014; Taylor et al., 2013), whereas ecosystems connected to deeper flow paths, mostly located downslope are likely concerned to greater resilience through sustained groundwater discharge.

In addition to groundwater flow and water availability, our outcomes have implications on weathering processes and solute export at a catchment scale. Indeed, high *K* and channelized flow in the shallow subsurface are expected to enhance weathering rates by increasing contact between groundwater and reactive minerals, thus facilitating dissolution and transport of solutes (Riebe et al., 2017). High-*K* zones at crests may provide preferential pathways for water with dissolved gases (e.g., O₂ and CO₂), enhancing oxidation reactions locally, while low-*K* zones, that is, valley bottoms, could promote reductive chemical processes. As shown by Bochet et al., 2020, the enhanced mixing between the oxic and anoxic waters induced by channeling could drive local hotspots for chemical reactions. These processes will shape the long-term development of the weathering profiles, ultimately shaping the critical zone and the evolution of the landscape (Rempe & Dietrich, 2014).

To transfer this knowledge to real-world systems requires the integration of interdisciplinary field-based observations. Environmental tracers, geophysical imaging, groundwater seepage mapping and its temporal dynamics, alongside stress information from boreholes, would create the required framework to predict the role of geomechanically constrained subsurface heterogeneity in shaping groundwater flow and solute export in the Critical Zone.

## Data Availability Statement

All geomechanical and groundwater model are publicly available in Figueroa (2024).

## Acknowledgments

We are grateful for the two anonymous reviewers and the editor, Hari Viswanathan, for their valuable suggestions that helped improve the quality of the manuscript. L.J.S.H. acknowledges the SNSF project RADMOGG (Grant 212622). R.A. and C. R. acknowledge the European project WATERLINE (ID CHIST-ERA-19-CES-006).

## References

Achtziger-Zupančič, P., Loew, S., & Mariéthoz, G. (2017). A new global database to improve predictions of permeability distribution in crystalline rocks at site scale. *Journal of Geophysical Research: Solid Earth*, 122(5), 3513–3539. https://doi.org/10.1002/2017JB014106

Akande, W. G., Gan, Q., Cornwell, D. G., & De Siena, L. (2021). Thermo-Hydro-mechanical model and Caprock deformation explain the onset of an ongoing seismo-volcanic unrest. *Journal of Geophysical Research: Solid Earth*, 126(3). https://doi.org/10.1029/2020JB020449

Ameli, A. A., Amvrosiadi, N., Grabs, T., Laudon, H., Creed, I. F., McDonnell, J. J., & Bishop, K. (2016). Hillslope permeability architecture controls on subsurface transit time distribution and flow paths. *Journal of Hydrology*, 543, 17–30. https://doi.org/10.1016/j.jhydrol.2016.04.071

Anderson, R. S. (2015). Pinched topography initiates the critical zone. *Science*, 350(6260), 506–507. https://doi.org/10.1126/science.aad2266

Bakker, M., Post, V., Langevin, C. D., Hughes, J. D., White, J. T., Starn, J. J., & Fienen, M. N. (2016). Scripting MODFLOW model development using Python and FloPy. *Ground Water*, 54(5), 733–739. https://doi.org/10.1111/gwat.12413

Bochet, O., Bethencourt, L., Dufresne, A., Farasin, J., Pédrot, M., Labasque, T., et al. (2020). Iron-oxidizer hotspots formed by intermittent oxic-anoxic fluid mixing in fractured rocks. *Nature Geoscience*, 13(2), 149–155. https://doi.org/10.1038/s41561-019-0509-1

Bresciani, E., Davy, P., & De Dreuzy, J. R. (2014). Is the Dupuit assumption suitable for predicting the groundwater seepage area in hillslopes? *Water Resources Research*, 50(3), 2394–2406. https://doi.org/10.1002/2013WR014284

Cappa, F., & Rutqvist, J. (2011). Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO₂. *International Journal of Greenhouse Gas Control*, 5(2), 336–346. https://doi.org/10.1016/j.ijggc.2010.08.005

Cardenas, M. B., & Jiang, X. W. (2010). Groundwater flow, transport, and residence times through topography-driven basins with exponentially decreasing permeability and porosity. *Water Resources Research*, 46(11). https://doi.org/10.1029/2010WR009370

Carroll, R. W. H., Niswonger, R. G., Ulrich, C., Varadharajan, C., Siirila-Woodburn, E. R., & Williams, K. H. (2024). Declining groundwater storage expected to amplify mountain streamflow reductions in a warmer world. *Nature Water*, 2(5), 419–433. https://doi.org/10.1038/s44221-024-00239-0

Chen, H., Niu, Q., McNamara, J. P., & Flores, A. N. (2024). Influence of subsurface critical zone structure on hydrological partitioning in mountainous headwater catchments. *Geophysical Research Letters*, 51(6). https://doi.org/10.1029/2023GL06964

Chen, H., Niu, Q., Mendieta, A., Bradford, J., & McNamara, J. (2023). Geophysics-informed hydrologic modeling of a mountain headwater catchment for studying hydrological partitioning in the critical zone. *Water Resources Research*, 59(12). https://doi.org/10.1029/2023WR035280
Chin, L. Y., Raghavan, R., & Thomas, L. K. (2000). Fully coupled geomechanics and fluid-flow analysis of wells with stress-dependent permeability. *SPE Journal*, 5(01), 32–45. https://doi.org/10.2118/58968-PA

Costigan, K. H., Jaeger, K. L., Goss, C. W., Fritz, K. M., & Goebel, P. C. (2016). Understanding controls on flow permanence in intermittent rivers to aid ecological research: Integrating meteorology, geology and land cover. In *Ecohydrology* (Vol. 9(7), pp. 1141–1153). John Wiley and Sons Ltd. https://doi.org/10.1002/eco.1712

Fan, Y., Clark, M., Lawrence, D. M., Swenson, S., Band, L. E., Brantley, S. L., et al. (2019). Hillslope hydrology in global change research and earth system modeling. *Water Resources Research*, 55(2), 1737–1772. https://doi.org/10.1029/2018WR023903

Figueroa, R. (2024). Geomechanical and hydrogeological models for different hillslopes and tectonic stresses [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.14014516

Floriancic, M. G., Abhervé, R., Bouchez, C., Jimenez-Martinez, J., & Roques, C. (2024). Evidence of groundwater seepage and mixing at the vicinity of a knickpoint in a mountain stream. *Geophysical Research Letters*, 51(17). https://doi.org/10.1029/2024GL111325

Freeze, R. A. (1972). Role of subsurface flow in generating surface runoff: 1. Base flow contributions to channel flow. *Water Resources Research*, 8(3), 609–623. https://doi.org/10.1029/WR008i003p00609

Freeze, R. A., & Witherspoon, P. A. (1967). Theoretical analysis of regional groundwater flow: 2. Effect of water-table configuration and subsurface permeability variation. *Water Resources Research*, 3(2), 623–634. https://doi.org/10.1029/WR003i002p00623

Gauvain, A., Leray, S., Marçais, J., Roques, C., Vautier, C., Gresselin, F., et al. (2021). Geomorphological controls on groundwater transit times: A synthetic analysis at the hillslope scale. *Water Resources Research*, 57(7). https://doi.org/10.1029/2020WR029463

Goderniaux, P., Davy, P., Bresciani, E., De Dreuzy, J. R., & Le Borgne, T. (2013). Partitioning a regional groundwater flow system into shallow local and deep regional flow compartments. *Water Resources Research*, 49(4), 2274–2286. https://doi.org/10.1002/wrcr.20186

Haitjema, H. M. (1995). On the residence time distribution in idealized groundwatersheds. *Journal of Hydrology*, 172(1–4), 127–146. https://doi.org/10.1016/0022-1694(95)02732-5

Harbaugh, A. W. (2005). MODFLOW-2005: The U.S. Geological survey modular ground-water model: The ground-water flow process. *Techniques and Methods*. https://doi.org/10.3133/tm6A16

Harman, C. J., & Cosans, C. L. (2019). A low-dimensional model of bedrock weathering and lateral flow coevolution in hillslopes: 2. Controls on weathering and permeability profiles, drainage hydraulics, and solute export pathways. *Hydrological Processes*, 33(8), 1168–1190. https://doi.org/10.1002/hyp.13385

Ingebritsen, S. E., & Gleeson, T. (2015). Crustal permeability: Introduction to the special issue. *Geofluids*, 15(1–2), 1–10. https://doi.org/10.1111/gfl.12118

Ingebritsen, S. E., & Manning, C. E. (2010). Permeability of the continental crust: Dynamic variations inferred from seismicity and metamorphism. *Geofluids*, 10(1–2), 193–205. https://doi.org/10.1111/j.1468-8123.2010.00278.x

Jiang, X. W., Wan, L., Wang, X. S., Ge, S., & Liu, J. (2009). Effect of exponential decay in hydraulic conductivity with depth on regional groundwater flow. *Geophysical Research Letters*, 36(24). https://doi.org/10.1029/2009GL041251

Keegan-Treloar, R., Irvine, D. J., Solórzano-Rivas, S. C., Werner, A. D., Banks, E. W., & Currell, M. J. (2022). Fault-Controlled springs: A review. In *Earth-science reviews* (Vol. 230, p. 104058). Elsevier B.V. https://doi.org/10.1016/j.earscirev.2022.104058

Kirchner, J. W., Benettin, P., & van Meerveld, I. (2023). Instructive surprises in the hydrological functioning of landscapes. *Annual Review of Earth and Planetary Sciences*, 51(1), 277–299. https://doi.org/10.1146/annurev-earth-071822-100356

Kløve, B., Ala-Aho, P., Bertrand, G., Gurdak, J. J., Kupfersberger, H., Kvarner, J., et al. (2014). Climate change impacts on groundwater and dependent ecosystems. *Journal of Hydrology*, 518(PB), 250–266. https://doi.org/10.1016/j.jhydrol.2013.06.037

Laouafa, F., Guo, J., & Quintard, M. (2021). Underground rock dissolution and geomechanical issues. *Rock Mechanics and Rock Engineering*, 54(7), 3423–3445. https://doi.org/10.1007/s00603-020-02320-y

Leray, S., Engdahl, N. B., Massoudieh, A., Bresciani, E., & McCallum, J. (2016). Residence time distributions for hydrologic systems: Mechanical foundations and steady-state analytical solutions. *Journal of Hydrology*, 543, 67–87. https://doi.org/10.1016/j.jhydrol.2016.01.068

Leray, S., Gauvain, A., & de Dreuzy, J. R. (2019). Residence time distributions in non-uniform aquifer recharge and thickness conditions – An analytical approach based on the assumption of Dupuit-Forchheimer. *Journal of Hydrology*, 574, 110–128. https://doi.org/10.1016/j.jhydrol.2019.04.032

Litwin, D. G., Tucker, G. E., Barnhart, K. R., & Harman, C. J. (2022). Groundwater affects the geomorphic and hydrologic properties of coevolved landscapes. *Journal of Geophysical Research: Earth Surface*, 127(1). https://doi.org/10.1029/2021JF006239

Liu, H.-H., Wei, M.-Y., & Rutqvist, J. (2013). Normal-stress dependence of fracture hydraulic properties including two-phase flow properties. *Hydrogeology Journal*, 21(2), 371–382. https://doi.org/10.1007/s10040-012-0915-6

Liu, J., Wang, J. G., Gao, F., Leung, C. F., & Ma, Z. (2019). A fully coupled fracture equivalent continuum-dual porosity model for hydro-mechanical process in fractured shale gas reservoirs. *Computers and Geotechnics*, 106, 143–160. https://doi.org/10.1016/j.compgeo.2018.10.017

Liu, L., & Zoback, M. D. (1992). The effect of topography on the state of stress in the crust: Application to the site of the Cajon Pass Scientific Drilling Project. *Journal of Geophysical Research*, 97(B4), 5095–5108. https://doi.org/10.1029/91JB01355

Manning, C. E., & Ingebritsen, S. E. (1999). Permeability of the continental crust: Implications of geothermal data and metamorphic systems. *Reviews of Geophysics*, 37(1), 127–150. https://doi.org/10.1029/1998RG900002

Maréchal, J. C., Dewandel, B., & Subrahmanyam, K. (2004). Use of hydraulic tests at different scales to characterize fracture network properties in the weathered-fractured layer of a hard rock aquifer. *Water Resources Research*, 40(11). https://doi.org/10.1029/2004WR003137

Martel, S. J. (2000). Modeling elastic stresses in long ridges with the displacement discontinuity method. *Pure and Applied Geophysics*, 157(6–8), 1039–1057. https://doi.org/10.1007/s000240050016

Martel, S. J., & Muller, J. R. (2000). A two-dimensional boundary element method for calculating elastic gravitational stresses in slopes. *Pure and Applied Geophysics*, 157(6–8), 989–1007. https://doi.org/10.1007/s000240050014

Masset, O., & Loew, S. (2010). Hydraulic conductivity distribution in crystalline rocks, derived from inflows to tunnels and galleries in the Central Alps, Switzerland. *Hydrogeology Journal*, 18(4), 863–891. https://doi.org/10.1007/s10040-009-0569-1

Miller, D. J., & Dunne, T. (1996). Topographic perturbations of regional stresses and consequent bedrock fracturing. *Journal of Geophysical Research*, 101(B11), 25523–25536. https://doi.org/10.1029/96jb02531

Moon, S., Perron, J. T., Martel, S. J., Holbrook, W. S., & St. Clair, J. (2017). A model of three-dimensional topographic stresses with implications for bedrock fractures, surface processes, and landscape evolution. *Journal of Geophysical Research: Earth Surface*, 122(4), 823–846. https://doi.org/10.1002/2016JF004155

Nielson, T., Bradford, J., Holbrook, W. S., & Seyfried, M. (2021). The effect of aspect and elevation on critical zone architecture in the Reynolds creek critical zone observatory: A seismic refraction study. *Frontiers in Water*, 3. https://doi.org/10.3389/frwa.2021.670524
Niswonger, R. G., Panday, S., Ibaraki, M., & Survey, U. S. G. (2011). MODFLOW-NWT, a Newton formulation for MODFLOW-2005. In *Techniques and methods*. https://doi.org/10.3133/tm6A37

Oda, M., Takemura, T., & Aoki, T. (2002). Damage growth and permeability change in triaxial compression tests of Inada granite. *Mechanics of Materials*, 34(6), 313–331. https://doi.org/10.1016/S0167-6636(02)00115-1

Oda, M., Takemura, T., & Suzuki, K. (2023). Depth-dependent hydraulic conductivity tensor for jointed granite. *International Journal of Rock Mechanics and Mining Sciences*, 168, 105405. https://doi.org/10.1016/j.ijrmms.2023.105405

Oestreicher, N., Lei, Q., Loew, S., & Roques, C. (2023). Bedrock fractures control groundwater-driven mountain slope deformations. *Journal of Geophysical Research: Earth Surface*, 128(11), 105406. https://doi.org/10.1029/2022JF006885

Pollock, D. W. (2016). User guide for MODPATH version 7—A particle-tracking model for MODFLOW. In *Open-file report*. https://doi.org/10.3133/ofr20161086

Prancevic, J. P., & Kirchner, J. W. (2019). Topographic controls on the extension and retraction of flowing streams. *Geophysical Research Letters*, 46(4), 2084–2092. https://doi.org/10.1029/2018GL081799

Rempe, D. M., & Dietrich, W. E. (2014). A bottom-up control on fresh-bedrock topography under landscapes. *Proceedings of the National Academy of Sciences*, 111(18), 6576–6581. https://doi.org/10.1073/pnas.1404763111

Riebe, C. S., Hahm, W. J., & Brantley, S. L. (2017). Controls on deep critical zone architecture: A historical review and four testable hypotheses. *Earth Surface Processes and Landforms*, 42(1), 128–156. https://doi.org/10.1002/esp.4052

Rocscience Inc. (2001). User's guide. 2D finite element program for stresses and estimating support around underground excavations. Retrieved from https://www.rocscience.com/software/rs2

Roques, C., Rupp, D. E., De Dreuzy, J. R., Longuevergne, L., Jachens, E. R., Grant, G., et al. (2022). Recession discharge from compartmentalized bedrock hillslopes. *Hydrology and Earth System Sciences*, 26(16), 4391–4405. https://doi.org/10.5194/hess-26-4391-2022

Rupp, D. E., & Selker, J. S. (2006). On the use of the Boussinesq equation for interpreting recession hydrographs from sloping aquifers. *Water Resources Research*, 42(12), 128–156. https://doi.org/10.1029/2006WR005080

Savage, W. Z., Swolfs, H. S., & Powers, P. S. (1985). Gravitational stresses in long symmetric ridges and valleys. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, 22(5), 291–302. https://doi.org/10.1016/0148-9062(85)92061-3

Shanafield, M., Bourke, S. A., Zimmer, M. A., & Costigan, K. H. (2021). An overview of the hydrology of non-perennial rivers and streams. In *Wiley interdisciplinary reviews: Water* (Vol. 8(2)). John Wiley and Sons Inc. https://doi.org/10.1002/wat2.1504

Slim, M., Perron, J. T., Martel, S. J., & Singha, K. (2015). Topographic stress and rock fracture: A two-dimensional numerical model for arbitrary topography and preliminary comparison with borehole observations. *Earth Surface Processes and Landforms*, 40(4), 512–529. https://doi.org/10.1002/esp.3646

St. Clair, J., Moon, S., Holbrook, W. S., Perron, J. T., Riebe, C. S., Martel, S. J., et al. (2015). Geophysical imaging reveals topographic stress control of bedrock weathering. *Science*, 350(6260), 534–538. https://doi.org/10.1126/science.aab2210

Stober, I., & Bucher, K. (2015). Hydraulic conductivity of fractured upper crust: Insights from hydraulic tests in boreholes and fluid-rock interaction in crystalline basement rocks. *Geofluids*, 15(1–2), 161–178. https://doi.org/10.1111/gfl.12104

Taylor, R. G., Scanlon, B., Döll, P., Rodell, M., Van Beek, R., Wada, Y., et al. (2013). Ground water and climate change. *Nature Climate Change*, 3(4), 322–329. https://doi.org/10.1038/nclimate1744

Tóth, J. (1962). A theory of groundwater motion in small drainage basins in central Alberta, Canada. *Journal of Geophysical Research*, 67(11), 4375–4388. https://doi.org/10.1029/JZ067i011p04375

Winter, T. C. (1999). Relation of streams, lakes, and wetlands to groundwater flow systems. *Hydrogeology Journal*, 7(1), 28–45. https://doi.org/10.1007/s100400050178

Zhang, L., Zhang, H., Yuan, Y., & Yin, S. (2020). Poroelastoplastic reservoir modeling by tangent stiffness matrix method. *Petroleum*, 6(4), 438–450. https://doi.org/10.1016/j.petlm.2020.07.004