Bedrock Fractures Control Groundwater‐Driven Mountain Slope Deformations

Seasonal deformation of mountain rock slopes can be driven by groundwater infiltration and depletion. Such processes could explain our field observation in the Aletsch Valley, Switzerland, where GNSS‐derived 3D annual displacement amplitudes reach 3.4 cm. However, the physical mechanisms behind such groundwater‐driven surface displacements are not well understood. Here, we develop a fully coupled hydromechanical model to simulate the relevant processes in a valley slope embedded with numerous fractures of variable sizes. The magnitude and orientation of transient annual slope surface displacement obtained from our model are in overall agreement with the field observations. The key geological factors controlling the type and magnitude of reversible mountain slope deformations are fracture network geometry, fracture aperture, and regional stress field. We show that the heterogeneity and anisotropy of bedrock hydromechanical responses, originating from depth‐dependent variations of fracture properties, play a critical role in groundwater recharge and valley slope deformation. During recharge events, pore pressure perturbations migrate downward from the groundwater table and toward the receiving stream and the deep subsurface. This process driven by pressure diffusion and poroelastic stressing develops in the subsurface with a great reach of up to a few kilometers, called critical hydromechanical response zone, and controls surface deformation patterns. During groundwater recession, this hydromechanical response zone expands downward and ground surface displacement vectors rotate upwards. Our results suggest that slope surface deformation can inform about subsurface permeability structures and pore pressure fluctuations, which have important implications for understanding groundwater flow in fractured bedrock slopes.

Within an alpine valley slope, groundwater flow is strongly influenced by topography (Gleeson & Manning, 2008;Goderniaux et al., 2013;Welch & Allen, 2012) and bedrock permeability structure that is often characterized by a high-permeability zone at shallow depth (Achtziger-Zupančič et al., 2017;Fu et al., 2022;Manning & Ingebritsen, 1999;Rapp et al., 2020;Welch & Allen, 2014).Moon et al. (2020) showed that stress conditions in the subsurface, influencing the closure and shear dilation of fractures, can explain a zone of high permeability in the near-surface in the first few hundred meters below the surface.In addition, fracture density is depth dependent with a high fracture density in the near-surface accommodating enhanced permeability, as shown in many places (e.g., Achtziger-Zupančič et al., 2017;Carlsson et al., 1983;Moon et al., 2020;Mortimer et al., 2011).However, the consequences of the rapidly decaying rock mass permeability in the shallow subsurface on pore pressure-related surface deformation magnitude and orientation remain unexplored.Due to the coupled hydromechanical processes involved and the multiscale nature of fracture networks associated with heterogeneous properties, a thorough understanding of the underlying physical mechanisms (e.g., poroelasticity, pore pressure diffusion) driving cyclic slope deformations in alpine valleys and their spatio-temporal dynamics remains elusive.

Key Observations in the Aletsch Valley
Improving the understanding of impacts caused by hydromechanical coupling in fracture networks on the spatial patterns of transient ground surface displacement orientations and magnitudes requires comprehensive field observations.Ground surface displacement patterns and structures have been monitored at an ETH Zürich in situ laboratory, located around the tongue of the Great Aletsch Glacier in the Swiss Alps by Glueer et al. (2021) and Oestreicher et al. (2021) since 2012, using continuous Global Navigation Satellite System (cGNSS) stations (Limpach et al., 2016) (Figure 1).The six cGNSS stations (i.e., AL01, AL02, AL03, ALTS, ALTD, and FIES), anchored on the granitic and gneissic bedrock of the Aar massif, exhibit long-term deformational trends and cyclic displacements on top of some random noise (Oestreicher et al., 2021).The observed 3D displacements at different elevations in the valley and topographic profiles, which pass through the six cGNSS stations, are illustrated in Figure 1.These show that most of the annual cyclic displacements captured at the stations, in general larger than a centimeter in amplitude, are constrained within the two-dimensional cross-sections.All stations except ALTD record outward displacement (positive in the longitudinal direction shown in green) from spring to early summer and inward displacement from late summer to winter.The annual cyclic vertical displacement is either small or below the noise level at most stations.This differs at station FIES, where the vertical displacement component 10.1029/2022JF006885 3 of 23 exhibits a significant uplift signal in spring and subsidence in summer and fall.We note that FIES is situated closest to the ridge of the slope (and in another valley), while other stations are closer to the valley bottom or the current glacier elevation (Figures 1a and 1b). Figure 1c shows the daily position of each cGNSS station through the year in the vertical plane of the profiles shown in Figure 1b.The monthly average positions of the cGNSS points (labeled with numbers) follow an annual hysteretic loop, with significant local variations between stations.Stations AL02, FIES, and ALTD are located on South-East facing slopes, while AL01, AL03, and ALTS are on North-West facing slopes.ALTD is the only station not displaying any significant annual cyclic displacement.Station AL02 moves approximately horizontally out of the slope in spring, with a fast motion in April and May.The peak displacement occurs in May.FIES starts its spring displacement in June and peaks in July before returning to its initial position at a lower speed during summer and autumn.Between April and July, FIES moves upward and outward, 2.4 cm.FIES is located in the Fiesch valley at an elevation of ∼2,350 m a.s.l., while AL02 is located in the Aletsch valley at ∼1,950 m a.s.l., where the temperature is higher, and the snow melts earlier in the year.The three stations on the North-West facing slope of the Aletsch valley all show a similar pattern of fast spring displacement in May and June outward, followed by a slower displacement back to its winter position (November-April).At station AL03, the displacement between February and July reaches 2.6 cm.The stations displaying annual cyclic motions have a hysteresis, close to the noise level of the cGNSS instruments outlined by the daily average positions (colored dots in Figure 1).

