Assessing the seasonal compartmentalization of water fluxes in the soil-plant-atmosphere continuum of a high-elevation mountain grassland

Improving our understanding of snow–groundwater connectivity remains a key challenge in high-elevation mountain environments. This calls for a multidisciplinary and multimethod research framework that integrates different types of field observations, including the collection of water samples from diverse sources for stable isotope analysis. However, in remote alpine areas, the limited frequency of sampling hinders the generation of robust, data-driven insights into ecohydrological processes. Therefore, accurately modelling water movement and stable isotope transport through soil, vegetation, and groundwater recharge is essential for advancing our understanding of the hydrological functioning of high-altitude ecosystems.

In this work, we combine a recently introduced snow isotope model with the HYDRUS-1D model to simulate water fluxes and isotope transport within the soil–plant–atmosphere continuum of a high-elevation mountain grassland located in the Aosta Valley, north-western Italy, over the period November 2017–February 2023. We use this modelling framework to:

1. investigate the seasonal origin of two key water fluxes, namely transpiration and deep drainage (the latter assumed to contribute to groundwater recharge)

2. clarify how seasonal water inputs and root water uptake patterns contribute to ecohydrological separation.

The results demonstrate the effectiveness of the proposed modelling framework in accurately simulating volumetric water content (mean absolute error, MAE ≈ 0.03 cm³ cm⁻³ at 10, 20, and 40 cm), actual evapotranspiration (MAE = 0.62 mm d⁻¹), soil (MAE = 1.6 ‰, 3.6 ‰, 3.7 ‰ at 10, 20, and 40 cm) and xylem (MAE = 1.85 ‰) isotope content at the study site. Based on the model outputs, a separation between the water used by plants and the water contributing to deep drainage is evident during the 2018–2020 growing seasons, when median Seasonal Origin Index (SOI) values of transpiration and deep drainage are positive and negative, respectively. Here, positive and negative SOI values indicate an overrepresentation of summer- and winter-derived water in the considered flux. This separation weakens in 2021, when drainage water exhibits a slightly positive median SOI value. Over intense snowmelt periods, meltwater (winter water) rapidly drains through the lower soil layers, whereas rainfall (summer water), which predominantly occurs after the snowmelt period, remains in the soil longer and sustains plant transpiration. The median SOI of transpiration remains consistently positive during the 2018–2021 growing seasons (SOI = 0.19–0.69) but shifts to a negative value in 2022 (SOI = −0.07), highlighting that winter-derived water becomes overrepresented in transpiration fluxes under snow drought conditions. This finding offers valuable insight into how mountain ecosystems may respond to projected increases in temperature and decreases in solid precipitation.

Overall, this work highlights the hydrological conditions that drive the seasonal compartmentalization of water resources in a high-elevation alpine environment, with potential implications for similar mountainous regions worldwide.

1 Introduction

Ecohydrological processes occurring in the subsurface, encompassing both the vadose zone and the phreatic zone, have been often conceptualized relying on physical intuition (Kirchner et al., 2023). However, advancing our understanding of these processes requires detailed insights into the sources and flow paths of subsurface water, which can be gained through measurement techniques (Bovier et al., 2025, 2026; Kirchner et al., 2023). In this regard, stable water isotopes (¹⁸O, ²H) have been widely used to trace the origin of water taken up by plants, and to distinguish among different geographic-sources that contribute to these fluxes (Ceperley et al., 2024; Orlowski et al., 2023; Sprenger et al., 2016). Although precipitation represents the primary input for plant water uptake, its isotopic composition is highly variable in both space and time. Moreover, only a fraction of precipitation infiltrates the soil and is available to roots (Allen et al., 2022); the remainder either runs off at the surface or drains beyond the root zone, contributing to groundwater recharge. As a result, the isotopic signatures (δ¹⁸O, δ²H) of plant water and groundwater recharge may diverge from that of the original precipitation (Allen et al., 2022), complicating efforts to trace water sources and flow paths (Bovier et al., 2025). Despite these challenges, stable isotopes remain a powerful tool in ecohydrological research (White, 1989). They continue to support investigations into questions such as “From where plants take up water?” (von Freyberg et al., 2020), and “What is the seasonal origin of water used by plants?” (Allen et al., 2019a; Floriancic et al., 2024, 2025).

Previous studies have highlighted the potential for ecohydrological separation between water used by plants and that contributing to groundwater recharge or supplying streamflow (McDonnell, 2014). This concept, known as the Two Water Worlds (TWW) hypothesis, proposes that plants access a portion of soil water that is partially disconnected from the water that recharges aquifers and supplies streams (McDonnell, 2017; Evaristo et al., 2015; McDonnell, 2014; Renée Brooks et al., 2010). In essence, this hypothesis posits the existence of two isotopically distinct water pools that exhibit varying degrees of ecohydrological separation bounded by two conceptual extremes (Radolinski et al., 2021; Sprenger and Allen, 2020):

• Mobile water moves rapidly through the soil profile and contributes to groundwater recharge and streamflow.

• Bound water is retained in the soil matrix and primarily used by plants for transpiration.

These isotopically distinct water pools are visible in the dual-isotope plot due to isotopic fractionation mechanisms (Dubbert et al., 2019; Zhou et al., 2021). Soil and plant water samples frequently plot below the Local Meteoric Water Line (LMWL), whereas groundwater and streamwater typically align with it (Evaristo et al., 2015). These pools could derive from different geographic-sources used by roots that would prefer bound soil water (Dubbert et al., 2019) than mobile water. However, the isotopic composition of extracted water can be affected by the choice of sampling technique and extraction conditions (Allen and Kirchner, 2022; Berry et al., 2018; Chen et al., 2020; Ellsworth and Williams, 2007; Millar et al., 2022; Orlowski et al., 2016a, b, 2018). Despite these uncertainties, there has been a longstanding call to test the TWW hypothesis across diverse climates and vegetation types (McDonnell, 2014). The ecohydrological community has actively pursued this challenge, but the hypothesis remains one of the most debated and controversial topics in current hydrological research (Dubbert et al., 2019).

For instance, Barbeta and Peñuelas (2017) found that groundwater constitutes a relevant source of water for plants, and that groundwater–plant connectivity may be both more extensive and quantitatively greater than previously reported. Penna et al. (2013), in a study conducted in a small forested catchment in the Italian pre-Alps, found that beech trees primarily relied on soil water rather than groundwater to support transpiration during late summer and early fall. Radolinski et al. (2021) investigated ecohydrological separation under variable preferential flow conditions using soil columns with varying macropore structures. Their results suggest that mobile water, moving through preferential flow paths, can remain separated from less mobile water and that such separation is most likely to occur following high-intensity precipitation events. Dubbert et al. (2019), based on results obtained by combining δ¹⁸O and δ²H of precipitation, groundwater, soil and xylem water of Quercus suber and Cistus ladanifer with observations of soil water contents and sap flow, stated that the differences in the isotopic composition between soil and plant water vs groundwater can be fully explained by spatio-temporal dynamics of soil-related hydrological processes. Finkenbiner et al. (2022) analyzed hundreds of model configurations of soil, climate, and mobile/immobile soil-water domain characteristics. They concluded that traditional soil physics alone may be insufficient to reproduce large ecohydrological separation. However, they also noted that previous findings of separation could be influenced by other unmodeled processes such as root water uptake dynamics and the interacting effects of seasonality. Overall, these studies highlight that plant water sources and the degree of ecohydrological separation can vary substantially across sites under different climatic and biological conditions. Nevertheless, additional unmodelled factors, including seasonality, may play a key role in modulating ecohydrological separation.

Regarding seasonality, numerous studies have demonstrated the value of δ¹⁸O and δ²H isotopes in investigating seasonal water dynamics within the soil-plant-atmosphere continuum. Seasonal variations in the isotopic composition of precipitation follow a sinusoidal cycle, with typical summer and winter values representing the upper and lower bounds of the cycle, respectively. These summer- and winter-derived isotopic inputs, characterized by markedly different δ¹⁸O and δ²H values, serve as endmembers to identify the seasonal origin of water (Allen et al., 2019b, a). Indeed, Allen et al. (2019a) introduced the Seasonal Origin Index (SOI) to assess the seasonal origins of water in soils and vegetation. The SOI ranges from −1 for water derived entirely from winter precipitation to +1 for water derived entirely from summer precipitation. The SOI = 0 implies that similar fractions of summer and winter precipitation contribute to a considered water flux. The scientific literature presents different results concerning the seasonal origin of water used by plants. Allen et al. (2019a) findings showed negative SOI values in plant water, indicating a substantial contribution from winter precipitation: a pattern also observed by Goldsmith et al. (2022) and Floriancic et al. (2024). Conversely, Zuecco et al. (2026), in the same catchment investigated by Penna et al. (2013), found that trees use predominantly summer water. These differences can be due to varying landscape, climate and plant species characteristics of different study sites (Allen et al., 2019a; Kirchner et al., 2023). Overall, these findings highlight that also the seasonal origin of plant water is highly variable across ecosystems and remains strongly dependent on site-specific climatic, hydrological, and biological conditions.

