Variations in Bedrock and Vegetation Cover Modulate Subsurface Water Flow Dynamics of a Mountainous Hillslope

Predicting the hydrological response of watersheds to climate disturbances requires a detailed understanding of the processes connecting hillslopes and streams. Using a network of soil moisture and temperature sensors, electrical resistivity tomography monitoring, and a weather station we assess the above and below‐ground processes driving the hydrological response of a hillslope during snowmelt and summer monsoon. The transect covers bedrock and vegetation gradients, with a steep upper part characterized by shallow bedrock, and gentle lower part underlain by colluvium. The main vegetation cover is conifers on the upper, and grass and veratrum on the lower part. Combined with a simplified hydrological model, we show that the thin soil layer of the steep slope acts as a preferential flow path, leading to mostly shallow lateral flow, interrupted by vertical flow, mostly at tree locations, and likely facilitated by flow along fractures and roots. Vertical flow and upstream‐driven groundwater dynamics are prevailing at the colluvium, presenting a very different hydrological behavior compared to the upper part. These results show that subsurface structure and features have a strong control on the hydrological response of a hillslope and that those can create considerably varying hydrological dynamics across small spatial scales.


Introduction
Predicting and managing the hydrological response of a watershed critically depends on an appropriate understanding of the climatic forcing and the hydrological connectivity of hillslopes and streams (Blume & van Meerveld, 2015;Bracken et al., 2013).While climatic forcing can be constrained through detailed weather observations, for example, high resolution ground-based radar installations (Sokol et al., 2021), understanding the subsurface characteristics and processes that control how water flows in the subsurface remains a methodological and modeling challenge, especially for non-diffusive processes (Beven & Binley, 2014), yet remains critical for predicting the hydrological response of watersheds to climate disturbances (Hartmann et al., 2017;Schreiner-McGraw & Ajami, 2020).Particularly in mountainous catchments that provide essential input to the water supply of lowland basins (Viviroli et al., 2007), subsurface characteristics have been shown to pose strong control on the discharge characteristics (Apurv & Cai, 2020;Carroll et al., 2018;Rapp et al., 2020;Somers & McKenzie, 2020), and on plant health and dynamics (M.G. Anderson & Ferree, 2010;Callahan et al., 2022;Jiménez-Rodríguez et al., 2022).
Subsurface flow pathways of mountainous hillslopes have been investigated for decades, yet understanding the controls and contributions of various flow pathways remains challenging.Early studies concluded that shallow subsurface flow in the soil layer is the predominant process controlling hillslope hydrology (Mosley, 1979).More recent studies show that deeper flow through the bedrock may provide a considerable contribution to the overall hydrological response.Yet, the results are site dependent and various processes control local and basin scale responses (Frisbee et al., 2017;García-Gamero et al., 2021).For example, preferential flow within the unsaturated zone comprising soil and weathered bedrock is known to assert strong control on the connectivity between hillslope and valley (A.E. Anderson et al., 2009), while deep storage has been shown to provide means to mitigate hydrological droughts (Miguez-Macho & Fan, 2021;Nardini et al., 2021).Bedrock fissures are known to provide preferential conditions for root growth and water infiltration (Stothoff et al., 1999), and bedrock controlled water storage and nutrient supply have been shown to control ecosystem performance and drought resilience (Callahan et al., 2022;Hahm et al., 2022;McCormick et al., 2021).Yet, there still remains a lack of detailed observations of hillslope hydrological processes of mountainous watersheds, and in particular on the interplay between variable soil and bedrock units, and how vegetation and subsurface hydrology may co-evolve, for example, how root growth may contribute to bedrock weathering and thereby alter the hydrological properties (Hasenmueller et al., 2017).
While tracer methods and borehole sampling can provide discrete and direct measurements of flow processes, they usually only provide an estimate of conditions at the sampling location, and upscaling to understand flow processes can be difficult (Blume & van Meerveld, 2015;Jackisch et al., 2017).Similar difficulties arise from capturing the hydrological inputs of mountainous hillslopes with heterogeneity in snowpack dynamics, vegetation cover, and convergence of surface and subsurface flow.Thus, complementing point-based measurements with geophysical methods can improve imaging of the hydrological processes taking place within hillslopes (Robinson et al., 2008), by providing spatially continuous estimates of subsurface properties and processes (Binley et al., 2015;Blume & van Meerveld, 2015).This is particularly important for understanding subsurface flow in heterogeneous conditions.Surface and borehole geophysical measurements have been used to estimate the water holding capacity of the critical zone, distinguishing storage in saprolite, and fractured and unfractured bedrock and their spatial variability (Flinchum et al., 2018).For shallow, fractured bedrock, timelapse 3D GPR measurements have shown that fracture characteristics determine the subsurface hydrological response of steep, soil-mantled hillslopes (Guo et al., 2019), and that when combined with tracer methods hydrologically important subsurface structures and their activation can be determined (Jackisch et al., 2017).Peskett et al. (2020) highlight the benefit of geoelectrical monitoring, which is sensitive to variations in soil moisture content, pore water resistivity, and temperature (Archie, 1942), and combined it with a network of soil moisture sensors to extrapolate findings from discrete, point-based measurements.Their results show that forest strips have a significant impact on the hillslope hydrology by lowering the groundwater table and reducing soil moisture, compared to adjacent grassland.Electrical resistivity tomography (ERT) monitoring also allows to assess soil moisture dynamics below single trees in detail (Fäth et al., 2022), and Rieder and Kneisel (2023) show that individual tree species show different subsurface water demands and dynamics.
Previous studies within the East River watershed (Hubbard et al., 2018), which is representative of many headwater catchments of the Rocky Mountains, have shown that watershed functioning is controlled by hillslope-scale processes that are expressed in terms of their vegetation cover, and their geomorphological and subsurface characteristics (Wainwright et al., 2022).Large-scale fracture zones have been identified that are expected to pose a strong control on the watersheds hydrological response (Miltenberger et al., 2021;Uhlemann et al., 2022).Following a system-of-systems approach and going to smaller granularity, here we focus on hillslope-scale hydrological processes, with a particular emphasis on the differences induced by varying bedrock type and vegetation cover.We aim to answer two questions: how do subsurface flow pathways differ between fracture-dominated and porous subsurface environments?And how do trees access and potentially redistribute water in water-limited, shallow bedrock environments?To provide insights on these questions, we equipped a hillslope characterized by strong variations in bedrock and vegetation types with environmental, point-based, and geophysical sensing, and studied the subsurface response at hourly to daily temporal resolution over more than a year.Our observations are in agreement with other studies, yet provide detailed insights on subsurface flow characteristics that are modulated by the variability in bedrock type and vegetation cover.Thus, the results shown here provide further understanding on mountainous hillslope hydrology that should be considered when assessing and predicting the hydrological response of catchments to climate disturbances.