Research Objectives and Approach
In this work, we investigate the effects of bedrock fractures on the spatio-temporal evolution of slope surface deformations, that is, the displacement magnitudes, orientations and annual hysteresis, associated with transient groundwater flow.We develop a fully coupled hydromechanical model capable of simulating the two-way coupling processes of groundwater flow and rock mass deformation within km-scale valley flanks embedded with a multiscale fracture system and assigned with geologically relevant properties.We focus on alpine recharge cycles (with dominant groundwater recharge from snowmelt in late spring to early summer and groundwater depletion during the rest of the year) and investigate the effects of depth-dependent fracture properties and in situ stress conditions on valley slope surface deformations.This model is tested based on the unique data set from the Great Aletsch Glacier valley, where we dispose of detailed knowledge about multi-annual reversible slope displacements at catchment scale (Oestreicher et al., 2021) as well as fracture patterns and recharge conditions (Hugentobler et al., 2020).The ultimate goal of this study is to advance our understanding of the hydromechanical mechanisms of groundwater-driven rock slope deformations, especially regarding the roles of bedrock fractures and hydromechanical couplings.We also aim to identify the key controlling geological and hydrogeological factors on spatially heterogeneous slope surface displacement patterns.

Numerical Methods and Model Setup
In order to assess the driving mechanisms and relationships between groundwater recharge, fracture properties, and surface displacements, we conduct 2D numerical simulations of coupled hydromechanical processes in an alpine valley.We represent fractured rock slopes in typical alpine mountainous environments (Figure 2) and set up the model using the finite element software Comsol Multiphysics v6.0 (COMSOL AB, 2021).We discretize the model (with a width of about 10 km) into a dense grid of blocks with block size in the hundred meter scale.We model the distribution of mesoscale fractures in crystalline bedrock (i.e., with lengths smaller than the grid block size) based on the discrete fracture network (DFN) approach (Lei et al., 2017) with their impact on the rock mass properties of each grid block modeled through a homogenization treatment (L.Wang & Lei, 2021).Note that we do not attempt to build a model to accurately reproduce the conditions at the Aletsch valley, but rather aim to investigate the problem from a generic perspective and to gain insights into the common physical mechanisms that are applicable to similar crystalline rock slopes.We test this generic model by using its results to qualitatively interpret the field observations from the Aletsch site.

Geomechanics Model
In our model, the mechanical equilibrium of fractured rock is governed by: where σ is the total stress, f b is the body force, and f ext are the external forces such as boundary loads.The stress-strain relationship of the rock material obeys the linear poroelasticity law (Jaeger et al., 2007) given as: where ϵ is the strain, C is the compliance tensor, σ′ is the effective stress, α is the Biot-Willis coefficient, p is the pore pressure, and I is the identity matrix.
The closure behavior of a fracture under normal compression is governed by Bandis et al. (1983) and Barton et al. (1985): where v m is the maximum allowable closure, κ n0 is the initial normal stiffness, and   ′ n is the effective normal stress (i.e.,   ′ n = n −  , where σ n is the normal stress).For a fracture subject to shearing, the shear stress-shear displacement relationship is governed by Coulomb's friction law with the peak shear stress  p =  ′ n tan f , where ϕ f is the friction angle of the fracture.The shear dislocation of the fracture walls τ s could induce dilational displacement v s which can be estimated as (Gan & Elsworth, 2016;Rahman et al., 2002;Rutqvist et al., 2013;Ucar et al., 2017;Willis-Richards et al., 1996): where κ s is the shear stiffness (Barton et al., 1985), and ϕ d is the dilation angle given by Ladanyi and Archambault (1969) and Lei et al. (2020): where ϕ d0 is the initial dilation angle, and σ c is the uniaxial compressive strength of the intact rock.The fracture aperture b f is then calculated as (Lei et al., 2016): where b f0 is the initial aperture of the fracture (at zero stress).
For each grid block embedded with fractures, we compute the compliance tensor as (Oda et al., 1993): where N c is the number of fractures within the given block, κ n and κ s are the normal and shear stiffnesses of the fracture, respectively, l is the fracture length, δ ij is the Kronecker's delta,   mat ijkl is the compliance tensor of the intact rock matrix, F ij and F ijkl are the so-called crack tensors derived as (Oda, 1986;L. Wang & Lei, 2021): where A b is the block area, and n i , n j , n k , and n l are the directional components of each fracture.

Groundwater Flow Model
We compute groundwater flow through deformable fractured rock masses with spatially and temporally varying permeability tensors.The continuity equation for single-phase fluid flow in fractured rocks is given by: where ρ w is the density of water, φ is the porosity, Q s is the source term, ϵ v is the volumetric strain of the rock mass, and q is the flux defined in Darcy's law as: where k is the local permeability tensor and μ w is the dynamic viscosity of water.Furthermore, the storage behavior is governed by Rutqvist and Stephansson (2003): where S is the storage coefficient.The above Equations 10-12 may be reduced to (Rutqvist & Stephansson, 2003): For each grid block, the rock mass permeability tensor is calculated as: where   mat ij is the permeability tensor of the intact rock matrix, and P ij is the crack tensor calculated as:

Model Setup and Calculation Procedures
We construct a representative alpine valley with two symmetrical valley flanks (highlighted by a dotted rectangle in Figure 2), bounded by two additional auxiliary valleys to minimize boundary effects (Text S2 and Figure S1 in Supporting Information S1).The model has an elevation difference between the top and bottom of the slope of 1,000 m and a width from crest to crest of 5,000 m.We consider a fluvial valley approximated as a smooth V shape with a slope angle of ∼21.8° (Figure 2).The model is subject to gravity and has the bottom constrained by a roller boundary condition, while the top is a free surface.The left and right boundaries are subject to a displacement control; zero horizontal displacement corresponds to a roller boundary condition, while a non-zero value defines a prestrained condition linearly increasing with depth (the maximum lateral strain takes a value of 0-1,000 μm/m); vertical displacement at the side and top boundaries is unconstrained.
We assume that there is little lateral flow in the subsurface below the lateral boundaries of the model, and therefore define the bottom as well as the sides of the model as no flow boundaries.At the surface, we model both the recharge process during periods of recharge and the seepage process at locations where the groundwater table reaches the ground surface (Chui & Freyberg, 2009;Grämiger et al., 2020): where n is the normal to the surface, u is the fluid velocity, k m is the permeability of the surface rock, ρ r is the rock density, g is the gravitational force, f inf is the recharge function and f s is the seepage velocity function, χ s and χ u are the smoothing functions for the saturated and unsaturated parts of the surface boundary, respectively.
Material properties are listed in Table 1.The parameters are chosen to represent bedrock consisting of gneiss or granite embedded with small-scale (<5 m) fractures, similar to bedrock found in the Aletsch valley (Berger et al., 2016;Steck, 1983).Properties of the larger fractures (<200 m) are listed in Table 1.In this study, we conduct sensitivity analyses of several key parameters, such as fracture density, aperture, orientation, and lateral confinement, in the ranges given in Table 1, using the one-factor-at-a-time approach (Razavi & Gupta, 2015).
Two different scenarios are explored when generating fracture networks.In the first scenario, we consider a system with a random orientation of fractures, and a depth dependency of fracture density, as observed in many  -Zupančič et al., 2017;Carlsson et al., 1983).The variation of fracture density follows an exponential function as (Jiang et al., 2010): where z is the depth below surface, d min the residual fracture density, d max the maximum fracture density, and ζ is the density decay rate.We define a base case with ζ = 7.4 × 10 −4 1/m, d min /d max = 0.1, and the total number of fractures in the model at 10 5 .Here, the value of ζ adopted (7.4 × 10 −4 1/m), which is calibrated based on the permeability calculation (see Section 4), is somewhat lower than the value of 2.8 × 10 −3 1/m reported in Jiang et al. (2010), due to the different geological systems under consideration.The distribution of fracture lengths follows a power-law, as fractures in rock often exhibit a wide range of lengths following a power-law distribution (Bonnet et al., 2001;Lei & Wang, 2016): where l min is the minimum length set at 5 m, l max the maximum length set at 200 m, and a the power-law exponent, assumed as 2 in the current study, a typical value as reported in the literature (see e.g., Bonnet et al., 2001).In the second scenario, we model a fracture system dominated by a systematic fracture set consisting of subparallel fractures with the same depth-density relationship as in the first scenario (see Equation 17).
We perform a series of simulations to explore the role of fractures on pore pressure variations, groundwater flow, rock mass deformation and surface displacements.The selection of parameters is informed by past studies of fractured rock aquifers, where such properties have been studied in the field (e.g., Masset & Loew, 2010) or from extensive data compilations (Achtziger-Zupančič et al., 2017).For the base case, we populate 10 5 fractures having an isotropic orientation and depth-dependent density (residual density ratio d min /d max of 0.1).The initial aperture of the fractures b 0 is set at 0.15 mm with a maximum allowable closure v m at 0.135 mm, and an initial normal and shear stiffnesses of 25 and 12.5 GPa/m, respectively.To test the importance of depth-dependent fracture density, we run simulations with various residual density ratios, from 0.01 (strong density decrease with depth) to 1 (uniform density in the model, no depth-dependency) (see Section 3.2).We also analyze the effect of fracture anisotropy in a case where there is only one vertical set of fractures, with a random dispersion of the dip angle of 30° and a case with a set of oblique 45° inclined fractures, with a random dispersion of the dip angle of 30° too (Section 3.3).Then, we run a sensitivity analysis on the initial fracture aperture (Section 3.3).We also apply prestrained conditions on the lateral boundaries to explore the effect of higher tectonic confinements (Section 3.4).
We discretize the study domain using an unstructured mesh of triangular elements with a maximum element size of 60 m and a refinement at the top surface (with a size of 5 m).We select mesh element sizes leading to stable results within reasonable run time for a simulation.The adaptive time stepping of the solver is subject to the constraint of a maximum time step of 1 d during periods of recession and 0.1 d during periods of recharge.Each simulation is run in three consecutive stages.In the first stage, the pore pressure and mechanical loads are progressively ramped up in a steady-state manner.In the second stage, the model is run in a transient manner and progressively brought to equilibrium through a total of 9 recharge-depletion cycles.Finally, a last, 10th, recharge-recession cycle is simulated and analyzed, such that the model is at equilibrium and does not exhibit inter-annual trends in pore pressure change or surface deformation.We model a single annual recharge step representing groundwater recharge in response to snowmelt.The other smaller recharge events during the rest of the year are considered by applying a constant minimum recharge rate of a third of the recharge rate during snowmelt.The snowmelt recharge pulse starts at 35% and stops at 60% of the year, corresponding roughly to the period between May and early August in the Northern hemisphere.The maximum recharge rate is taken as ∼1.1 mm/d, corresponding to a total yearly water amount of 0.1 m.Considering a recharge rate of ∼10%, the precipitation needed is around 1 m/y which is a reasonable amount for typical alpine environments (Markovich et al., 2019).

Simulation Results and Analysis
In this section, we first show the simulation results of groundwater table changes, slope surface displacements and their spatio-temporal variations in the base case.We then show the results of the one-factor-at-a-time sensitivity analysis with respect to the base case, where we change each of the parameters (fracture network geometry, 10.1029/2022JF006885 9 of 23 initial fracture aperture and rock mass permeability, and regional confinement) to test their effects on the pore pressure variations and slope surface deformations.