A complementary perspective involves examining the seasonal origin of streamflow. Since actual evapotranspiration (AET) and discharge (Q) together close the water balance, the seasonal origin of AET must, to some extent, be complementary to that of Q (Allen et al., 2019b). For example, in Swiss catchments having at least four years of streamwater isotope measurements previously analyzed by Seeger and Weiler (2014), von Freyberg et al. (2018) and Bovier et al. (2025), SOI in streamflow (SOI_Q) was found to be ≈ 0 (Allen et al., 2019b), reflecting that streams are sustained by nearly equal fractions of summer and winter precipitation. Thus, the SOI of AET (SOI_AET) fluxes approximates zero, i.e., AET is also sustained by nearly equal fractions of summer and winter precipitation (Allen et al., 2019b). Despite this somewhat counterintuitive result, variations in the amount and timing of annual precipitation, as well as differences in geology or soil and plant characteristics, may lead to different seasonal contributions to both Q and AET.

Thus, understanding how seasonal water inputs and temporal dynamics of root water uptake shape the degree of ecohydrological separation between mobile and bound water pools remains a key challenge in ecohydrology, particularly in high-elevation environments where data scarcity and inherent complexity remain substantial challenges (Gisolo et al., 2024, 2025; van Tiel et al., 2024). Accordingly, the main aims of this paper are to:

• i.

  investigate the seasonal origin of transpiration and deep drainage fluxes across multiple years with contrasting hydrological conditions ranging from very wet to extremely dry

• ii.

  clarify how contrasting seasonal water inputs, typical of high-elevation environments, interact with root water uptake to assess the degree to which the TWW hypothesis holds. Accordingly, we test the following null hypothesis (H₀): “Winter precipitation (i.e., snowmelt) rapidly transits the soil profile recharging groundwater and streams, while summer precipitation (i.e., rainfall) remains available to sustain transpiration fluxes”.

To this end, meteorological data, soil-related data and field measurements of stable isotopes in precipitation, soil water, plant water, and spring water of a high-elevation mountain grassland have been used. Starting from these data, a modelling framework that integrates a recently introduced snow isotope model and the modified version of HYDRUS-1D (Nasta et al., 2023; Stumpp et al., 2012) has been employed, thus enabling a process-based assessment of seasonal water flux partitioning and providing new insights into alpine grassland ecohydrology by partially disentangling the processes underlying ecohydrological separation.

2 Material and Methods

2.1 Study site and datasets

Dora del Nivolet (DOR) is a 16.99 km² catchment located in the western Italian Alps (Valsavarenche (AO), Aosta Valley) with elevations ranging from 2390 to 3430 m a.s.l. (Fig. 2b). The DOR catchment (outlet coordinates: 45°31^′16.66^′′ N; 7°10^′47.44^′′ E) exhibits, under current climatic conditions, a snow-dominated hydroclimatic regime with snow typically accumulating from November and persisting until May, when snowmelt begins (Gentile et al., 2023; Painter et al., 2023). In contrast, the summer period is marked by the onset of rainfall events. Daily winter precipitation was estimated at the study site from snow height measurements from an SR50AT sonic sensor (Campbell Scientific, Inc.), while summer rainfall was reconstructed using data from nearby meteorological stations within the Valsavarenche valley. From summer 2022 onward, summer rainfall was directly measured at the study site using a CS125 weather sensor (Campbell Scientific, Inc.). Between 2017–2018 and 2021–2022 hydrologic years (each hydrological year spans from 1 October to 30 September of the following year), the alpine grassland within the catchment experienced varying meteorological patterns, reflecting a clear trend (Fig. 1).

Figure 1 Cumulative values of rain, snowmelt, and actual evapotranspiration (AET) measurements (derived from the eddy-covariance station installed at the study-site) from hydrological years 2017–2018 to 2021–2022. An exception is made for the year 2017–2018, where measurements begin on 1 November 2017. Coloured bars represent yearly totals, and lines connect the top of each variable's bars across years to highlight temporal trends. Periods without available AET measurements were filled using linear interpolation prior to annual aggregation. Rain and snowmelt values shown in the figure are derived using the model by Ceperley et al. (2020), as described in Sect. 2.4.

The 2017–2018 hydrological year was exceptionally wet and snowy, followed by two relatively wet years with reduced snowmelt. A notable shift occurred in 2020–2021, when the catchment received decreasing amounts of solid precipitation, resulting in lower snowmelt totals. These conditions intensified in 2021–2022, largely due to a snow drought affecting the Italian Alps during the winter of 2021/22 (Koehler et al., 2022). Over the 2017–2018 to 2021–2022 hydrological years, rainfall remained relatively stable with only minor interannual variability. The same was observed for AET, which also remained relatively consistent, although it exhibited a slight increase during the 2022 snow drought. This pattern is consistent with the “drought paradox” described by Mastrotheodoros et al. (2020).

The bedrock geology of the DOR catchment differs between the two sides of the mainstream. The right bank is primarily composed of orthogneisses, granites, metagranites, metagranophyres, porphyroids, and lamprophyric dikes, whereas the left bank features paragneisses, micaschists, and metaconglomerates. Between 2400 and 2600 m a.s.l., talus deposits dominate the right bank, while shallow Dystrict Cambisol are prevalent on the left bank (D'Amico et al., 2020b, a). Around 2400 m a.s.l., alpine meadow is characterized by saturated Fluvisols and peat substrates, through which the mainstream flows (D'Amico et al., 2020b, a). The study catchment (Fig. 2b) hosts a diverse range of plant species, including Gentiana lutea L., Juniperus communis L., Vaccinium myrtillus L., Salix breviserrata Flod., and Trifolium alpinum L. For isotopic analysis of plant water, lignified twigs of Juniperus communis L., Salix breviserrata Flod, Vaccinium myrtillus L., along with roots of Gentiana lutea L. and Trifolium alpinum L., were collected at monthly intervals from June to October. Despite this botanical diversity, the dominant genus across the catchment is Festuca spp. Plant samples were collected on a hillslope located on the left bank of the Dora del Nivolet River, along with soil samples extracted at depths of 10, 20, and 40 cm for isotopic analysis of soil water (Fig. 2a). Both plant and soil waters were extracted via cryogenic vacuum distillation (CVD), performed at the Laboratory of the Faculty of Science and Technology, Free University of Bozen-Bolzano (Italy). The CVD system follows the setup described by Koeniger et al. (2011) and further detailed in Zuecco et al. (2022) and Amin et al. (2021). Isotopic analyses were conducted in the same laboratory. Soil water isotopic composition was measured using cavity ring-down spectroscopy (model L2130-i, Picarro Inc., California, USA), while plant water isotopes were analyzed via isotope ratio mass spectrometry (IRMS; Delta V Advantage Conflo IV, Thermo Fisher Scientific Inc., Waltham, MA, USA). The IRMS was coupled with a Thermo Scientific Gas Bench II to determine δ¹⁸O values, following the methodology of Zuecco et al. (2022).

Figure 2 (a) Location of the study site within the Dora del Nivolet (DOR) catchment. Soil samples and lignified twigs/roots are collected nearby the eddy-covariance station where also volumetric water content and matric potential (10, 20, 40 cm depth) are measured. Blue arrows indicate the flow direction. (b) Geographical framework of the DOR catchment. (c) Photo (2 July 2019) of the monitored source (SOU). (d) Photo (23 August 2019) of the eddy-covariance station and of the precipitation isotope sampler. In the background it is possible to observe the solar panel that supplies the eddy-covariance station.

An eddy-covariance station (45°31^′8.97^′′ N; 7°10^′17.07^′′ E) was installed in November 2017 on a small, flat plateau (Fig. 2d) with a south–southeast aspect at an elevation of 2555 m a.s.l. (Gisolo et al., 2022). The station is equipped with a 3D sonic anemometer (CSAT3B, Campbell Scientific), and a gas analyzer (LI-7500A, Li-Cor). Close to the station (Fig. 2d), a precipitation isotope sampler has been installed. In the same location, the soil profile is instrumented at 10, 20, and 40 cm depths with sensors for volumetric water content (10HS, Meter, accuracy of ±0.05 cm³ cm⁻³) and matric potential (TEROS 21, Meter, accuracy of ±10 % of reading +2 kPa). To provide site-specific pedological characterization, a soil pit was opened on 20 September 2017 near the Eddy-Covariance station and soil samples were collected for physical characterization. The measured percentages of sand, silt and clay are reported in Table 1.

Table 1 Site-specific sand, silt, and clay percentages measured near the Eddy-Covariance station.

Approximately 400 m downslope from the eddy-covariance station, on the same hillslope, a spring (45°31^′8.27^′′ N; 7°10^′36.21^′′ E) named “Source” (SOU) (Fig. 2c) is monitored. Spring water levels are recorded at 10 min intervals using piezoresistive pressure sensors (DL.OCS/N/RS485, STS Sensors), and water samples for isotopic analysis are collected monthly. These samples are analyzed at the Department of Land, Environment, Agriculture and Forestry, University of Padova (Italy), using a liquid water isotope analyzer based on off-axis integrated cavity output spectroscopy (model DLT-100, Los Gatos Research, California, USA).

2.2 Isotopic fractionation correction

The isotopic composition of the water samples is expressed using the delta notation (δ, in ‰), representing the relative deviation from the Vienna Standard Mean Ocean Water (either V-SMOW or R_standard):