Study Site
The study site is located at Chicken Bone meadow, on a south-east facing slope in the middle elevations of Mt.Snodgrass (38.92584, −106.97998,WGS84, Figure 1), at about 3,150 m elevation, close to Crested Butte, CO.Mt.Snodgrass is a 3,397 m tall tertiary laccolith, formed by intrusion of granitic magma from the lower crust into the Mancos Shale dominated watersheds of the East River and Washington gulch (Gaskill et al., 1991).Its flanks (between 2,900 and 3,200 m) are covered mostly by Aspen trees (Populus Tremuloides), while the high elevations, like our study site, are covered by spruce-fir and lodgepole pine communities (Langenheim, 1962).The south-eastern side of the mountain is characterized by landslide features, with borehole records indicating landslide deposits varying in thickness between 10 and 30 m. Airborne geophysical data and geological mapping indicate that these landslides usually occur at the interface between quartz monzonite and Mancos shale (Gaskill et al., 1991;Zamudio et al., 2021).

10.1029/2023WR036137
3 of 19 At the Chicken Bone meadow (south-eastern part of Mt.Snodgrass) we established a 235 m long monitoring transect, which was selected to cross the interface between two geological units; in the lower elevations (3,135-3,150 m) colluvium overlays Mancos shale, and quartz monzonite porphyry characterizes the upper part of the transect (elevations of 3,150-3,215 m).Soil sampling showed thin soils (0-1 m thick) with no significant variability across the site.Based on grain size analysis, the soils can be classified as sandy loam, with an average clay content of 13.7% ± 3.7%, a silt content of 18.6% ± 4.9%, and a sand content of 39.0% ± 8.2%.The fraction of coarse grained particles (>2 mm) shows more variability, accounting for 28.7% ± 15.0% of the analyzed samples, with higher coarse content in the western side of the transect.The study site features a change in vegetation, with pine trees (predominantly pinus flexilis) covering the quartz monzonite porphyry dominated part of the transect, and meadow plants and veratrum covering the lower part, except for the end of the profile where spruce trees (mostly picea engelmannii) are present (Falco et al., 2020; Figure 1).The canopy height, which was used to determine the location of trees, is based on high resolution LiDAR data (Falco et al., 2020) and shows distinct trees with canopy heights of up to 20 m at the upper part of the transect and at its lower end.The transect follows the general flow direction as defined by the surface topography and ends at a perennial stream (usually referred to as "Unnamed drainage"), that feeds into Washington Gulch, which is a tributary to the Slate River that drains in the East River, the Gunnison River, and eventually flows into the Colorado River.