Base Case: Hydromechanical Response of Fractured Aquifers
In the base case, the rock mass has an exponentially decaying fracture density with depth (Figure 3a), varying from around 20 m/m 2 at the near-surface to around 12 m/m 2 in the deep subsurface.Permeability also decreases from a value of around 2 × 10 −14 m 2 in the near-surface to a permeability of ∼1 × 10 −17 m 2 in the deep subsurface (Figure 3b).The principal directions of stress tensors as shown in Figure 3c exhibit a rotation from parallel to the topography (with the horizontal-to-vertical stress ratio between 1 and 5) in the near-surface to vertically oriented in the deep subsurface (with the stress ratio between 0.3 and 1). Figure 3d shows the variation of groundwater table in the rock mass between the low water level period (before the start of the main recharge pulse) and the high water level period (at the end of the main recharge period).The phreatic surfaces of both water level conditions are marked by blue curves; the area that is permanently saturated is in dark blue, while the area within which the water table fluctuates is colored in light blue (Figure 3d).
The displacement field is shown in Figure 3d as well, where black arrows illustrate the orientation and magnitude of ground displacement between the start and end of the recharge period.Deformation is oriented toward the center of the valley and upwards during recharge events.The displacement amplitude at the mountain crest is up to 3.4 and 0.02 cm near the valley bottom.At the mountain crests, deformation is subvertical (Figure 3d).For locations at mid-slope position, we observe a progressive rotation of the displacement vector, oriented toward the center of the valley around the points where the slope is permanently fully saturated, and oriented upwards for locations higher in the slopes (Figure 3d).
Displacements change over space and time in response to the recharge-related pore pressure front diffusing from the top of the mountain crests toward the subsurface region and laterally to the valley center, where discharge occurs (Figure 4).During recharge, most deformation is concentrated in the top region of the model, and surface displacement is greater in the top half of the slope than near the valley bottom (Figure 4a).The orientation of displacement is vertically upwards close to the crests, but rotates to horizontal around the bottom of the slope and is oriented downwards at the center of the valley.At the end of the recharge period, surface displacement reaches its maximum values (Figure 4b).Points in the bottom half of the slope start to rotate and point upwards.We see a shift of the deforming zone of the rock mass to greater depth (>2 km below mountain crests), together with diffusion of the pore pressure front.After the end of recharge, the pore pressure front diffuses to the deepest parts of the model (up to 3 km below the crests) and laterally to the center of the valley (Figure 4c).The differential pore pressure is rapidly reduced.However, the pore pressure difference is preserved at the end of recharge in the bottom part of the system.During the recession period, surface displacement magnitude is strongly reduced, and the displacement direction progressively rotates upward (Figure 4d).
Figure 5 shows displacement data at the three control points P1, P2, and P3, for a full cycle of recharge and recession.Displacement at the crest (P1) and valley bottom (P3) is predominantly vertical.However, we identify minor hysteresis caused by a small random asymmetry in the fracture network on both sides of the control points (Figure 5).Hysteresis is greater at the mid-slope (P2), with displacement of the control point out of the slope during recharge and an upward rotation during recession (Figure 5).

Effect of Fracture Network Geometry
The model variant with a uniform fracture density distribution (i.e., without depth dependency) shows a marked decrease of rock mass permeability with depth too, varying from around 5 × 10 −14 m 2 near the surface to 1 × 10 −16 m 2 at depth, caused by the variation of stresses with depth.In comparison to the base case, the model with uniform fracture density has larger permeability at depth.This results in smaller (∼100 m) groundwater table variations during a recharge-recession cycle (Figure S8 in Supporting Information S1), but the deforming part extends further down below the mountain crests and reaches the bottom of the model at the end of the recharge (i.e., 3 km below the crests).The smaller pore pressure fluctuations observed in the rock mass lead to smaller magnitude seasonal deformations.Maximum displacement observed during recharge for the model with a uniform fracture density distribution is 6 mm, while it is 36 mm for a model with density ratio d min /d max = 0.1, and 42 mm for a model with d min /d max = 0.05 (Figure S8 in Supporting Information S1).The density decrease with depth influences the orientation of displacement too.We observe that vertical displacement is more sensitive to variations of the density ratio than horizontal displacement, as captured by points in the middle of the slopes.This effect is more pronounced for a density ratio less than 0.2, where the changes in displacement magnitudes are also greater.For a density ratio above 0.2, the horizontal displacement magnitude is consistently around 1.5 times larger than the vertical (Figure S8 in Supporting Information S1).
In conditions where a single fracture set dominates the fracture system (e.g., Figures S2-S5 in Supporting Information S1), the orientations of hydraulically active fractures result in anisotropy of permeability.This leads to smaller groundwater table elevation changes, less pore pressure variations, and higher groundwater table elevations in the slope.In parallel, the anisotropic elasticity of the rock mass affects the slope's displacement.With a set of vertical fractures, as shown in Figure S2 of the Supporting Information S1, there is a zone (∼200 m deep below the valley bottom and ∼1,200 m deep below the crests) where the slope deforms around 6%-12% more in the vertical direction than in the horizontal direction.