$$\begin{array}{}\text{(1)}& \mathit{\delta }\left(\mathrm{‰}\right)={\displaystyle \frac{R_{\mathrm{sample}}-R_{\mathrm{standard}}}{R_{\mathrm{standard}}}}\mathrm{1000}\end{array}$$

Where R_sample is the isotopic ratio (either ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}$O or ${}^{\mathrm{2}}\mathrm{H}{/}^{\mathrm{1}}$H) in the sample, while R_standard is the corresponding ratio in the reference standard. Positive δ values indicate enrichment in heavy isotopes relative to the standard, while negative values indicate depletion. Isotopic fractionation involves changes in the relative abundance of isotopes, such as ¹⁸O, ¹⁶O, ²H and H, due to differences in the physical behavior of heavy and light isotopes during phase transitions (Scandellari and Penna, 2018). The modified version of HYDRUS-1D (Stumpp et al., 2012) does not account for fractionation processes. Therefore, to ensure comparability between measured and simulated isotopic values, a correction based on the Craig and Gordon (1965) model is applied. This correction accounts for both equilibrium and kinetic isotopic fractionation during the liquid-to-vapor phase transition. Due to evaporation, the isotopic composition of residual liquid water diverges from the Local Meteoric Water Line (LMWL), forming a so-called evaporation line (Benettin et al., 2018). Using monthly mean air temperature, vapor pressure, and precipitation δ¹⁸O as inputs, the monthly evaporation line slopes are computed. These slopes are then used to reproject the isotopic composition of fractionated soil and plant water samples back to the Local Meteoric Water Line (LMWL). Technical details of the Craig and Gordon (1965) model are provided in Benettin et al. (2018), together with the MATLAB implementation adopted in this study. To determine whether a water sample is affected by fractionation, the line-conditioned excess^* (lc-excess^*) is calculated, accounting for uncertainty in the isotopic analysis (Landwehr and Coplen, 2004):

$$\begin{array}{}\text{(2)}& {\displaystyle }{\text{lc-excess}}^{*}={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}-a{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-b}{S_u}}\text{(3)}& {\displaystyle }{S}_u=\sqrt{{\mathrm{SD}}_{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}^{\mathrm{2}}+{\left(a\cdot {\mathrm{SD}}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}\right)}^{\mathrm{2}}}\end{array}$$

where δ²H or δ¹⁸O refer to the isotopic composition of the sample, while a and b are the slope and intercept of the LMWL, obtained by linear regression of precipitation isotope data in dual-isotope space (see Sect. 3.2). SD${}_{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}$ and SD${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ denote the standard deviation associated with the isotopic measurement method. For δ²H and δ¹⁸O measured with the Picarro L2130-i analyzer, SDs are 1.0 ‰ and 0.2 ‰, respectively (Marchina et al., 2020). For measurements conducted using the Isotope Ratio Mass Spectrometer (IRMS) Delta V Advantage (Thermo Fisher Scientific), SDs are 2.5 ‰ and 0.1 ‰, respectively. Samples with negative lc-excess^* values are classified as fractionated. For these samples, the intersection points between the LMWL and the corresponding evaporation line are taken to represent their unfractionated isotopic composition. No correction is applied to samples with positive lc-excess^* values.

2.3 Estimation of equivalent precipitation (P_eq) and its isotopic composition ($C_{P_{\mathrm{eq}}}$)

To estimate the timing, quantity, and isotopic composition of liquid water inputs, we employ the snow accumulation, melt, and isotope model developed by Ceperley et al. (2020). While this model involves simplifications compared to more physically detailed snow models, such as Snowpack (Lehning et al., 1999), Crocus (Brun et al., 1989, 1992; Vionnet et al., 2012), GEOtop (Rigon et al., 2006; Endrizzi et al., 2014), Amundsen (Strasser et al., 2024) and FSM (Essery, 2015), its objective in this study is to simulate the timing and amount of snowmelt so that the snow model could be considered a statistically acceptable proxy for the snow energy balance. In this context, liquid water input refers to rainfall (P_R) and snowmelt (SM) whose combined value is termed equivalent precipitation ($P_{\mathrm{eq}}=P_{\mathrm{R}}+$ SM). While HYDRUS-1D includes a snow routine similar to that of Ceperley et al. (2020), it does not explicitly simulate the isotopic composition of infiltrating water from both snowmelt and rainfall, nor their potential mixing (Stumpp et al., 2012). Therefore, the estimation of P_eq and its associated isotopic composition ($C_{P_{\mathrm{eq}}}$) is conducted externally to HYDRUS-1D using the Ceperley et al. (2020) model. A brief description of this snow accumulation, melt, and isotope model is provided below. For full technical and methodological details, readers are directed to the original publication, which includes a MATLAB implementation (used in this study) made publicly available by the authors.

2.4 Snow accumulation and melt model

The snow accumulation and melt model proposed by Ceperley et al. (2020) computes snow dynamics at the catchment scale using an elevation band approach. In this study, the DOR catchment is discretized into 100 elevation bands. Snow accumulation and melt are calculated at a daily resolution from November 2017 to February 2023. A key input for these calculations is air temperature (T), which is linearly interpolated across elevation bands using a fixed temperature lapse rate of −3.44 °C per 1000 m. This temperature lapse rate was derived from air temperature data of both the monitoring station described in Sect. 2.1 (2555 m a.s.l.), and three additional monitoring stations managed by the Centro Funzionale Valle d'Aosta located in the Valsavarenche valley – Valsavarenche-Orvieille (2170 m a.s.l.), Valsavarenche-Pont (1951 m a.s.l.) and Valsavarenche Eaux-Rousses (1651 m a.s.l.) – thus accounting for local site effects (Eeckman et al., 2022). Precipitation is held spatially constant across the catchment, as in the original implementation by Ceperley et al. (2020). Since equivalent precipitation (P_eq) is computed separately for each elevation band, we used the P_eq time series corresponding to the elevation band containing the full experimental setup – eddy-covariance station, soil sensors, and isotope sampling sites (Fig. 3a) – as input for HYDRUS-1D.

Figure 3 HYDRUS-1D input and related data obtained as described in Sect. 2.3–2.6. (a) Equivalent precipitation (P_eq) and corresponding isotopic composition ($C_{P_{\mathrm{eq}}})$, (b) maximum (T_max), mean (T_mean) and minimum (T_min) daily air temperature, (c) Leaf Area Index (LAI), (d) Potential evapotranspiration (ET_P) computed with HSM, potential evaporation (E_P) and potential transpiration (T_P) obtained by applying the Beer's law. Panels (a) and (d) show data of the time-variable boundary condition inputs used in HYDRUS-1D, while panels (b) and (c) present the intermediate variables used to compute the inputs in panel (d).

Within each elevation band, snow accumulation is modeled using a linear temperature-based transition between liquid and solid precipitation (Harpold et al., 2017; Jarvis, 1994). According to Hock (2003), precipitation is classified as snow (P_S) when temperatures are below a lower threshold (T_S), and as rain (P_R) when temperatures exceed an upper threshold (T_R). Snowmelt begins once air temperature surpasses a defined melting threshold (T₀), and is simulated using a degree-day method (Schaefli et al., 2014):

$$\begin{array}{}\text{(4)}& \mathrm{SM}\left(t\right)=\left\{\begin{array}{ll}\max \left[\mathit{\xi }\left(T\left(t\right)-T_{\mathrm{0}}\right),\mathrm{SWE}\left(t\right)\right],& \mathrm{if}T\left(t\right)>T_{\mathrm{0}}\\ \mathrm{0},& \mathrm{if}T\left(t\right)\le T_{\mathrm{0}}\end{array}\right.\end{array}$$

Where SM(t) is the daily snowmelt (mm d⁻¹) at the time t, SWE(t) is the snow water equivalent (mm) at the time t, and T(t) is the mean daily air temperature (°C) at the time t. ξ (mm °C⁻¹ d⁻¹) is the degree-day factor. The parameters used in the model are summarized in Table 2.

Table 2 Parameters used in the snow accumulation and melt model. The listed parameters were selected based on values reported in the scientific literature (Ceperley et al., 2020; Jarvis, 1994; Schaefli et al., 2014).

Ideally, these parameters should be calibrated using direct snowmelt observations, which were unavailable in this study. Therefore, the default values from HYDRUS-1D's internal snow routine were adopted, which have shown good consistency with values reported in the literature. For instance, the degree-day factors (ξ) typically range from 1.6 to 6 mm °C⁻¹ d⁻¹ (Van Mullem et al., 2004), and Ceperley et al. (2020) used ξ values between 2.7 and 5 mm °C⁻¹ d⁻¹ across three alpine catchments. The HYDRUS-1D default value (ξ=4.3 mm °C⁻¹ d⁻¹) falls within this range and was also used by Stumpp et al. (2012). The temperature thresholds used in snow accumulation models vary across studies. The T_R and T_S values proposed by Jarvis (1994), consistent with the HYDRUS-1D implementation, were adopted. The melting threshold T₀ was set to 0 °C, following Schaefli et al. (2014) and Ceperley et al. (2020).

2.5 Snow isotope model

The isotopic composition of water stored in the snowpack (C_S) is computed for each elevation band based on three key assumptions:

1. Complete water mixing within the snowpack.

2. Any rainfall falling on an existing snowpack mixes with the water stored within it.

3. Rainfall mixing with the snowpack is assumed to exit the system within the same time step (t).

The first assumption is generally valid during peak snowmelt periods, when the snowpack is isothermal, but less accurate during early or intermittent melt events (Ceperley et al., 2020). The third assumption simplifies the process by neglecting the snowpack's water holding capacity and the potential for temporary refreezing (Schaefli et al., 2014).

Under these assumptions, the snowpack isotopic mass balance equation can be defined as follows (Ceperley et al., 2020):

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}(\mathrm{SWE}\left(t\right)\cdot C_{\mathrm{S}}\left(t\right))}{\mathrm{d}t}}=C_{\text{P-fit}}\left(t\right)\cdot P\left(t\right)-C_{\mathrm{S}}\left(t\right)\cdot P_{\mathrm{eq}}\left(t\right)\end{array}$$