Monitoring Installations
Monitoring instrumentation was deployed to observe both above-ground and subsurface dynamics and processes.The below ground monitoring includes an electrical resistivity tomography (ERT) profile, a 30 m deep borehole equipped with a water depth, temperature, and conductivity sensor, five shallow piezometers, and six stations with depth resolved soil moisture and temperature sensors.
The ERT monitoring profile comprises 112 electrodes, with 3 m spacing between electrodes 1 and 49 (0-144 m on the steep, upper part [west] of the profile) and 1.5 m spacing from electrode 50 to the end (at the lower, eastern part of the profile, Figure 2).Data were acquired using a solar-powered, 10 W PRIME ERT monitoring system (Chambers et al., 2022), measuring dipole-dipole measurements with dipole lengths n of 1-7 electrode spacings, and dipole spacings a of 1 to 7n.The survey protocol comprised a total of 4,446 electrical resistance measurements with a full set of reciprocal measurements to estimate the measurements errors.With this electrode layout and survey design, we are sensitive to variations in subsurface properties in the upper 1 m (Knödel et al., 2007), and are able to image to about 35 m depth.We estimated the depth-of-investigation using the sensitivity distribution, and removed all parts of the model that had a sensitivity smaller than 10 −6 .Surveys were conducted daily  (Falco et al., 2020).Coordinates are given in UTM Zone 13N, aerial photograph was obtained from Google Earth.
starting 11 June 2021 until 22 February 2022, from which date survey frequency was increased to twice a day (until 2 March 2022), then four times a day (until 28 March 2022), reaching six times a day from 29 March to 31 May 2022.From 1 June 2022, frequency was reduced to four times a day, and to daily from 10 October 2022.The higher measurement frequency during the spring and summer time was chosen to be able to capture short term dynamics related to snowmelt and summer monsoon processes.ERT data were filtered to remove negative values and measurements with reciprocal errors >25%.Time-lapse inversion of the data, and the optimization function that needs to be minimized, requires that every survey has the same number of measurements (Johnson et al., 2015).Hence, measurements that were removed from the data set were interpolated using an inverse distance approach to account for varying sampling in time; interpolated data were assigned a high error of 100%; details on the temporal variation of measurement errors and the number of removed data points can be found in Supporting Information S1.Data were inverted using E4D (Johnson et al., 2010), with a spatial to temporal constraints ratio of 2, putting larger constraints on smooth spatial variability.Inversions converged to a χ 2 = 1, indicating that the inverted resistivity model fits the data within their errors.
A borehole, SNI1B, was drilled in September 2021 to a depth of 30 m below ground level (bgl).The hole was completed with PVC casing and screened from 30 to 28.5 m depth.Geophysical borehole logging was conducted in the borehole to obtain fluid electrical conductivity and temperature, and nuclear magnetic resonance (NMR) to estimate water content and hydraulic conductivity.The principles of NMR measurements are explained in detail in Dunn et al. (2002).Here we want to highlight that NMR measurements are non-invasive measurements that exploit the magnetic moments of hydrogen atoms to estimate total water content.Based on relaxation times, the distribution of water into clay and capillary-bound, and free water can be estimated, and the relaxation times can be used to estimate hydraulic conductivities (Knight et al., 2016).An In-Situ Inc. Aqua TROLL 200 data logger was installed on 12 November 2021 and recorded water level, temperature, and conductivity at 30 min intervals.Additionally, 5 shallow piezometers (up to 3 m depth) were installed in October 2021, with four being located within the lower, eastern part (BH-A to D), and one in the upper, western part of the transect (BH-E).
At six locations along the transect (at 18, 62, 124, 148, 195, and 231 m, Figure 1) soil moisture, temperature and bulk EC were measured at 0.1, 0.3, 0.6, and 0.8 m depth at 15 min sampling interval (using TEROS 12, Meter) and averaged to provide hourly data (Uhlemann et al., 2024).These depths were chosen to capture the shallow soil moisture dynamics within the soil layer.Deeper dynamics were captured using the piezometers.Co-located with the soil moisture sensors (except for x = 124 m), additional, high-resolution, vertically resolved temperature sensors (Dafflon et al., 2022) were installed measuring soil temperatures to a depth of up to 1.0 m bgl, and air/snow temperatures up to 1.8 m above ground.These sensors measure temperature at 5-10 cm intervals with an accuracy of ±0.075°C.Above-ground temperatures were processed to obtain snow thickness (Dafflon et al., 2022), and subsurface temperatures were used to fit the diffusive heat equation to model seasonal ground temperature variations throughout the transect, from the surface to the bottom of the resistivity modeling domain.Since the electrical resistivity is sensitive to temperature variations, in order to use it to assess soil moisture variations, temperature effects have to be removed, and hence shallow and deep ground temperature variations need to be known.The temperature model can be described as a function of depth z and day of year t by Brunet et al. (2010) and Chambers et al. (2014): with the mean annual ground surface temperature T mean , the amplitude of the annual temperature variation ΔT, the characteristic depth d at which ΔT decreased by 1/e, and a phase offset Φ to ensure that surface and air temperatures are in phase.Note that this equation assumes a sinusoidal annual variation of the ground surface temperature, and is hence not accounting for snow insulation effects, which lead to a relatively constant ground surface temperature during the winter months.To account for spatially varying subsurface thermal parameters, we fitted Equation 1 to each sensor location and interpolated the estimated subsurface temperature to obtain a spatially resolved subsurface temperature model to correct the ERT data for seasonal temperature variations.A ratio model (Hayashi, 2004;Ma et al., 2011) was used to correct the inverted resistivity models ρ to a mean annual temperature of T mean = 6.5°C.The temperature corrected resistivity is defined as employing a temperature correction factor c t = 2.0°C (Hayashi, 2004).Note that we used modeled temperatures for the entire ERT model domain from shallow to deep depths.We also tested whether using the actual measured temperatures at depths ≤1.0 m bgl provides a considerable improvement (i.e., a clear reduction of temperature related shallow resistivity variations), but we found only negligible improvements.Given the expected variability in daily shallow subsurface temperature variations across the transect, the modeled temperatures provide a first order approximation that we deem sufficient for correcting the inverted resistivity data.
Weather and environmental conditions were recorded every 30 min by a multi-sensor station located about 30 m north of the monitoring transect.This station recorded mean horizontal wind speed, precipitation, atmospheric pressure, and relative humidity (ATMOS 41, Meter).A four-component net radiometer (SN-500, Apogee) was installed 1.9 m above the soil surface.This net radiometer includes two pyranometers and two pyrgeometers to measure the upward and downward radiation, respectively, from 0.3 to 3 μm and from 4.5 to 100 μm; this data were used to calculate net radiation.A Campbell Scientific HFP01SC sensor was used to measure the soil heat flux.Snow thickness was measured using a Campbell Scientific SR50AT sonic distance sensor, mounted 1.9 m above ground.The data from this station were used to estimate potential evapotranspiration (PET) based on the FAO-56 Penman-Monteith equation (Allen et al., 1998).

Coupled Hydrological Modeling
To assess the hydrological dynamics of the study site, we created a subsurface flow model using PFLOTRAN (Hammond et al., 2012).The domain was discretized using a refined unstructured mesh based on the discretization used for the ERT inversion, on which Richards equation was solved to model unsaturated flow conditions.We assigned a hydrostatic pressure at an elevation of 3,140 m as initial condition, and used the annual precipitation and PET time-series, that is net infiltration (P-ET), of 2022 as a spatially constant, cyclic boundary condition at the surface of the model, and assigned a constant liquid pressure at the surface nodes representing the stream toward the end of the model domain.Detailed measurements of local precipitation and PET only started in summer 2021, and hence we use the calculated net infiltration of 2022 as cyclic boundary condition, intrinsically assuming that 2022 is a representative year.Given that we want to understand general subsurface flow characteristics, variations in net infiltration from 2022 are expected to not change these general observations.Also note that root uptake was not explicitely modeled, but is represented through the surface boundary condition.Van Genuchten parameters typical for the soils and bedrock present at the site were taken from literature values (Ghanbarian-Alavijeh et al., 2010).The model was run for 10 years, with a maximum time step of one day.The first 5 years acted as spin up (Schäfer et al., 2023), where groundwater levels are increasing to constant seasonal dynamics from year 6 onward.The calculated flow velocities of years 6-10 were used to perform particle tracking, with particles being released at 0.5 m intervals across the surface of the modeling domain.
For each time step, the velocity vectors and the time sampling were used to propagate the particle to an updated position.This was repeated until all time steps from year 6 onward were considered.We then calculated the mean particle flow path for 1 year.
Since there is no control on the potential input from uphill locations, which is not known, this model reflects a simplified and by no means realistic representation of the subsurface hydrological conditions.Hence, we are focusing on the near surface flow pathways, which we assume to be represented appropriately by the surface boundary conditions, that is, net infiltration, which is known relatively accurately through measurements of precipitation, snow accumulation, and atmospheric parameters controlling evapotranspiration.We note that groundwater conditions and dynamics, as observed at borehole SNI1B are captured qualitatively by the model.