Effect of Fracture Aperture
We evaluate the impact of initial and residual fracture aperture and stiffness, which exert major control on the rock mass permeability (e.g., for the initial aperture variation of 0.1-1 mm, the permeability varies by around three orders of magnitude; Figure S6 in Supporting Information S1).The increase in permeability is inhomogeneous in space and is particularly pronounced close to the surface, where permeability is greater (Figure S7a in Supporting Information S1).Normalizing the permeability change shows that the largest relative permeability difference is at depth below the valley, and that permeability change is minimal near the mountain crests (Figure S7b in Supporting Information S1). Figure S7 in Supporting Information S1 shows that the relationship between the initial fracture aperture and the rock mass permeability is not trivial.To better compare the results of the models with different initial fracture aperture values, we analyze the time series of pore pressure and displacement at the five control points as marked in Figure 2.
In Figure 6, we show that the groundwater table follows a similar trend for all control points if the initial aperture is large (e.g., 1 mm), indicating an increase during recharge followed by a gradual decrease during the recession period.The maximum groundwater level is reached at the end of recharge (red dashed line), and the minimum level is reached at its start (blue dashed line).Decreasing the initial fracture aperture (and thus, the rock mass permeability) strongly affects the groundwater table and its fluctuations.We observe that the groundwater table elevation difference between control points at a given time increases with reduced fracture aperture.
Horizontal displacement shown in Figures 6f-6k is negligible in a system with large initial fracture aperture.The average change in hydraulic gradient between low and high water conditions is relatively small (0.09 for an initial fracture aperture of 1 mm).With the decrease of fracture aperture, we find that horizontal displacement becomes greater (>1.5 cm) for points in the middle of the slope (P2 and P4) while it remains small (<0.1 cm) for points at the crests (P1 and P5) and valley bottom (P3).The hydraulic gradient also varies more between low and high water conditions (0.12 for an initial fracture aperture of 0.15 mm).Displacement between the start and end of the recharge period is positive at P2 and negative at P4, that is, oriented toward the center of the valley (Figures 6g and 6i).Vertical displacement shows uplift at all control points during recharge for models with large initial fracture aperture (Figures 6l-6p), but the uplift is reduced in systems with small initial fracture aperture at the mid-slope and below.At the valley bottom (P3), for a particularly small initial fracture aperture (say below 0.2 mm), the effect is reversed, and the control point exhibits minor subsidence during the recharge period.
In Figure 7a, the groundwater table elevation at the end of the recharge (red dashed line in Figure 6) shows an increase for systems with lower initial fracture aperture.Pore pressure changes (Figure 7b) reach a maximum of 250 m below the mountain crests for the scenario with an initial fracture aperture of 0.2 mm.For scenarios with a larger initial aperture, the groundwater table changes converge toward a value of 100 m; for smaller initial apertures, the magnitude of changes decreases at all control points (Figure 7b).Horizontal and vertical displacements are shown in Figures 7c and 7d.Most of the horizontal deformation occurs over a small range of initial fracture apertures.At low values (when the groundwater head changes decrease in Figure 7b), both the horizontal and vertical components of deformation are negligible (<0.05 cm).At large initial fracture apertures (when the groundwater head changes converge to 100 m, and the hydraulic gradient is small, say <30 m, in Figure 7a), the orientation of the surface displacement is subvertical, and the entire valley is uplifted during recharge.

Effect of Horizontal Confinement and Stress Levels
We investigate the effect of horizontal stress levels and stress ratios by progressively increasing horizontal strain at the lateral model boundaries (from a base case with fixed roller boundaries to 400 and 800 μm/m prestrained  boundaries).For these modified boundary conditions, the stress distribution and regime in the rock mass undergo substantial changes from a normal faulting to a strike-slip/thrusting regime, and strong changes of the ratio between the horizontal and vertical stress components occur (Figure 8).We observe a decrease in permeability (particularly in the valley center) for models with increased lateral confinement compared to the base case.
Groundwater pressure changes during recharge decrease when confinement (i.e., horizontal stress) is added, with a larger change for small strain values (Figure 9b).The observed displacement decreases, with the vertical component more affected compared to the horizontal one (Figures 9c and 9d).Therefore, the ratio between the horizontal and vertical displacements increases with increasing lateral confinement (Figure 9e).Also, the ratio between displacement magnitude and groundwater head change generally decreases with increasing lateral strain and stress.A system under higher confinement (strike slipe or thrust faulting), therefore, produces less surface deformation for a given head change.