Where C_P-fit(t) is the isotopic composition of precipitation (P(t)) at the time t, derived by fitting a sine curve to observed data, and C_S(t) is the isotopic composition of the water stored in the snowpack at the time t. During time steps in which SWE$\left(t\right)+P_{\mathrm{eq}}\left(t\right)>\mathrm{0}$, the input isotopic composition is C_S(t), which is computed by numerically solving Eq. (5) through a time-stepping approach, initialized with: C_S(t₀)=C_P-fit(t₀). Else, the input isotopic composition is C_P-fit. Since the isotopic composition of equivalent precipitation ($C_{P_{\mathrm{eq}}}$) is computed for each elevation band, we used the $C_{P_{\mathrm{eq}}}$ time series corresponding to the elevation band containing the full experimental setup – eddy-covariance station, soil sensors, and isotope sampling sites (Fig. 3a) – as input for HYDRUS-1D.

2.6 Potential evapotranspiration (ET_P), evaporation (E_P) and transpiration (T_P)

To compute the potential evapotranspiration (ET_P), the approach proposed by Ravazzani et al. (2012), who calibrated a modified version of the Hargreaves–Samani (HS) equation using potential evapotranspiration data derived from the FAO-56 Penman–Monteith method, was adopted. Their modification introduces a correction factor that incorporates two calibration parameters along with the elevation of the meteorological station (Elev_S). The calibration was performed for the Upper Po and Rhone River basins, thereby including our study area within the Alpine domain considered by Ravazzani et al. (2012). The resulting modified Hargreaves–Samani equation (HSM) is used to compute daily ET_P as follows:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathrm{ET}}_{\mathrm{P}}\left(t\right)& =(\mathrm{0.817}+\mathrm{0.00022}\cdot {\mathrm{Elev}}_{\mathrm{S}})\cdot \mathrm{HC}\cdot R_{\mathrm{a}}\cdot \left(T_{max}\left(t\right)\right.\\ & {\left.-T_{min}\left(t\right)\right)}^{\mathrm{HE}}\cdot \left({\displaystyle \frac{T_{max}\left(t\right)+T_{min}\left(t\right)}{\mathrm{2}}}+\mathrm{HT}\right)\end{array}\end{array}$$

Where ET_P(t) (mm d⁻¹) is the daily ET_P computed with HSM (Fig. 3d), Elev_S (m a.s.l.) is the station elevation (2555 m a.s.l. in our study), HC is an empirical coefficient equal to 0.0023, R_a (mm d⁻¹) is the extraterrestrial radiation, HE is an empirical exponent equal to 0.5, HT is a factor used to convert units from Fahrenheit to Celsius and equal to 17.8, T_max (°C) is the daily maximum air temperature (Fig. 3b) and T_min (°C) is the daily minimum air temperature (Fig. 3b).

ET_P is then partitioned into potential evaporation (E_P) and potential transpiration (T_P) using Beer's law that partitions the solar radiation component of the energy budget via interception by the canopy (Ritchie, 1972) as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{E}_{\mathrm{P}}\left(t\right)={\mathrm{ET}}_{\mathrm{P}}\left(t\right)e^{-k\cdot \mathrm{LAI}\left(t\right)}\text{(8)}& {\displaystyle }{T}_{\mathrm{P}}\left(t\right)={\mathrm{ET}}_{\mathrm{P}}\left(t\right)-E_{\mathrm{P}}\left(t\right)\end{array}$$

Here, k=0.463 is the extinction coefficient for global solar radiation within the canopy (Ritchie, 1972), and LAI is the Leaf Area Index (Fig. 3c). The LAI time series (November 2017 to February 2023) was derived using a Google Earth Engine script applied to the MODIS/061/MCD15A3H image collection (spatial resolution: 500 m; temporal resolution: 4 d) for the pixel encompassing the eddy-covariance station. To obtain daily values, linear interpolation was performed between 4 d intervals. The resulting daily time series of E_P and T_P (Fig. 3d) are used as input to HYDRUS-1D. No further partitioning of potential evaporation into soil evaporation and snow sublimation is made: this partitioning is beyond the scope of the present study. During winter, this term predominantly reflects sublimation due to the presence of snow cover at the study site.

2.7 HYDRUS-1D: main equations and model set up

The methodological sections above describe how the input data (Fig. 3) were obtained for simulating water flow and isotope transport within the soil profile using HYDRUS-1D (Stumpp et al., 2012). Prior to importing these inputs into HYDRUS-1D, all fluxes were converted from mm d⁻¹ to cm d⁻¹ and a positive constant was added to the isotopic compositions of equivalent precipitation to ensure all values were positive to run HYDRUS-1D (Stumpp et al., 2012). The same constant will then be subtracted during post-processing to restore the isotopic composition in ‰. HYDRUS-1D (Šimůnek et al., 2018) simulates variably saturated water flow by solving Richards' equation:

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}=-{\displaystyle \frac{\partial q}{\partial z}}={\displaystyle \frac{\partial }{\partial z}}\left[K\left(h\right){\displaystyle \frac{\partial h}{\partial z}}+K\left(h\right)\right]-S\end{array}$$

Where θ is the volumetric water content (cm³ cm⁻³), z is the depth below soil surface (positive upward, in cm), q is the water flux (cm d⁻¹), t is time (d), S is a sink term representing root water uptake (d⁻¹), and K(h) is the soil hydraulic conductivity (cm d⁻¹), which is a function of pressure head h (cm) and θ. The soil water retention curve θ(h) and the hydraulic conductivity function K(h) are described by the van Genuchten (1980) and Mualem (1976) models, respectively:

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathit{\theta }\left(h\right)=\left\{\begin{array}{ll}{\mathit{\theta }}_{\mathrm{r}}+\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}{(\mathrm{1}+{\left|{\left.\mathit{\alpha }h\right|}^n)\right.}^m}& h<\mathrm{0}\\ {\mathit{\theta }}_{\mathrm{s}}& h\ge \mathrm{0}\end{array}\right.\text{(11)}& {\displaystyle }K\left(h\right)=K_{\mathrm{s}}{\displaystyle \frac{{\left[\mathrm{1}-{\left(\mathit{\alpha }h\right)}^{n-\mathrm{1}}{\left[\mathrm{1}+{\left(\mathit{\alpha }h\right)}^n\right]}^{-m}\right]}^{\mathrm{2}}}{{\left[\mathrm{1}+{\left(\mathit{\alpha }h\right)}^n\right]}^{\frac{m}{\mathrm{2}}}}}\end{array}$$

Where θ_r and θ_s (cm³ cm⁻³) are the residual and saturated water contents, respectively; α (cm⁻¹), n and m ($=\mathrm{1}-n^{-\mathrm{1}}$) are empirical shape parameters; and K_s is the saturated hydraulic conductivity (cm d⁻¹). Root water uptake is simulated using the stress response function of Feddes and Zaradny (1978) with parameters for grass selected from the HYDRUS-1D internal database. Root depth was observed to extend down to 60 cm (equal to the simulated soil profile depth), and roots distribution is assumed to be homogeneous throughout the profile (Stumpp et al., 2012). Isotope transport is modeled using the advection–dispersion equation:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial \left(\mathit{\theta }C\right)}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\theta }D{\displaystyle \frac{\partial C}{\partial z}}\right)-{\displaystyle \frac{\partial \left(qC\right)}{\partial z}}-SC\end{array}$$

Where C is the tracer concentration (‰), D (cm² d⁻¹) is the dispersion coefficient and S (d⁻¹) is the root water uptake. The dispersion coefficient D is calculated based on Bear (1972) for one-dimensional transport:

$$\begin{array}{}\text{(13)}& D={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{L}}q}{\mathit{\theta }}}+D_{\mathrm{w}}{\mathit{\tau }}_{\mathrm{w}}\end{array}$$

Where λ_L is the longitudinal dispersivity (cm), D_w is the molecular diffusion coefficient in free water (10⁻⁹ m² s⁻¹) and τ_w is the tortuosity factor of Millington and Quirk (1961).

Two assumptions are made in the isotope transport modeling:

• evaporation does not fractionate water, meaning water and isotopes exit at the same rate at the upper boundary (Stumpp et al., 2012)

• root water uptake does not induce isotopic fractionation.

Model parameters for water flow and isotope transport are calibrated using the inverse modeling tool embedded in HYDRUS-1D, which employs the Marquardt–Levenberg optimization algorithm. Observations used for calibration include volumetric water content, pressure head, and isotopic composition at depths of 10, 20, and 40 cm, whereas model validation relied on observed AET and xylem water isotopic composition.