Characterization of the Study Site
The baseline ERT data (acquired on 10 October 2021) and borehole geophysical logging data provide detailed insights into the subsurface structure of the study site (Figure 2).The upper, western part of the transect shows comparably high resistivities ranging from 250 to more than 1,300 Ωm.The highest resistivities are recorded close to the start of the transect, where quartz monzonite porphyry is outcropping, and no trees are growing.Further downslope (x = 25-90 m), shallow resistivities show high variability, with low resistivity anomalies occurring at the vicinity of large trees.From x ≈ 100 m a shallow low resistivity anomaly with values around 250 Ωm extends to larger depths and thins out toward the end of the transect.This layer of intermediate resistivities overlays a less resistive layer (<150 Ωm).Drilling records and geophysical borehole logging define the upper intermediate resistivity layer as colluvium, and the lower, less resistive layer as Mancos shale.Additional ERT profiles acquired at the interface between the Mancos shale and quartz monzonite porhpyry at Chicken Bone meadow show similar resistivity pattern that indicate potential landsliding at the interface between the two bedrock units.Whether Chicken Bone meadow is formed of landslide deposits or glacial colluvium is still debated (Gaskill et al., 1991;U.S. Forest Service, 2006), but our data suggests that paleo-landslide deposits formed of weathered quartz monzonite porphyry are likely overlaying the Mancos shale at this meadow.
The geophysical borehole logging data is shown in Figures 2b-2d.The fluid temperature decreases from the water table to about 7 m depth (showing the influence of seasonal temperature variations), and then generally increases with depth, representing the local geothermal gradient.More interestingly, the fluid electrical conductivity shows a considerable increase from about 100 to 125 µS/cm at about 16 m depth, which is within the Mancos shale unit.This may be indicative of a change in groundwater contribution, where the upper part likely has a higher contribution of less electrically conductive snow melt, and hence presents a higher connectivity with surface water, while this seems to be reduced at deeper depths.We derived estimates of the hydraulic conductivity and the water content from the NMR logging data (Figures 2c and 2d).The hydraulic conductivity shows values in the order of 10 −6 -10 −5 m/s in the upper 5 m, which are in agreement with shallow (0.5 m) constant head infiltration test that were conducted throughout Chicken Bone meadow.At about 10 m depth, the hydraulic conductivity decreases to values <10 −8 m/s, before rising to values >10 −6 m/s at depths of more than 15 m.Thus, this layer of low hydraulic conductivity may limit surface water-groundwater interactions with deeper depths.The comparably high hydraulic conductivity values at depth may be indicative of a fractured Mancos Shale, which has been observed in other parts of the East River Watershed (Uhlemann et al., 2022), and geological mapping has identified extensive fracture zones particularly in the vicinity of the granitic intrusions of Mt.Crested Butte and Mt.Gothic (Gaskill et al., 1991).Hence, fractured Mancos Shale is likely also present close to the intrusion of Mt.Snodgrass.
Characteristic values of the simplified heat equation (Equation 1) were obtained by fitting the subsurface temperature monitoring data at six locations along the monitoring transect (Figure 3).The mean annual surface temperature, as well as the data-derived temperature at 0.1 and 0.8 m depth show considerably higher values at the western (Figure 3a), steep portions of the transect.The aspect and change in slope angle result in a reduced direct exposure to solar radiation on this steeper part (see hillshading in Figure 3), thus the higher temperature is driven by a variety of processes.The western part is expected to have larger thermal conductivity and storage due to the presence of shallow quartz monzonite porphyry bedrock than the colluvium of the eastern part, where annual temperatures are more than 1°C colder.In this eastern part, where groundwater levels are shallower and soil moisture is higher, sustained subsurface flow of cold melt water, and elevated evapotranspiration may contribute to lower summer soil temperatures (Alkhaier et al., 2012).The dominating vegetation type, that is, veratrum, is also known to cause lower soil temperatures due to strong shading of the soil (Dafflon et al., 2023).Similarly, tree shading, and possibly different subsurface characteristics, such as variations in the hydraulic properties at the tree locations in the center of the western part reduce the ground temperatures both close to the surface and at shallow depths.A larger amount of leaf litter, limiting evaporation and shading the ground surface may also contribute to these lower temperatures.This shows that a range of processes can contribute to the lower observed temperatures in the gentle and steep parts of the profile, and more data is required to determine the contribution of each.The fitted characteristic depth d, at which the annual temperature anomaly decreases by 1/e, shows similar trends (Figure 3b); the western part has considerably higher values than the eastern part, and because d is independent of the amplitude of the temperature anomaly this indicates a larger thermal conductivity in the western part of the transect.

Overview of Environmental Monitoring Data
Precipitation, snowpack, air temperature, PET, and groundwater depth, resistivity, and temperature for a period from June 2021 to December 2022 are shown in Figure 4.The data highlight the seasonal characteristics of the study site, with warm summer temperatures relating to elevated PET, while cold winters relate to small PET.Precipitation is characterized by intense summer rainfall events with daily maxima of >10 mm, and a snowpack of >1 m covering the surface from mid December to early May.Groundwater levels at the end of summer are around −5 m bgl, and reach two peaks at −2.5 and −2.0 m in early summer.These peaks of the groundwater levels are associated with a drop in water temperature from 4.45°C to 3.95°C and 4.28°C, respectively.Similarly, the rising groundwater in May 2022 is associated with a decreasing groundwater resistivity from 115 to 102 Ωm, while for the second peak no change was observed.Due to the decreasing temperature and groundwater resistivity and its timing only days after snow melted at the study site, we are confident that the first peak is caused by relatively local snowmelt.For the second rise in groundwater levels, the decreasing temperature also suggests snowmelt as a major contributing factor.We note here, that a second peak that is originating from a single input event is often associated to a delayed response from subsurface flow processes (McMillan, 2020); yet the timing and hydrological dynamics, which will be discussed in the following, suggest that the two peaks originate from separate events.