Mechanisms of Groundwater-Driven Slope Deformation
In Section 3, we showed (a) how displacement signals monitored at the ground surface relate to groundwater table variations and pore pressure changes in the subsurface, and (b) how fracture properties like density, orientation and aperture as well as regional stress conditions control the groundwater flow and surface displacement dynamics in alpine mountains.Transient ground surface displacements follow a hysteretic path related to the diffusion of pore pressure front following a recharge event.The pressure front diffuses from the main recharge location at the mountain crest toward the seepage zone in the lower part of the slope.In our model, the typical depth range of this hydromechanical response zone is about two times the valley relief (say up to 2-3 km into the subsurface rock), reaching substantially below the valley bottom.Thus, the surface displacement orientation and magnitude depend on the location and shape of the hydromechanical response zone, and evolve through time (Figure 4).Our results indicate that the key factor controlling the shape and size of the hydromechanical response zone, and thus the direction and amplitude of surface displacements, is the subsurface rock mass permeability and anisotropy.
Several studies have defined a "hydraulically active region" with increased permeability and thus groundwater flow down to 200-500 m depth (Gleeson & Manning, 2008;Markovich et al., 2019;Ofterdinger et al., 2014;Welch & Allen, 2012).We would like to emphasize the fundamental difference between the hydraulically active region defined in the previous studies and the hydromechanical response zone identified in our study.The former is dominated by pressure diffusion and thus characterized by a shallow depth (i.e., less than 500 m), while the latter is additionally boosted by poroelastic stressing which can reach beyond the region of pressure diffusion  and is therefore characterized by a great depth (i.e., up to a few kilometers).The presence of this poroelastic stressing effect can be justified through an estimation of the diffusion timescale using the simple 1D diffusion equation t = z 2 /(2D) for water vertically reaching a depth of z = 2 km, where D is the hydraulic diffusivity (taking a geometric mean of 0.06 m 2 /s for the upper 2 km in our model).This gives a time of ∼400 days, which is much larger than the deformational response timescale of ∼100 days or less as observed in our model.This poroelastic stressing effect is mainly caused by the groundwater recharge-induced extra gravitational forces that load the deep subsurface rock via stress transfer, which can also accommodate pressure perturbations in subsurface rock at great depth subject to a partially drained condition.
Clearly, the capture of these processes requires fully coupled hydromechanical simulation as conducted in our current work.Similar phenomena that poroelastic stressing can cause deformations beyond the pressure diffusion zone have been reported in many studies on injection-induced earthquakes (Atkinson et al., 2020;Goebel et al., 2017;Schultz et al., 2020;Segall & Lu, 2015) and tunneling-induced ground deformations (Zhao et al., 2023).
In our model, we imposed a decrease in permeability with depth which results from the combined effect of increasing stress (Figure 10) and decreasing fracture density (Achtziger-Zupančič et al., 2017;Manning & Ingebritsen, 1999).Furthermore, the change of permeability depends on topography, with a stronger permeability decrease with depth below the valley (2.7 × 10 −17 m 2 at 2 km depth in Figure 10) than below the mountain crests (1.3 × 10 −16 m 2 at 2 km depth in Figure 10), due to local topography-driven stress variations, for example, Figure 3 (Achtziger-Zupančič et al., 2017;Fu et al., 2022;Manning & Ingebritsen, 1999;Moon et al., 2020;Rapp et al., 2020;Welch & Allen, 2014).Permeability is strongly affected by local stress conditions, which control fracture normal closure and shear dilation behavior (Lei et al., 2016(Lei et al., , 2017)).The depth-permeability relationship exerts a strong control on the shape and extent of the hydromechanical response zone.A stronger decrease in permeability with depth results in substantially larger surface deformations due to larger groundwater table variations, but also greater horizontal deformation at the mid-slope (P2 and P4 in Figure S8 of the Supporting Information S1).
Low permeability of less than 1 × 10 −15 m 2 (with initial fracture aperture <0.09 mm in Figure 7) in the near-surface region (say <150 m) keeps groundwater close to the surface and prevents large pore pressure changes at depth.This permeability is typical for aquitards where potential recharge is greater than permeability, such as in intact crystalline rocks or sediments with high clay/shale contents.High permeability of more than 1 × 10 −13 m 2 (e.g., with initial fracture aperture >0.3 mm in Figure 7) in the near-surface region drives rapid groundwater flow at depth and a flatter groundwater table, associated with strong vertical deformation and little horizontal displacement at the surface.Significant variations in phreatic groundwater level elevation are observed only between these two regimes, that is, in a narrow permeability range close to the magnitude of annual recharge.Here, the magnitudes of induced ground surface displacements are observable by conventional monitoring systems, and the hydromechanical response zone takes the shape of an elongated zone below the mountain crests (Figure 4b) with the pressure front diffusing downwards and laterally during and after recharge.Our model shows that the magnitude of groundwater table variations is not linearly related to permeability changes in slopes, as exemplified by varying the initial aperture for fractures in the model.The change of regime is relatively sharp at an initial fracture aperture between ∼0.08 and 0.3 mm (Figure 7b) and leads to larger groundwater table variations and larger horizontal displacement at the mid-slope.
For models with a preferred fracture orientation, we observed strong anisotropy in both the rock mass permeability and surface displacement pattern (see Section 3.2).Such observations of strong anisotropy in permeability with preferred orientation of fractures are consistent with small-scale (≤10 m) studies (Noorian Bidgoli & Jing, 2014;Ren et al., 2015;Tsang et al., 2007).We can formulate two different mechanisms to explain the observations.First, anisotropic deformation of the rock mass may be solely due to the anisotropic rock mass elasticity.Indeed, the presence of discontinuities strongly impacts the rock mass deformation modulus (Hoek & Diederichs, 2006), and the organization of the fracture network dominated by a single set naturally results in an  2017), and the dashed black line is the equation from the data set of Manning and Ingebritsen (1999).Equations link permeability k, in m 2 , with depth z, in km.Variability in permeability between valley flanks is caused by local effects of fracture network geometry.anisotropic modulus (Barton, 2013).This is because the shear stiffness of fractures is typically lower than their normal stiffness (Bandis et al., 1983).Second, anisotropic deformation of the rock mass may be the result of hydromechanical forcing.In this case, the presence of the fracture set induces anisotropy in permeability (e.g., Rutqvist & Stephansson, 2003), modifying the hydromechanical response zone at depth.The change in the shape of the hydromechanical response zone then impacts the rock mass deformation field and observed displacement at the surface, together with the anisotropic elasticity.To test which of the two mechanisms is dominant, we run the model with a set of vertical fractures (Figure S2 in Supporting Information S1), but imposing the same pore pressure variations as in the base case (Figure 3). Figure 11a shows the deformation during recharge for the model with random fracture orientation and can be compared with deformation in the model with vertical fractures (Figure 11b).The difference is less than a millimeter, while there is a significant increase in deformation (up to 0.8 cm) for the fully coupled hydromechanical model with a set of vertical fractures (Figure 11c).The effect of anisotropic permeability leads to greater deformation below the crest of the mountains (where recharge occurs) and at depth, and slightly smaller deformation in the lower half of the slope (Figure 11).Thus, our model suggests that the primary impact of the preferential fracture orientation on the valley-scale deformation is the redistribution of groundwater flow in the subsurface, rather than the anisotropic modulus of the rock mass.
The impact of fracture normal stiffness-even when varied over a broad range-on horizontal and vertical displacements is relatively minor (Figure S9 in Supporting Information S1) and smaller than that of fracture closure (Figure S10 in Supporting Information S1) and stress ratio (Figure 9).Therefore, we expect significant recharge-related surface deformations in a wide range of tectonic settings, as long as the rock mass permeability lies in the critical range (approximately 1 × 10 −16 -1 × 10 −14 m 2 at shallow depth for an annual precipitation rate around 1 m/y), typically observed for fractured crystalline rocks in Alpine settings (Achtziger-Zupančič et al., 2017;Masset & Loew, 2010;Ofterdinger et al., 2014).
Our results suggest that we can possibly infer information about the rock mass anisotropy in permeability and extent of the hydromechanical response zone from distributed surface displacement observations, that are much easier to access than having to install deep piezometers in alpine rock slopes.The direction of surface displacement following recharge enlightens the dominant fracture set for groundwater flow, that is not always possible to determine from fracture mapping only.Site-specific models including more detailed geometry of the slope surface and fracture network could inform about local stress conditions.
Previous studies showed that larger transmissive fault zones may exert strong impact on groundwater-related slope deformation, for example, in the Aletsch valley (Oestreicher et al., 2021).In the current modeling study, we showed that hydraulically active fractures have a strong control on groundwater flow and surface deformation.Fault zones are at much larger scale than the mesoscale fractures modeled here, and can significantly modify the nearby hydraulic and mechanical conditions.We expect that some fault zones contain densely fractured damage zones and may act as a conduit, allowing fast and channelized pressure diffusion to depth.Such channelized effects along discrete features will introduce singular deformation patterns that depend on fault orientation with respect to the slope and the contrast in hydraulic and mechanical properties with the surrounding host rock.The impact of large-scale fault zones, together with mesoscale fractures, should be investigated in future studies.