At the soil profile upper boundary, we use the “atmospheric boundary condition (BC) with Surface Runoff” condition thus leading both external and soil conditions to control the water flux across the upper boundary. Moreover, we set up a “concentration flux BC”, thus specifying liquid phase concentration of the infiltrating water. At the soil profile lower boundary, we set up “free drainage (zero gradient)” BC that well describes water flow (solute transport) in the vadose zone field studies. A schematic overview of the HYDRUS-1D model setup is provided in Fig. 4.

Figure 4 HYDRUS-1D setup. Red squares indicate observation points in which isotope and soil probes measurements are available. Blue square indicates an observation point in which measurements are not available.

2.8 The Seasonal Origin Index (SOI)

To assess whether a seasonal separation exists between the water used by plants and the water contributing to groundwater recharge and streamflow, the Seasonal Origin Index (SOI), following Allen et al. (2019a), was calculated:

$$\begin{array}{}\text{(14)}& \mathrm{SOI}=\left\{\begin{array}{ll}\frac{{\mathit{\delta }}_x-{\mathit{\delta }}_{\mathrm{annP}}}{{\mathit{\delta }}_{\mathrm{summerP}}-{\mathit{\delta }}_{\mathrm{annP}}},& \mathrm{if}{\mathit{\delta }}_x>{\mathit{\delta }}_{\mathrm{annP}}\\ \frac{{\mathit{\delta }}_x-{\mathit{\delta }}_{\mathrm{annP}}}{{\mathit{\delta }}_{\mathrm{annP}}-{\mathit{\delta }}_{\mathrm{winterP}}},& \mathrm{if}{\mathit{\delta }}_x<{\mathit{\delta }}_{\mathrm{annP}}\end{array}\right.\end{array}$$

where δ_x denotes the fractionation-compensated isotopic composition of the considered flux, while δ_winterP, δ_summerP, and δ_annP correspond to the isotopic compositions of typical winter, typical summer, and volume-weighted annual precipitation, respectively. The values of δ_winterP and δ_summerP are defined as the minimum (−18.12 ‰) and maximum (−7.85 ‰) of the sinusoidal fit (C_P-fit) describing the seasonal precipitation isotope cycle, whereas δ_annP (−13.28 ‰) is derived directly from the observational dataset. The SOI ranges from −1 for water derived entirely from winter precipitation to +1 for water derived entirely from summer precipitation. When considering a given water flux (e.g., transpiration, evaporation, or drainage), SOI = 0 implies that similar fractions of summer and winter precipitation contribute to that flux. In this study, the SOI is calculated using the isotopic compositions (δ_x) of water fluxes simulated with HYDRUS-1D. The SOI is adopted because, as noted by Allen et al. (2019a), it is specifically designed to assess whether winter or summer precipitation is overrepresented in the considered flux. In other words, the SOI accounts for site-specific seasonality in precipitation, recognizing that, at the study site, we should expect a greater proportion of water fluxes deriving from winter inputs since precipitation is unevenly distributed over the year and predominantly occurs during the winter season. Moreover, this metric has proven to be relatively insensitive to several sources of uncertainty commonly affecting isotope-based rooting depth analyses, particularly those related to sampling and extraction of soil water that accurately represents the water taken up by roots (Allen et al., 2019a; Goldsmith et al., 2019; Orlowski et al., 2018; Penna et al., 2018).

3 Results and Discussion

3.1 Soil hydraulic and solute transport parameters

The optimized soil hydraulic and solute transport parameters are summarized in Table 3.

Table 3 Optimized soil hydraulic and transport parameters.

The physical plausibility of the optimized parameters is supported by both site-specific measurements and values reported in the literature. The optimized θ_S resulted higher than the maximum volumetric water content (0.43 cm³ cm⁻³) measured by soil probes and aligns to the upper limit of typical porosity values for loamy sand textured soils. Indeed, according to Clapp and Hornberger (1978), a porosity of approximately 0.410±0.068 cm³ cm⁻³ is a representative value for this soil texture and defines the upper limit of volumetric water content (Nimmo, 2013). The optimized θ_r is consistent with both the wilting point value at 1500 kPa (0.142 cm³ cm⁻³) and the average measured water content for pressure heads exceeding 10⁵ cm (0.14 cm³ cm⁻³). A non-zero residual water content is further supported by field observations, which showed that even during the extreme drought of 2022, the lowest recorded volumetric water content remained around 0.10 cm³ cm⁻³. By considering the solute travel distance, i.e., the 60 cm soil profile depth, the optimized longitudinal dispersivity λ_L falls within the range of 0.9–20 cm indicated by Vanderborght and Vereecken (2007) for the 31–80 cm travel distance class. The comparison between measured and optimized soil water retention functions is illustrated in Fig. 5.

Figure 5 Volumetric water content (θ) as a function of pressure head ($\left|h\right|$). Points represent measurements, while the solid red line denotes the optimized van Genuchten water retention function. Performance metrics (RMSE, MAE, and Pearson's correlation coefficient, r) are reported.

3.2 Isotopic fractionation correction

Using the Craig and Gordon (1965) model as implemented by Benettin et al. (2018), monthly evaporation slopes (ES) were derived to identify the original isotopic signatures of soil and plant water that had undergone evaporation-driven fractionation. The resulting mean evaporation slope (MES) was 3.39±0.05, with monthly ES values varying modestly from 3.33 in July to 3.49 in May (Fig. 6a). The isotopic composition ranges of soil and plant water samples, before and after fractionation correction, are presented in Table 4.

Table 4 Range of isotopic composition of soil and plant water pre- and post- the Craig and Gordon (1965) model application.

Figure 6 (a) Results of the monthly evaporation slopes (ES) obtained by applying the Craig and Gordon (1965) model being implemented by Benettin et al. (2018). The monthly ES have been used to reproject the fractionated water samples on the local meteoric water line (LMWL). The numbers next to the residual liquid points indicate the month (1: January; 12: December); (b) Isotopic fractionation correction for fractionated (line-conditioned excess${}^{*}<\mathrm{0}$) plant/soil water samples. Precipitation and Source (SOU) water samples have been also reported. The median isotopic composition of each reservoir is indicated by the dashed lines.

When compared to the median isotopic composition of precipitation (δ¹⁸O: −11.5 ‰, δ²H: −84.6 ‰), the median isotopic composition of plant water (δ¹⁸O: −10.5 ‰, δ²H: −70 ‰) clearly shows enrichment in heavy isotopes, a pattern typical of summer precipitation (Fig. 6b). Conversely, the median isotopic composition of spring water at SOU (δ¹⁸O: −14.5 ‰, δ²H: −105.8 ‰) indicates marked depletion, characteristic of winter precipitation (Fig. 6b). The δ¹⁸O and δ²H values of SOU water are also confined within a narrow range (δ¹⁸O: −15.6 ‰ to −12.8 ‰; δ²H: −116.6 ‰ to −95.8 ‰). These observations provide a first line of evidence for a seasonal partitioning between water sources used by vegetation and those contributing to streamflow. Specifically, the dual-isotope plot suggests that transpiration fluxes are predominantly supported by summer precipitation, whereas discharge at the SOU spring is mainly sustained by winter precipitation inputs. Interestingly, the median isotopic composition of soil water (δ¹⁸O: −11.9 ‰; δ²H: −84.1 ‰) closely matches that of precipitation, and both span similar isotopic ranges. In this regard, Radolinski et al. (2021) showed that δ¹⁸O of soil water is more sensitive to changes in precipitation signature than drainage water. A physical explanation of these empirical observations will be provided in Sect. 3.5, where the HYDRUS-1D results are presented and discussed.

3.3 Model performance and Compartment-based calibration/validation approach

Model performance metrics for all evaluated variables are reported in Table 5. The adopted performance metrics are Pearson's linear correlation coefficient (r), the mean absolute error (MAE) and the Root Mean Square Error (RMSE).

Table 5 Model performance metrics for all evaluated variables. The adopted performance metrics are Pearson's linear correlation coefficient (r), the mean absolute error (MAE) and the Root Mean Square Error (RMSE).

These HYDRUS-1D performance are comparable to those achieved in other mountainous environments, such as those reported by Bertoldi et al. (2014) using the GEOtop model and Gisolo et al. (2024) using HYDRUS-1D configured for double vegetation. The estimated δ¹⁸O satisfactorily describe the dynamics of measured soil and plant water, but with lower accuracy at 20 and 40 cm depths. A good correspondence was found between simulated and observed AET: performance metrics are consistent with those reported by Gisolo et al. (2024) who simulated AET dynamics in an abandoned alpine grassland, also accounting for shrub encroachment, using HYDRUS-1D configured with a double vegetation.

The calibration/validation strategy adopted in this study was designed to support a process-based reconstruction of ecohydrological dynamics over the 2017–2023 period. While predictive applications typically require temporally independent out-of-sample validation, the objective here is to reproduce and interpret the processes that occurred during the monitored years. Owing to the intrinsic complexity of high-elevation mountain environments, the temporal resolution of isotope sampling is necessarily limited, which allows only preliminary insights into the ecohydrological processes occurring within the system. For instance, seasonal compartmentalization of water fluxes was already suggested by the dual-isotope analysis (Fig. 6); however, a quantitative assessment of the processes leading to this compartmentalization requires a model capable of consistently reproducing the occurring dynamics both temporally (at daily resolution) and spatially (along the soil profile) throughout the entire study period.

To this end, a compartment-based validation strategy was adopted. As reported in Sect. 2.7, the model was calibrated using observations from the soil compartment and subsequently evaluated against observations from the soil-plant compartment related to outgoing fluxes toward the atmosphere, specifically AET measurements (derived from the eddy-covariance technique) and xylem water isotopic composition. The model reproduces actual evapotranspiration and the xylem water isotopic composition satisfactorily (Table 5), despite not being calibrated directly against variables related to outgoing fluxes toward the atmosphere. This consistency reduces the likelihood that good performance is achieved for the wrong reasons and supports a good description of the processes occurring in the soil-plant-atmosphere continuum under study. Furthermore, in Sect. 3.1 the calibrated parameters were evaluated against ranges reported in the scientific literature to ensure physical plausibility and, where possible (e.g., for θ_r and θ_S), estimated parameter values were also directly compared with field measurements, providing additional validation and strengthening confidence in the robustness of the parameterization.

The analyzed hydrological years exhibit marked interannual variability, including substantial differences in snowmelt dynamics (Fig. 1). Accordingly, calibrating the model on a limited temporal subset and validating it on a subsequent period could bias parameter estimation toward the calibration years and reduce robustness under contrasting hydroclimatic conditions. Thus, a traditional split-sample approach might have introduced differences in model reliability between calibration and validation periods, potentially leading to lower confidence in process-based interpretations during the validation phase relative to the calibration phase. By calibrating the model across multiple years characterized by highly variable conditions, the resulting parameterization aims to represent a very broad range of system behaviour.

3.4 Simulated volumetric water content, soil/plant water isotopic composition and actual evapotranspiration

Simulated volumetric water content at the three observation depths (10, 20, and 40 cm) is shown in Fig. 7. A strong agreement is observed between simulated and measured values (Table 5). However, some discrepancies between model output and observations are noted, particularly during winter. The peaks in volumetric water content simulated by HYDRUS-1D during winter correspond to simulated snowmelt events. In snow-dominated catchments like DOR, significant winter snowmelt events contributing to runoff are rare compared to hybrid catchments (Gentile et al., 2023, 2024). Still, minor snowmelt can occur. These events, identified by the degree-day model, coincide with observed declines in snowpack depth (Fig. 7d), but they do not translate into measurable increases in volumetric water content. This is likely due to refreezing processes within the snowpack, which are not represented in the model and which inhibit water infiltration during winter (Lundberg et al., 2016; Hirashima et al., 2017; Leroux and Pomeroy, 2017). Accurately capturing the influence of refreezing would require explicitly modelling its effects on both the isotopic composition of snowpack and meltwater over successive melt–freeze cycles: a complexity beyond the scope of this study. Zhou et al. (2008) revealed that the refreezing process would inevitably result in a refrozen snowpack characterized by a line on the dual-isotope plot with a decreased slope compared to solid phase of the initial melting snowpack. Consequently, the line representing the refrozen snowpack on the dual-isotope plot would diminish progressively its slope with each diurnal melt-freeze cycle. On the other side, the line representing the liquid phase on the dual-isotope plot shows an overall slight decrease in the melting period (Zhou et al., 2008).

Figure 7 Measured (line with markers) and simulated (line) volumetric water content at (a) 10 cm, (b) 20 cm and (c) 40 cm. Snowmelt and rain have been also indicated. (d) Snow depth, snowmelt and rain. θ indicates the volumetric water content.

Further uncertainty arises from estimating the amount and timing of snowmelt at each timestep (Stumpp et al., 2012). Indeed, this is linked to the uncertainty of the parameters used in the degree-day model of Ceperley et al. (2020). The degree-day model relies on air temperature to trigger snowmelt, but in some conditions, it could poorly include the effect of other factors such as topographic shading and proximity to snow-free areas which are identified as further drivers of snowmelt by a stochastic cellular automaton model applied at this site (Painter et al., 2023). In this regard, Bertoldi et al. (2010) highlighted the role of topography, variable precipitation, and solar radiation in shaping volumetric water content patterns, which ultimately affect AET. As in our study, also Stumpp et al. (2012), by assessing the effects of land cover and fertilization on water flow and solute transport of five lysimeters using HYDRUS-1D, found the main discrepancies between simulated and measured values because of the uncertainties related to infiltration during snowmelt. Therefore, for future studies, a more accurate estimation of these fluxes could be achieved by directly using instrumental methods, as reported in Eeckman et al. (2025).

Overall, the degree-day approach used in this work has the limitation of providing a simplified representation of snowmelt processes, in which processes such as compaction, refreezing, stratification, and snow sublimation are not explicitly included. Regarding the latter, a further partitioning of potential evaporation into soil evaporation and snow sublimation is beyond the scope of the present study. The main interest here is to adequately represent the outgoing flux toward the atmosphere. Moreover, explicitly accounting for the isotopic composition of sublimation would require additional model developments, as this process remains subject to ongoing research (Beria et al., 2018). Therefore, no further distinction is made between soil evaporation and snow sublimation, as both processes are included within the potential evaporation term computed using Eq. (7). During winter, this term predominantly reflects sublimation due to the presence of snow cover. Accordingly, we compared the simulated AET with eddy-covariance measurements, which do not distinguish between sublimation, evaporation, and transpiration.

Notably, at the end of snowmelt periods (June), measured volumetric water content is often slightly higher than simulated values. This discrepancy may be explained by the impact of rain-on-snow events. Indeed, most of these conditions (i.e., precipitation falling as rain while the SWE is greater than 0) occur in June (Table S1 in the Supplement, Fig. S1). During these events, the Ceperley et al. (2020) model simulates melt fluxes (Fig. S1). However, the magnitude of snowmelt is likely underestimated, as rain-on-snow events are known to generate more intense and short-lived melt pulses compared to melt driven solely by temperature (Myers et al., 2023). The lack of an explicit representation of rain-on-snow processes therefore constitutes an additional limitation of the modeling approach. Nevertheless, it is worth noting that such conditions account for only 5.6 % of the total study period, mainly occurring at the end of the melting season.

Despite these limitations, our results indicate that the degree-day approach provides overall reliable estimates of snowmelt dynamics at the Alpine study site where the melt is largely driven by latent heat transfers (Ceperley et al., 2020; Ohmura, 2001).

Simulated isotopic composition of plant water and soil water at the three observation depths (10, 20, and 40 cm) is shown in Fig. 8b. It is important to consider the uncertainties in isotopic measurements introduced by water extraction techniques. Millar et al. (2022) reviewed the accuracy (expressed as standard deviation, SD) of various extraction methods. For cryogenic vacuum distillation (CVD), the SD of δ¹⁸O can range from 0.09 ‰ to 2.3 ‰. To reflect this methodological uncertainty, the maximum SD values reported by Millar et al. (2022) are displayed as error bars for each measured δ¹⁸O value in Fig. 8b.

Figure 8 (a) Snowmelt, rain and isotopic composition of equivalent precipitation ($C_{P_{\mathrm{eq}}})$. (b) Comparison of simulated soil/plant water isotopic composition at the three observation nodes (10, 20 and 40 cm) and measured soil/plant water isotopic composition. The maximum standard deviation (SD) values reported by Millar et al. (2022) are displayed as error bars.

The Ceperley et al. (2020) snow model used in this study includes simplifying assumptions, such as neglecting the snowpack's water holding capacity and temporary refreezing. These assumptions may influence the isotopic composition of equivalent precipitation ($C_{P_{\mathrm{eq}}}$) which in turn affects the simulated δ¹⁸O values of soil and plant water which may therefore show discrepancies compared to the observed δ¹⁸O values. The agreement between simulated and measured isotopic composition declines during 2022, particularly at 20 and 40 cm depths. This may be attributed to the assumption of complete mixing within the snowpack, which likely does not hold under conditions of a more ephemeral snowpack, such as those observed in 2022.

The scarcity of wintertime field data in high-elevation environments makes the simulated δ¹⁸O of soil water particularly valuable, as it provides insight into soil hydrological processes that cannot be directly observed during this season. Following the growing season, soil water shows isotopic compositions closer to those of summer precipitation, remaining relatively stable in the absence of early-season snowmelt events that could otherwise modify soil isotopic dynamics during winter. Interestingly, during winter, δ¹⁸O values are less depleted at 40 cm and become progressively more depleted toward the surface. This pattern is reversed during the growing season, with more depleted signatures at deeper soil layers. Such behavior supports the hypothesis that isotopically depleted snowmelt filtrates vertically and contributes to deep drainage during the growing season (Gentile et al., 2023; Cochand et al., 2019; Du et al., 2019; Flerchinger et al., 1992).

To further evaluate the performance of the calibrated HYDRUS-1D parameters for simulating water flow and isotope transport, we compared simulated actual evapotranspiration (AET H-1D) against AET derived from eddy-covariance (AET Eddy) measurements (Fig. 9a, b).

Figure 9 (a) Measured (AET Eddy) and simulated (AET H-1D) actual evapotranspiration. (b) AET H-1D versus AET Eddy.

3.5 Variable degrees of ecohydrological separation driven by time-variable seasonal water inputs

The comparison between model outputs and observations was used to evaluate the ability of HYDRUS-1D to reproduce volumetric water content, soil and plant water isotopic composition, and AET in the high-elevation grassland under study. As shown in Sect. 3.4 the model performs well and produces results consistent with those reported in similar alpine contexts (Bertoldi et al., 2014; Gisolo et al., 2024). Building on this validation, we next analyze model outputs to explore the hydrological processes occurring within the soil-plant-atmosphere continuum, with a specific focus on the seasonal partitioning of winter and summer precipitation between plant water uptake and deep drainage (assumed to recharge groundwater). To gain a deeper understanding of this topic, we calculate the SOI (as described in Sect. 2.8) starting from the simulated isotopic composition of the water fluxes in the soil-plant-atmosphere continuum under study (Figs. 10, 11).

Figure 10 Seasonal Origin Index of transpiration (SOI_T) and bottom fluxes (SOI_Bot), simulated with HYDRUS-1D, during the growing seasons from 2018 to 2022.

Figure 11 (a) Equivalent precipitation (P_eq). (b) Evaporation flux. (c) Transpiration flux. (d) Actual flux across the bottom of the soil profile. (e) Monitored spring (SOU) discharge. The red color indicates SOI > 0 (summer water is overrepresented in the flux), while the blue color indicates SOI < 0 (winter water is overrepresented in the flux). In panel (a) equivalent precipitation and the isotopic composition of equivalent precipitation ($C_{P_{\mathrm{eq}}}$), from which the SOI has been calculated, are obtained with the Ceperley et al. (2020) model. In panels (b), (c) and (d) the fluxes and their isotopic composition, used to retrieve the SOI, are simulated by using HYDRUS-1D. In panel (e) the SOU isotopic composition is derived from measurements: in order to have a continuous isotopic composition (and consequently continuous SOI) at all time-steps we fit a sine function on data as described in von Freyberg et al. (2018).

Figure 11a presents the simulated equivalent precipitation (P_eq) and its SOI. During snowmelt period (mid-April to mid-June), SOI of P_eq (SOI${}_{P_{\mathrm{eq}}}$) exhibits value close to −1, indicating a dominant contribution from the snow accumulated during winter. In contrast, during late summer (July to September), SOI${}_{P_{\mathrm{eq}}}$ exhibits value close to 1, indicative of summer rainfall-dominated inputs.

Evaporation is sourced from winter precipitation (i.e., snowmelt) during the snowmelt period and from summer rainfall during the later months (Fig. 11b). Indeed, snowmelt recharges the soil predominantly between mid-April to mid-June, while rainfall inputs prevail from July onward.

Transpiration, on the other hand, is largely fed by summer rainfall, as reflected by SOI of transpired water (SOI_T) greater than 0 during the core of the growing season (Figs. 11c, 10). This is also evident from the sink term (S), reported in Fig. 12e, which peaks from July to August when summer precipitation dominates soil profile inputs, except in 2022 (Fig. 12a). These findings are in line with empirical evidence from the Matsch/Mazia catchment in the eastern Alps, where springtime snowmelt does not coincide with peak vegetation activity, leading plants to rely primarily on summer rainfall (Zuecco et al., 2024). Furthermore, our results confirm the finding by Nehemy et al. (2022) that snowmelt can contribute to transpiration early in the growing season, albeit over a brief period, as evidenced from negative SOI_T at the end of May (Fig. 11c). This is also evident in the early increase of the S at the end of May (Fig. 12e), coinciding with winter-sourced recharge (Fig. 12a), though values remain below the S peak seen in July–August.

Figure 12 (a) Seasonal Origin Index (SOI). (b) Volumetric water content (θ). (c) Hydraulic conductivity (K). (d) Darcy velocity (V, the sign convention is positive upwards and negative downwards). (e) Sink term over time and soil depth. Please, note that only 138 print times from 1 November 2017 to 6 February 2023 with a 14 d time step are reported. In panel (a) the red color indicates SOI > 0 (summer water is overrepresented in the soil water), while the blue color indicates SOI < 0 (winter water is overrepresented in the soil water). All the variables over time and soil depth have been simulated with HYDRUS-1D.

The intense and sustained meltwater inputs lead to the saturation of the soil profile (Fig. 12c) and correspond to the highest modeled Darcy velocities (Fig. 12d), suggesting enhanced deep filtration of the snowmelt. This is supported by previous findings revealing that snowmelt is generally more effective than rainfall in filtrating beyond the root zone (Earman et al., 2006). Indeed, the SOI of the bottom flux (SOI_Bot), representing water that contributes to groundwater recharge, exhibits values clearly below zero between May and June (Fig. 11d), indicating a substantial influence of snowmelt during this period. This is consistent with findings from other snow-dominated catchments, where groundwater typically reflects the isotopic composition of snowmelt (Michelon et al., 2023; Pavlovskii et al., 2018). Moreover, past studies revealed that seasonally snow-covered catchments resulted in a snowmelt pulse that enables high groundwater recharge (Ajami et al., 2012; Harrison et al., 2021; Hotovy et al., 2025; Winograd et al., 1998) during summer (Cochand et al., 2019; Du et al., 2019; Flerchinger et al., 1992; Hayashi, 2020). Further support comes from comparing the SOI_Bot to that of the monitored spring (Fig. 11d, e). Despite the bottom flux shows summer signatures (SOI_Bot>0) during wintertime (October to February) events (Fig. 11d), the SOI of spring water remains slightly lower than 0 (Fig. 11e) in this period, suggesting subsurface mixing with winter-sourced storage, likely recharged by snowmelt during the preceding summers.

The previous observations point to a vertical connectivity within the soil profile during snowmelt peaks, where infiltrating water rapidly fills available pore space and microtopographic storage, leading to sudden increases in vertical subsurface flow: a process consistent with the fill-and-spill conceptual model at the plot scale (McDonnell et al., 2021). In this regard, we infer a possible fill-and-spill mechanism at the plot scale (McDonnell et al., 2021). This insight is supported by the timing of peak bottom fluxes, which align with elevated water content and hydraulic conductivity across all soil depths (Fig. 12b and c), and may be further facilitated by preferential flow pathways in microporous grassland soils (Mohammed et al., 2019).

It should be noted that, in steep slope contexts, infiltration occurring upslope could be redistributed laterally at shallow depths (<60 cm) and subsequently re-emerge downslope supplying again the root zone. In these specific contexts, infiltration at the bottom of the soil column does not necessarily lead to groundwater recharge or streamflow contribution. Accordingly, we cannot entirely exclude the presence of such process. Nevertheless, this mechanism can reasonably be considered less relevant at the study site for three main reasons.

First, although the hillslope hosting the monitoring station has an average slope of approximately 32° (Gisolo et al., 2022), the instruments are installed on a small, relatively flat plateau. This local topographic configuration supports the assumption that water fluxes within the monitored soil profile are predominantly vertical and can thus be reasonably represented using a 1D modeling approach.

Second, independent evidence at the catchment scale suggests that vertical processes dominate over lateral transfers. In Gentile et al. (2023), the fraction of young water (i.e., water younger than 2–3 months) in both streamflow and spring water was estimated to be low (0.18 and 0.11, respectively). These values are consistent with findings from mountainous catchments worldwide (Jasechko et al., 2016), suggesting that, although seemingly counterintuitive, steeper catchments tend to favour deeper vertical infiltration rather than shallow lateral flow.

Third, the comparison between the isotopic composition of the modeled bottom flux and that of spring water (Fig. 11d–e) further supports the dominance of vertical processes over lateral transfers. This interpretation is consistent with isotope-based conceptualizations presented in Gentile et al. (2023) and with global-scale findings (Jasechko, 2019), which show that event water generally represents only a minor fraction of streamflow, while discharge is largely dominated by pre-event water. These results suggest that equivalent precipitation (rainfall + snowmelt) primarily infiltrates into the subsurface, displacing pre-event water that is subsequently released to the stream. This appears to be the dominant process in the study area.

From Figs. 11c, d and 10, in which is reported the SOI_T and SOI_Bot during the growing season (1 May–30 September), it is possible to observe a high degree of ecohydrological separation of this high-elevation grassland in the years 2018 to 2021. Winter precipitation (e.g., snowmelt) mainly constitutes a mobile water pool that rapidly recharges the groundwater storage (Earman et al., 2006), which in turn supplies streams, so that snowmelt is generally poorly available for plant transpiration. The latter is mainly supplied by a less mobile water pool constituted by summer rainfall that remains available in the soil profile during the core of the growing season. Indeed, the sustained snowmelt pulse was sufficient to saturate the soil, thus explaining the development of vertical connectivity among soil pores. In contrast, summer rainfall events are typically intermittent and, if not intense enough to saturate the soil, they do not generate vertical pore connectivity. This pattern is consistent with Radolinski et al. (2021) asserting that the TWW can occur after intense events.

Interestingly, deviations from this pattern were observed during the 2022 snow drought (Fig. 10). The average SOI_T resulted lower (and negative) than the previous years, thus highlighting a greater snowmelt contribution to these fluxes (Fig. 11c). This finding partially aligns with Mastrotheodoros et al. (2020) who observed enhanced AET from earlier snowmelt during the 2003 Alpine drought, thus implying a reduction of groundwater recharge. In our study site, this can be explained by considering an early (i.e., during winter/spring) and more intermittent (i.e., less concentrated and intense) meltwater input with a consequent low hydraulic conductivity K and Darcy velocity along the soil profile during the summer of 2022 (Fig. 12c, d). Most likely, these conditions strongly limit the vertical pore connectivity along the soil profile with a consequent drastic reduction of the bottom flux and with consequent low groundwater storage recharge. Thus, we can observe an enlargement of the time window in which the infiltrated snowmelt is retained in the soil for supplying transpiration along with summer rainfall (Fig. 12a). This finding, together with the evidence that snowmelt contributes to transpiration during the early stages of the growing season (Fig. 11c), confirms that ecohydrological separation should not be viewed as strict duality, but rather as a matter of “degree of separation” that depends on time-variable and site-specific hydrological dynamics (Kirchner et al., 2023).

Concluding, in light of previous studies suggesting that plants can access soil water disconnected from groundwater and streamflow, we have tested the following null hypothesis (H₀): “Winter precipitation (i.e., snowmelt) rapidly transits the soil profile recharging groundwater and streams, while summer precipitation (i.e., rainfall) remains available to sustain transpiration fluxes” which describes a seasonal nature of the TWW hypothesis. Considering our results, we cannot conclusively accept or reject this hypothesis, as it is framed in a strongly dichotomous manner. In this regard, our findings support the view that it is more appropriate to refer to a degree of ecohydrological separation, which, at our study site, appears more evident in years when snowmelt input is concentrated and continuous, and less pronounced during dry years – when snowmelt input is earlier and intermittent.

4 Conclusions

This study provides new insights into the seasonal partitioning of water resources in a high-elevation alpine grassland, with a particular focus on the degree of ecohydrological separation between seasonal water pools supplying plant transpiration and groundwater recharge. By integrating a snow isotope model with the HYDRUS-1D model, key water fluxes and their isotopic compositions within the soil–plant–atmosphere continuum were simulated under contrasting hydrometeorological conditions.

The results demonstrate that, during years with concentrated and sustained snowmelt inputs, a pronounced ecohydrological separation emerges: the snowmelt, due to soil saturation and the possible generation of a vertical pore connectivity, rapidly drains beyond the root zone, contributing to recharge, while summer rainfall is retained in the soil and primarily used by vegetation. However, during the 2022 snow drought, reduced and more intermittent snowmelt inputs led to lower soil saturation and possibly limited the vertical pore connectivity, enabling winter-sourced snowmelt water to remain accessible to plants for a longer period. This shift resulted in a reduced degree of separation between seasonal water pools, as also indicated by the SOI values of transpiration and bottom fluxes.

While the adopted modelling framework enabled a process-based interpretation of seasonal compartmentalization of water fluxes, some limitations should be acknowledged. The degree-day model used in this work provides a simplified representation of snowmelt processes: compaction, refreezing, stratification, and snow sublimation are not explicitly included. These simplifying assumptions propagate into the estimation of both the timing and magnitude of meltwater infiltration. In addition, the adopted snow isotope model does not account for the effects of refreezing or sublimation processes on the isotopic composition of meltwater. Addressing these processes would require dedicated efforts in model development. Further sources of uncertainty arise from the relatively low temporal resolution of isotope sampling, which is largely constrained by the intrinsic challenges of complex high-elevation mountain environments, as well as from potential biases associated with the cryogenic vacuum distillation, both of which may affect the estimation of longitudinal dispersivity and the representation of subsurface transport processes.

Despite these limitations, the findings suggest that the Two Water Worlds (TWW) hypothesis, while useful as a conceptual framework, may oversimplify the dynamic nature of subsurface water partitioning. Rather than a strict duality, our results support the interpretation of ecohydrological separation in this mountain environment as a continuum, with its magnitude modulated by time-variable seasonal water input and root water uptake patterns.

Given projected changes in snow regimes under climate warming, including an increasing frequency of snow droughts, the variability in the degree of ecohydrological separation observed in this study has important implications for anticipating future shifts in ecohydrological functioning within similar mountain ecosystems worldwide.

In this regard, incorporating climate change projections into the modelling framework of this high-elevation grassland would be crucial to provide new insights into future shifts in the relative contribution of winter- versus summer-derived water to transpiration and groundwater recharge.

Appendix A: List of Symbols

Symbol	Description
a	Slope of the Local Meteoric Water Line
AET	Actual evapotranspiration (mm d−1)
AET H-1D	HYDRUS-1D derived actual evapotranspiration (mm d−1)
AET Eddy	AET derived from Eddy-covariance measurements (mm d−1)
b	Intercept of the Local Meteoric Water Line
BC	Boundary condition
C	Isotopic composition (‰)
$C_{P_{\mathrm{eq}}}$	Isotopic composition of equivalent precipitation (‰) obtained with Ceperley et al. (2020) model
CP-fit	Isotopic composition of precipitation derived by fitting a sine curve to observed data
CS	Isotopic composition of the snowpack water (‰) obtained with Ceperley et al. (2020) model
CVD	Cryogenic Vacuum Distillation
D	Dispersion coefficient (cm2 d−1)
Dw	Molecular diffusion coefficient in free water (m2 s−1)
DOR	Dora del Nivolet catchment
ElevS	Monitoring station elevation (m a.s.l.)
EP	Potential evaporation (mm d−1)
ES	Evaporation slope (–)
ETP	Potential evapotranspiration (mm d−1)
HC	Empirical coefficient of the Modified Hargreaves–Samani equation
HE	Empirical exponent of the Modified Hargreaves–Samani equation
HS	Hargreaves–Samani equation
HSM	Modified Hargreaves–Samani equation
HT	Factor used to convert units from Fahrenheit to Celsius in the HSM
IRMS	Isotope Ratio Mass Spectrometer
k	Extinction coefficient for global solar radiation within the canopy (–)
K	Unsaturated hydraulic conductivity (cm d−1)
Ks	Saturated hydraulic conductivity (cm d−1)
LAI	Leaf Area Index (m2 m−2)
MAE	Mean Absolute Error
MES	Mean Evaporation Slope (–)
PR	Rainfall (mm d−1)
PS	Snowfall (mm d−1)
Peq	Equivalent precipitation (Rainfall + Snowmelt, mm d−1)
Ra	Extraterrestrial radiation (mm d−1)
Rsample	Isotopic ratio (${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}$O) of the water sample
Rstandard	Isotopic ratio (${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}$O) of the reference standard (V-SMOW)
RMSE	Root Mean Square Error
S	Root water uptake sink term (d−1)
Su	Correction factor for including measurement uncertainty in line-conditioned excess calculation
SD	Standard deviation
SD${}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$	Standard deviation of δ18O associated with the isotopic analysis method (‰)
SD${}_{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}$	Standard deviation of δ2H associated with the isotopic analysis method (‰)
SM	Snowmelt (mm d−1)
SOI	Seasonal Origin Index (–)
SOIQ	Seasonal Origin Index of streamflow (–)
SOIAET	Seasonal Origin Index of evapotranspiration (–)
SOIBot	Seasonal Origin Index of bottom flux (–)
SOIPeq	Seasonal Origin Index of Peq (–)
SOIT	Seasonal Origin Index of transpiration flux (–)
SOU	“Source”: spring within the Dora del Nivolet catchment

Symbol	Description
SWE	Snow water equivalent (mm)
T	Air temperature (°C)
Tmax	Daily maximum air temperature (°C)
Tmean	Daily mean air temperature (°C)
Tmin	Daily minimum air temperature (°C)
T0	SM begins once T surpasses a defined melting threshold T0 (°C)
TR	Precipitation is classified as PR when T exceeds an upper threshold TR (°C)
TS	Precipitation is classified as PS when T is below a lower threshold TS (°C)
TS	Trendline Slope (–)
TP	Potential transpiration (mm d−1)
TWW	Two Water Worlds
h	Pressure head (cm)
lc-excess*	Line-conditioned excess that accounts for uncertainty in the isotopic analysis
LMWL	Local Meteoric Water Line
m	van Genuchten shape parameter (–), $m=\mathrm{1}-\mathrm{1}/n$
n	van Genuchten shape parameter (–)
q	Soil water flux (cm d−1)
t	Time (days)
z	Depth below soil surface (positive upward, cm)
α	van Genuchten shape parameter (cm−1)
δ2H	Deuterium isotopic composition (‰)
δ18O	Oxygen-18 isotopic composition (‰)
δx	Fractionation-compensated isotopic composition of the considered flux (‰)
δannP	Volume-weighted annual precipitation isotopic composition (‰)
δsummerP	Isotopic composition of typical summer precipitation (‰)
δwinterP	Isotopic composition of typical winter precipitation (‰)
ξ	Degree-day factor (mm °C−1 d−1)
r	Pearson's Linear Correlation Coefficient
θ	Volumetric water content (cm3 cm−3)
θr	Residual volumetric water content (cm3 cm−3)
θs	Saturated volumetric water content (cm3 cm−3)
λL	Longitudinal dispersivity (cm)
τw	Tortuosity factor (–)