Monitoring of Snowmelt Event
In the following we focus on a single snowmelt event and analyze the change in resistivity and associated changes in soil moisture and temperature (Figure 5).During the monitoring period, snowmelt commenced on 26 March 2022, as can be seen by the entire snowpack reaching 0°C and the shallow soil moisture sensor at 0.1 m depth showing increasing values (Figures 5b and 5c), and lasted until 16 May 2022.The start of snowmelt occurred at The ERT monitoring data allow to assess the subsurface dynamics across the entire monitoring transect (Figure 5a; note that the entire data set is shown in Data Set S1).From 26 March to 14 April changes are highly localized and heterogeneous, with some areas showing increasing and some decreasing resistivities.At TMC 2, where the soil moisture sensors showed increasing values, a decrease in resistivity can be observed.At TMC 5, where no changes in soil moisture below 0.1 m depth were recorded during this period, the resistivity shows no notable change.With progressing snowmelt (after 14 April), we observe a decrease in resistivity at the location of TMC 2, indicating a progressing downward wetting front.Similarly, around TMC 5 a decrease in resistivity can be observed, which agrees with the increase in moisture content that was recorded from 20 April.At the end of snowmelt, 16 May, almost the entire transect shows shallow decreasing resistivity, indicative of a progressing wetting front, except for a positive change in resistivity, at a profile distance of about 120 m.
We also observe a subhorizontal feature of decreasing resistivity, which becomes particularly evident at the change image of 16 May, appearing between x = 140 m and the end of the transect.Additionally, images not shown here, indicate that this feature is moving from deeper to shallower depths.We associate this to a rise in groundwater level, as recorded in the water levels of SNI1B and shown in Figure 5.
The data show some distinct spatial patterns when analyzing the variability of mean resistivity changes of the upper 1 m across the monitoring transect (Figure 6).Isolated areas in the eastern part show negative changes in resistivity from 26 March.From 20 April, the eastern part of the transect shows spatially continuous decrease, first progressing from the surface to deeper layers, and later with an additional change from deeper to shallower layers.In contrast, the western part of the transect shows a more heterogeneous response with isolated areas initially showing mostly negligible changes in resistivity, followed by notably decreasing resistivities from 20 April (Figure 6e).The areas of the western transect not covered with trees show the first response to snowmelt on 20 April, while the area covered by trees show decreasing resistivity more than 2 weeks later on 6 May.We associate this delay to the shading of the trees, and hence a lower incoming radiation, and a lower mean and maximum temperature within the forested part of the transect, and hence slower snowmelt, which is confirmed by lower shallow subsurface temperatures at TMC5 compared to TMC6.At this time almost the entire upper, western part shows a decreasing trend with notably more heterogeneity (i.e., areas of strong negative change next to areas of negligible change) compared to the eastern part of the transect.We associate this pattern with a spatially homogeneous vertical infiltration of snowmelt in the eastern part, and mostly overland, or rapid flow along the soil-bedrock interface in the western part, with deeper infiltration in areas with higher fracture density, where fractures and roots may provide conduits for moisture to reach into the deeper subsurface.This preferential flow is also expected to create a more disconnected infiltration pattern.This can also be seen by the soil moisture response, which shows a flashy behavior with rapid rise and fall of moisture content, particular at depths of 0.6 and 0.8 m, with no response at the shallower sensors (see Figure 6c).The observed behavior in the western part of the transect is in agreement with observations at other sites dominated by shallow, fractured bedrock (García-Gamero et al., 2021;Guo et al., 2019;Mosley, 1979), with rapid lateral flow along the soil-bedrock interface, and percolation along vertical fractures and roots (Ghestem et al., 2011).
During the most intense snowmelt period, between 18 and 23 April, and 5 and 12 May, both the soil moisture sensors and the change in resistivity show clear daily variations.Considering a shallow depth profile at TMC 2 (Figure 7), we can see that peaks in soil moisture are reached almost daily between 18 and 22 April at 18:00 hr at 0.1 m depth, and progress to the deepest sensor at 0.8 m depth 6 hr later.This progressive wetting can also be seen in the change in resistivity, which shows a daily pattern of increasing amplitude with depth and time.This is particularly evident for the snowmelt period in April, but the pattern repeats itself also in May, before being overprinted by upwelling groundwater conditions.This is clear from the resistivity data, which shows a large amplitude negative change progressing from depths of >2 m to the surface, and is also reflected in the soil moisture data, which shows the same trend.This shows that snowmelt has a clear daily pattern, where most of the water seems to infiltrate into the subsurface in the late afternoon within only a few hours; at 0.1 m depth we observe a rise in soil moisture from 16:00 to 18:00 hr, followed by slow decrease in soil moisture until rising again the following day at 16:00 hr.This late afternoon peak is caused by cooling of the snowpack during the night time, which must be overcome before melting can continue (Woelber et al., 2018).