Interpretation of Field Observations at Aletsch
The hydromechanical model parameter ranges explored in Section 3 cover many observations at our Aletsch study site.There, subvertical joints parallel to the Alpine foliation form the primary set of persistent and closely spaced fractures, dipping to NW (Figure 12) (see also Grämiger et al., 2017).These subvertical fractures are suggested to strongly influence the infiltration and flow of groundwater and participate in rock mass deformation (Oestreicher et al., 2021).At the Aletsch study site, the valley slopes deform annually by about 1-3 cm in response to snowmelt-related recharge occurring during spring (Oestreicher et al., 2021), and groundwater flows through joints and faults in granites and gneisses (Figure 1 and Oestreicher et al. (2021)).The results of our model scenario dominated by vertical fractures exhibit a similar order of magnitude of ground surface displacement (i.e., 1-3 cm).
Most of the cGNSS stations are situated close to the valley bottom or the current glacier margin, and the orientation of displacement during recharge in spring is subhorizontal.In our model, subhorizontal displacements are also observed for points situated low on the slope in conditions where the shape of the hydromechanical response zone follows the topography at depth, conditioned by a large permeability gradient in the near subsurface.The station FIES is higher on the slope and shows a greater plunge angle during recharge (Figure 1), similar to our model calculations (Figure S2d in Supporting Information S1).This could indicate a deeper groundwater table below FIES, with a larger vertical pressure-change gradient and hence significant vertical displacement components at the ground surface.Spring line mapping in the area of the crest South-West of FIES confirms the deep elevation of the groundwater table and the importance of local recharge (Alpiger, 2013).
The station ALTD exhibits very small annual cyclic displacements, suggesting a different hydromechanical situation than that at the other cGNSS stations in the valley.According to our model, small ground surface displacement occurs in conditions when pore pressure variations at depth are small, for example, at low permeability and with a groundwater table close to the surface (Figure 7).ALTD is situated behind the head scarp of the Driest instability (Glueer et al., 2021), on stable gneissic bedrock.In the upper slope above ALTD, the Drietsch glacier provides water to a surface stream (Drieschtbach), which could limit the groundwater table fluctuations in slopes nearby by contributing to its recharge.In this case, reduced pore pressure changes during snowmelt may explain the absence of cyclic displacement of the station.
The observations of annual hysteresis in the cGNSS stations correspond to that described in our model during the diffusion of pore pressure in the subsurface rock following recharge from the snowmelt.All stations have a relatively constant position during winter (November-March, see Figure 1).Then the stations move outward from the slope during recharge with a plunge angle corresponding to the station's elevation, similar to that shown in Figure 3d.During recession, cGNSS stations move back to their winter position.At the start, they tend to move horizontally and then vertically, inducing the counter-clockwise hysteresis pattern shown in Figure 1.Such a hysteresis is expected for a pore pressure plume migrating to depth, as shown in Figures 4 and 5.
One feature that is present at Aletsch but not included in our simulations is the glacier in the valley bottom, whose dynamics could influence the slope behavior too.As shown in Figure 1, the glacier at the Aletsch site only covers small parts of the lowest slope sectors.In addition, as shown by Hugentobler et al. (2022), the superposition of englacial water pressure fluctuations on the annual slope pressures corresponds to a minor signal the glacier-slope interface.The elastic effect of the snow on the ground, surface water pounding, changes in glacier load due to glacier advance/retreat, and thermoelasticity of the near-surface rock layer are factors that could influence elastic slope surface displacements too.However, Oestreicher et al. (2021) showed that these factors have a minor impact on slope surface displacements at Aletsch and thus we did not include them in our current model.
Note that the generic model presented in this paper does not aim for a quantitative comparison or validation against the observations at the Aletsch site.For a quantitative comparison, a site-specific study is needed, which would have to include small-scale topographic surface roughness variations around monitoring stations, major changes in slope angle as observed on the left valley flank, and large-scale fault zones that can influence recharge, groundwater flow, and rock mass deformation.This will be explored in our future work.Clearly, the generic model presented in the current paper still provides important insights into the underlying mechanisms of observed groundwater-driven slope deformation at Aletsch and the key controlling geological/hydrogeological factors.

Conclusions
We analyzed the role of fractures in groundwater-driven deformation of an alpine valley using a fully coupled hydromechanical model.We showed that fracture-induced heterogeneity and anisotropy significantly influence groundwater flow and ground deformation in valley slopes.Fractures affect ground deformations mostly through their influence on rock mass permeability structure, while groundwater flow modifies the region of the slopes that deforms, that is, the hydromechanical response zone as a result of groundwater recharge-related pressure diffusion and poroelastic stressing.For a slope with height around 1 km and a V-shape, and fractured crystalline rock mass, the hydromechanical response zone is situated in the middle part of the mountain, below the crest, and extends up to a few kilometers into the subsurface.
Our findings bring new insights into understanding groundwater flow and related slope deformations in alpine mountainous environments.We show that after long periods of recharge, such as spring snowmelt in alpine regions, the hydromechanical response zone expands downwards and triggers a hysteretic rotation of the displacement vectors at the surface.Large horizontal surface displacements following recharge occur only at mid-slope for shallow rock mass permeability around 1 × 10 −16 -1 × 10 −14 m 2 with an annual precipitation rate around 1 m/y, and only with a significant gradient of fracture density decrease with depth.Finally, we demonstrate that the presence of preferential fracture orientation has a strong control on the anisotropy of rock mass permeability.However, the mechanical effect of the preferential orientation of fractures is minor.The slope surface displacement patterns are mostly influenced by the anisotropy in pore pressure diffusion.We suggest that measured surface deformations also inform the permeability structure (e.g., the orientation of main hydraulically active fractures) and pore pressure fluctuations at depth, valuable information for understanding groundwater flow in fractured alpine mountain rock slopes.Further studies could include applying our model to detailed alpine valley profiles, including fault zones, and applying the model presented to other study sites with similar groundwater-driven slope deformation.This project is funded by the Swiss National Science Foundation (project 172492).Q. L. is grateful for the support of the Swiss National Science Foundation (Grant 189882) and the National Natural Science Foundation of China (Grant 41961134032).We thank Nadir Dazzi for help tracing lineaments and during fieldwork, Franziska Glueer and Reto Seifert for installing most of the monitoring system, and Philippe Limpach for the processing of the GNSS data.We thank the many field assistants who gave their time for this project through the years.There would not be enough space here to name them all.We thank Jasmin Maissen, Bernadette "Berni" Würsch and Corinne Bischof, as well as colleagues and friends for their support and all the great discussions.We also appreciate the very detailed and constructive review comments from Jeffrey Moore, Anne Voigtländer and one anonymous reviewer.

Figure 1 .
Figure 1.(a) Topography of the study site showing the 6 cGNSS stations (yellow triangles) transected by (b) the 6 topographic profiles (green lines).The Great Aletsch Glacier is shown in blue; the dashed line beneath the ice is the extrapolated bedrock surface to the center of the valley; the positions of the cGNSS stations are at the intersection between the red and green arrows.(c) 3D displacement data from the cGNSS stations, with red being the vertical displacement, V, blue being the horizontal in-plane displacement (in the downhill direction, Hi), and black being the out-of-plane displacement (toward the reader, Ho) in cm.Blue periods are mid-March to mid-June when snowmelt is generally stronger.(d) Annual motion in the V-H cross-section in cm.Monthly average positions are points with a black outline and connected with a black arrow pointing to June, while colored points are daily averages.

Figure 2 .
Figure 2. Schematic model setup for the alpine valley simulation.The model is discretized into a mesh of grid blocks assigned equivalent properties taking into account the distribution of meso-scale fractures.The central dashed rectangle marks the study domain.P1 to P5 are control points used in following figures to compare displacement between model realizations.

Figure 3 .
Figure 3. Model with random fracture orientations and depth-dependent fracture density.(a) Generated discrete fracture network.(b) Permeability field, with permeability ellipses scaled with the logarithm of the permeability magnitude.(c) Stress modeled in the valley (the length and orientation of each line represent the magnitude and orientation of principal stresses, respectively).(d) Observed change in groundwater table before and after the spring snowmelt recharge and corresponding displacement observed in the rock mass between these two instants.

Figure 4 .
Figure 4. Selected times for comparison with low water table conditions before recharge begins.The first row shows the displacement field relative to the start of the recharge (black arrows, whose length is the magnitude of the displacement) and the groundwater table (continuous blue lines).The second row shows the pressure head difference relative to the low water table conditions.The pore pressure plume is in purple (dashed line is 20 m contour, delimited to help the reader visualize the propagation of the pore pressure front).The last row shows the corresponding deformation field in (μm/m).

Figure 5 .
Figure 5. Displacement during a cycle of recharge and recession at the three control points: P1, at crest position; P2, at mid-slope; and P3, at the valley bottom, see Figure 2 for location.The start and end of the recharge are indicated by the blue and red bars respectively.Both scales are logarithmic, and the data are transformed by removing the minimum value of the time series and adding 1 cm.

Figure 7 .
Figure 7. Effect of initial fracture aperture on groundwater table fluctuations and displacement at the five control points from Figure 2 between start and end of recharge (red and blue dashed lines in Figure 6).(a) Groundwatertable elevation in high water table conditions; (b) groundwater table elevation changes between start and end of recharge; (c) horizontal and (d) vertical displacement between start and end of recharge.

Figure 8 .
Figure 8. Stress ratio at depth below the control points from Figure 2, for cases of an additional lateral strain of 0 μm/m (a), 4,000 μm/m (b), and 8,000 μm/m (c), simulating increasing horizontal stress levels.

Figure 9 .
Figure 9.Effect of lateral confinement on groundwater table fluctuations and displacement at the five control points from Figure 2 between start and end of recharge.(a) Groundwatertable elevation in high water table conditions; (b) groundwater table elevation changes between start and end of recharge; (c) horizontal and (d) vertical displacement between start and end of recharge; (e) ratio of horizontal and vertical displacement in (c) and (d); (f) displacement magnitude (from c to d) over groundwater table elevation changes (b) in logarithmic scale.

Figure 10 .
Figure 10.Vertical permeability profiles below the five points displayed in Figure 2. The continuous black line is the equation from field measurements of permeability in crystalline rocks worldwide from Achtziger-Zupančič et al. (2017), and the dashed black line is the equation from the data set ofManning and Ingebritsen (1999).Equations link permeability k, in m 2 , with depth z, in km.Variability in permeability between valley flanks is caused by local effects of fracture network geometry.

Figure 11 .
Figure 11.Comparison between recharge-induced deformation for three different models.(a) Base case, with random orientation of fractures (Figure 3); (b) set of vertical fractures for mechanical model with imposed pore pressure variations from (a); (c) fully coupled hydromechanical model with set of vertical fractures (Figure S2 in Supporting Information S1).Colors and contours are deformation magnitudes, arrows show the deformation directions too.

Figure 12 .
Figure 12.(a) Geologic map of the study area, modified from Steck (2022) with lineaments longer than 200 m (red), the orientation of Alpine foliation with dip angle, and cGNSS stations (yellow triangles); (b) the rose diagram shows the distribution of the lineaments' orientation; (c) the histogram shows the lineaments' length (note the logarithmic scale).

Table 1 Material
Properties of Groundwater and Fractured Rocks geological settings worldwide (Achtziger