Subsurface Dynamics During Summer Monsoon
To better understand the origin of the second peak observed in the groundwater levels, we analyzed the subsurface dynamics between 15 June and 31 August 2022 (Figure 8; note that the entire data set is shown in Data Set S1). 15 June was used as a reference, as it is just prior to a prolonged period with frequent precipitation events exceeding 10 mm/day (Figure 8e).The changes in resistivity show distinct pattern.First, the shallow subsurface generally shows increasing resistivities (red colors) progressing from the very near surface to deeper layers, indicating a loss of soil moisture and hence a progressive drying front, despite the frequent rainfall events.This indicates that during this period the received rainfall is not sufficient to overcome evapotranspiration and thus provides no recharge to the system, leading to an effective loss of moisture in the shallow soil layers.This has also been concluded from multi-depth soil moisture measurements, obtained at a location down-gradient of the study site (Tokunaga et al., 2019), and from stable isotope data obtained across nine sub-catchments of the East River (Sprenger et al., 2022), indicating water-limited conditions.
The lower, eastern part of the transect is split into two domains, where between x = 180 and 210 m the otherwise dominant surface drying is interrupted at the surface by a subhorizontal feature of decreasing resistivity that extends at depth from x = 140 m to the end of the transect.While the shallow increase in resistivity is indicative of the shallow drying of the soil, the decreasing resistivity shows an upwelling of deeper groundwater related to the second groundwater peak observed in SNI1B.The location of this upwelling coincides with the part of the profile that is covered by veratrum.Outside of the veratrum the vegetation is mostly grass, which usually dies back in early summer, while the veratrum remains healthy and the soil remains at high moisture content for a longer period of time.This can also be seen from the moisture contents recorded at TMC 1 (dashed lines in Figure 8f), which are consistently lower than at TMC 2 throughout the summer.These observations highlight that in the eastern part of the transect, the plant dynamics show more dependence on the snowmelt and related groundwater upwelling processes than on the summer monsoon rainfall, which is consistent with observations at other sites in the East River Watershed that show that the seasonal dynamics of veratrum and other small plants are controlled by snowmelt dynamics (Dafflon et al., 2023).Also note that it is well known that veratrum preferentially occupies wet areas (Langenheim, 1962).
In the upper, western part changes are negligible until early August, at which time decreasing soil moisture trends become prevailing (Figure 8g), particularly at depths of up to 5 m (Figure 8d).Outside the tree covered area of the western part of the transect, resistivities increase earlier, highlighting the shading effect of the trees that reduces evaporation, but also a likely effect of increased surface plant litter in the tree covered area that increase infiltration and lower evaporation (Deutsch et al., 2010), and may store water providing a buffering effect to the underlying soil (Guevara-Escobar et al., 2007).Another explanation could also be that capillary rise at the soil-bedrock interface is transferring water from the weathered bedrock into the soil domain, thereby providing another means of buffering (Kitajima et al., 2013).Trees may also contribute to a more constant soil moisture regime due to a buffering effect and the redistribution of water through the roots tapping into rock moisture.This has been shown to be the case in water-limited environments, where surface soils are dryer than deeper layers during the summer (see, e.g., TMC6, Figure 8h).In such environments, trees take up water from deeper layers, and their root networks aid in the redistribution of water from deep to shallow layers (Kitajima et al., 2013;Prieto et al., 2012).Rempe and Dietrich (2018) have shown that for fractured bedrock rock moisture (water held in weathered bedrock) can provide water storage that may provide plant-available water during prolonged dry conditions, thus buffering the hydrologic response.Rock moisture has been shown to provide more than 50% of plant-available water in shallow bedrock areas, particularly in the Western US (McCormick et al., 2021).Note that in the western part of the transect soil is usually less than 0.5 m thick, and hence weathered and fractured bedrock are close to the surface, and thus our data is representative of such bedrock dynamics.
While during snowmelt the decrease in resistivity related to groundwater upwelling extended toward the end of the transect, and hence the stream, in late summer, this feature ends when entering the forested area at x = 220 m.We interpret this as a reduction in soil moisture due to water uptake by the trees at this location (Figure 8f).Significant reduction in soil moisture within forested areas have also been observed in other studies with soil and bedrock conditions comparable to the eastern part of the transect (Peskett et al., 2020;Rieder & Kneisel, 2023).
Soil moisture recordings and depth-resolved changes in resistivity at the soil moisture monitoring locations show good agreement (Figures 8f-8h).While TMC 2, located in the eastern part of the transect, shows increasing soil moisture and decreasing resistivity from the deep to the shallow layers, indicating a second pulse of upwelling groundwater, TMC 6, located on the top of the transect in an area not covered by trees, shows generally decreasing soil moisture and increasing resistivity in the shallow subsurface that is deepening from about 2 to 3 m over the summer period.This is indicating a reduction in rock moisture, which is accelerating and deepening from the beginning of August, at which time the soil moisture sensors at TMC 5 also show notably decreasing values.This decrease in soil moisture is associated with an increase in resistivity that first is bound to the upper 1 m, and later extends to 2 m depth.Note that for both locations TMC 5 and TMC 6, depths of more than 4 m show a decrease in resistivity, which can be interpreted as a recharge of rock moisture.While for TMC 6 no correlation with rainfall events can be observed, negative changes at deeper depths of TMC 5 seem to relate with strong rainfall events, thereby suggesting water bypassing the shallow soil layers and infiltrating through preferential flowpaths (i.e., fractures) into these deeper layers.This correlation between rainfall and increasing soil moisture also suggests that at forest covered locations rainfall bypasses the canopy, thus reaches the ground, where litter and shading enable some water input into the soil.Hence, in this part of the transect, summer monsoon may actu ally provide some mitigation of drought conditions, which contrasts with the eastern part of the transect, where summer monsoon does not seem to elevate declining soil moisture.We will discuss this further in the following by integrating data and hydrological modeling.

Modeling Subsurface Flow Dynamics
Based on borehole logs (NMR and borehole ERT) at various locations throughout the East River Watershed, a petrophysical transfer function was developed to link the subsurface resistivity to hydraulic permeability and porosity (Uhlemann et al., 2022).Using those relationships, we transformed the inverted resistivity model into subsurface distributions of hydraulic permeability and porosity (Figures 9a and 9b, respectively).We acknowledge that this transformation is likely overestimating both parameters at depth due to the smoothness constraints applied to the resistivity data during the inversion and a lack of resolution (Day-Lewis et al., 2005), and that site specific variations due to changes in lithology may not be thoroughly addressed by this transfer function.Generally, our data acquired across the watershed suggests that bedrock with high fracture density is related to lower resistivity, which is likely due to enhanced weathering in those fractures (Gu et al., 2020;Stolze et al., 2023).Hence, we relate lower resistivity with higher permeability and porosity.Using the recorded precipitation data and derived PET as forcing, we ran the PFLOTRAN model for 10 years, with the first 5 years sufficient for spin-up.We performed particle tracking based on the estimated velocity field and show the flow paths and their average annual shape and length in Figure 9c.The results show that flow at the western part of the transect is mostly lateral and close to the surface, and comparably rapid (shown by long flow paths).Considerable vertical flow seems to occur only at or close to tree locations, indicating that fractures may provide both a place at which trees can relatively easily access rock moisture, but that those fractures may also provide a preferential flowpath for recharge of rock moisture storage.Similar processes have been observed in other critical zone observatories (García-Gamero et al., 2021;Guo et al., 2019).We note that fractures were not modeled as distinct features, but that bedrock areas of higher fracture density are known to have higher porosity and permeability in this area (Uhlemann et al., 2022).Interestingly, rapid and deep lateral flow is observed at the interface between the weathered bedrock and the colluvium, which may provide a recharge pathway to the deeper groundwater.This may also cause the observed increase in resistivity at about x = 120 m, due to a diversion of uphill water into those deeper layers, and hence missing surface water at this location.From x = 150 m, flow patterns are predominantly vertical, highlighting that snowmelt is likely contribution to groundwater recharge.Toward the end of the profile, flow paths turn mostly lateral, indicating drainage into the stream.
We also investigated the flow paths of simpler models assuming a shallow soil layer and no bedrock variation (Figure 10a) and assuming a shallow layer and a bedrock variation as can be derived from the geological map and large scale airborne electromagnetic data (Figure 10b).The results show that a single high porosity and high permeability layer would cause mostly vertical flow across the entire profile, therefore ruling out a mostly topographic control on the flow paths.Assuming a shallow porous and permeable layer overlying a change in bedrock type with a more porous and permeable layer in the gentle part of the slope results in flow pattern that are comparable to the geophysical parameterization of the hydrological model (Figure 9).Shallow, lateral flow is dominating the upper, western part of the transect, while mostly vertical flow dominates the lower, eastern part.Rapid and deep lateral to vertical flow can also be observed at the bedrock boundary.Yet, critical observations, such as deeper flow at the location of trees and lateral flow toward the stream are missing.Hence, the geophysical parameterization allows for a more detailed understanding of the flow processes and how bedrock and vegetation variability are impacting upon the subsurface flow dynamics.Similar conclusions have recently been drawn from a seismic refraction study that used the geophysics derived subsurface property distribution to inform and calibrate a hydrological model of a mountainous headwater catchment (Chen et al., 2023).Note that this model is not accounting for overland flow.Hence, the model may underestimate the infiltration in the lower, eastern part of the transect, where water would accumulate, and may overestimate infiltration in the western part, where some of the precipitation will be lost due to overland flow.
The hydrological model only accounts for the local snowmelt, and recorded rain and evapotranspiration pattern, and shows a single peak in groundwater level followed by slowly declining levels.Hence, the local snowmelt and precipitation patterns, and the associated hillslope-scale subsurface flow cannot explain the double peak observed at SNI1B.Other studies have shown that double peaks are controlled by storage and that different landscape compartments can contribute to the second peak (Martínez-Carreras et al., 2016).Here, we hypothesize that the second peak is caused by snowmelt uphill of the study site.Assessing the snowmelt pattern at and uphill from the study site shows that snowmelt starts 24-41 days and ends 9 to 28 later for mid-mountain and summit locations, respectively, compared to the timing at the site (Figure 11).The second groundwater peak occurs 56 days after the first peak, and 30 days after snowmelt finished on the summit of Mt.Snodgrass (56 days after the end of snowmelt at the mid-mountain site).Assuming that the mid-mountain site controls the time of the second peak (since it is the closest location), this would require a flow velocity of 1.2 ⋅ 10 −4 to 3.2 ⋅ 10 −4 m/s, which falls within the range of observed hydraulic conductivities at the site.The slight but notable decrease in groundwater temperature with the arrival of the second peak further supports snowmelt as the origin of this peak.If this second pulse would be mostly storage controlled, we would expect to not see a significant change in the groundwater temperature (due to longer travel times and mixing with old water and thus attenuation of the temperature signal), and hence we suspect that this second peak is mostly fed by rapid, shallow lateral flow, with limited mixing, as observed in the upper part of the monitoring transect.

Conclusions
Hydrological processes in mountainous watersheds, and how soil, bedrock, and plants interact are still poorly understood.In this study we combine geophysical and environmental monitoring with groundwater flow modeling to image and better understand how variability in bedrock and vegetation properties modulate groundwater flow at a snow-dominated mountainous hillslope.We observe several important hydrological processes that show how variations in bedrock, topography, and vegetation type change subsurface flow pattern, allowing us to answer how subsurface flow pathways differ between shallow and deep bedrock units, and to assess the linkages between vegetation and subsurface flow dynamics across the various bedrock types.
We observe mostly shallow lateral flow on the steep part, which is characterized by thin soil and shallow bedrock.On the gentle slope, which is underlain by colluvium, vertical flow is prevailing that is expected to contribute to groundwater recharge.During snowmelt, vegetation on the shallow bedrock part seems to provide shading, and thus slows snowmelt.Time-lapse resistivity pattern indicate that fractures and collocated tree roots may provide preferential flow pathways into deeper bedrock units.This is clear when comparing processes between areas covered by trees and areas that have no vegetation, where the latter relates to rapid drying and no pattern of deep moisture recharge, while areas covered by trees show heterogeneous infiltration pattern, with deeper changes linked to tree locations.At the steep, shallow bedrock area subsurface conditions seem to mitigate summer drought conditions.In the gentle part of the profile hydraulic dynamics seem to be driven by snowmelt and upwelling groundwater conditions, implying that vegetation will depend more strongly on water supply through these pulse events.This is expected to have implications for the climate change response and their drought resilience at this site, where conifers on fractured bedrock may evolve differently to the conifers on colluvium, which currently experience an abundance of water, and where the trees may draw water from different subsurface compartments.
These observations show that bedrock and vegetation gradients pose a strong control on hillslope hydrology, creating spatially complex flow patterns.We show that bedrock topography has a strong control on groundwater flow pattern, and that the impact of trees on the hydrological response of a hillslope varies with bedrock type.These results highlight the spatial heterogeneity of hydrological processes in mountainous watersheds, which need to be understood to predict how watersheds respond to disturbances.

Figure 1 .
Figure 1.Overview map of the study site.(a) Location of the study site within the East River watershed, shown as hillshaded digital elevation model, with overlain aerial photograph.(b) Map of the study site and its installations, including electrical resistivity tomography monitoring line and the location of boreholes, weather stations, and soil and snow temperature sensors.The green shading indicates the canopy height estimates from LiDAR measurements(Falco et al., 2020).Coordinates are given in UTM Zone 13N, aerial photograph was obtained from Google Earth.

Figure 2 .
Figure 2. Subsurface characteristics of the study site.Interpreted electrical resistivity distribution showing the quartz monzonite porphyry of the upper, western part of the transect, and the landslide deposits overlying Mancos shale in the lower, eastern part.Shown are also the locations of soil moisture, temperature, and snow temperature sensors, as well as the location of a weather station and the groundwater monitoring borehole, and locations of trees along the Electrical resistivity tomography (ERT) profile.The shown size of the tree indicates their maximum canopy height.Borehole logs are shown on the right, showing fluid electrical conductivity (blue) and temperature (black) (b), nuclear magnetic resonance (NMR) derived hydraulic conductivity (K est ) (c), and NMR derived water content (d).Note that the K est was calibrated using pumping tests from nearby boreholes(Uhlemann et al., 2022).

Figure 3 .
Figure 3. Subsurface temperature characteristics.(a) Variability of mean temperature, showing fitted annual mean ground surface temperatures (T mean of Equation 1, circles) and measured annual mean temperature at 0.8 m depth (squares).(b) Shows the variability of the characteristic depth d as estimated through fitting Equation1.Note that the hillshading is representative of the incoming solar radiation during late afternoon of a summer day.

Figure 4 .
Figure 4. Weather and groundwater monitoring data.(a) Daily and cumulative rainfall, and snow accumulation.(b) Mean daily air temperature (2 m above ground), with the blue shaded area indicating the daily range of minimum and maximum temperature.Shown is also daily potential evapotranspiration calculated using the FAO-56 Penman-Monteith method (Allen et al., 1998).(c) Groundwater data, showing the water depth (blue), water temperature (orange), and water electrical resistivity (black) recorded in the borehole at x = 195 m.

Figure 5 .
Figure 5. Geophysical and environmental monitoring data for 2022 snowmelt period.(a) Changes in resistivity with regards to a measurement on 26 March 2022, just before snowmelt commenced.(b, c) Soil moisture measured at depths of 0.1, 0.3, 0.6, and 0.8 m depth, as well as recorded above (i.e., snow) and below-ground temperatures for a station on the eastern (c) and western (b) part of the transect.Black lines indicate the dates of the resistivity data shown in panel (a).

Figure 6 .
Figure 6.Monitoring data of the 2022 snowmelt event.(a-d) Multi-depth soil moisture and above-and below-ground temperatures for locations on the (a, b) eastern and (c, d) western part of the transect.(e) Mean changes in electrical resistivity of the uppermost 1 m along the entire profile.Note that the first snowmelt at the end of March only caused isolated changes in the eastern part of the transect.Green bars on the right side indicate the location of tree cover.

Figure 7 .
Figure7.Soil moisture and changes in resistivity for TMC2 (eastern part of the transect).Shown is the change in resistivity with depth.This shows a cyclic pattern in both the soil moisture and resistivity data.The first event cycle pattern from 17 April to 24 April show a decreasing resistivity from the surface to deeper layers, while for the second event evapotranspiration seems to reduce moisture content in the shallow subsurface.Note that from 05 to 07 resistivity decreases and soil moisture increases from the deeper toward the shallow layers, indication upwelling groundwater.

Figure 8 .
Figure 8. Geophysical and environmental monitoring data for 2022 summer monsoon period.(a) Baseline resistivity distribution.(b-d) Changes in electrical resistivity showing increasing resistivity (loss of moisture) in the shallow subsurface, and isolated and more continuous decreasing resistivity (increase in soil moisture) for the western and eastern part of the transect, respectively.(e) Daily mean air temperature (blue shading indicating daily minimum and maximum temperature), and daily precipitation.(f-h) Changes in resistivity with depth over the summer period for three locations indicated in panel (a).Shown are also the recorded multi-depth soil moisture values.Dashed lines in panel (f) show the recorded soil moisture content at TMC 1.

Figure 9 .
Figure 9. Flow modeling results.(a, b) Showing the permeability and porosity used to parameterize the model.(c) Simplified ground-model showing particle tracking results; colored units indicate the different lithologies, but are not used for parameterization of the model.The western part of the transect shows rapid, mostly lateral flow, with some vertical flow at tree locations.The eastern part shows more homogeneous downward flow toward the water table.

Figure 10 .
Figure 10.Flow modeling results using simplified ground model, assuming (a) a shallow soil layer above a single geological unit, and (b) a shallow soil layer above two geological units.The left panel shows the model parameterization, and the right panel shows the particle tracking results.

Figure 11 .
Figure 11.(a-d) Above ground temperatures with the black line indicating the snow thickness for four locations, with (a, b) being on the eastern and western part of the monitoring transect, respectively, and (c) being halfway between the transect and the summit (d).(c, d) Are at a distance of 570 and 830 m from the study site, respectively.(e) Shows the water depth and temperature recorded at the borehole at the monitoring transect, as well as the recorded daily precipitation.Indicated are also the start and end of snowmelt at the various sites, based on the data in panels (a-d).