Isotopic insights into the dynamics of soil water pools along an elevation gradient

Recent intensive research on the soil–plant–atmosphere continuum has introduced novel methodological approaches. These include new in-situ extraction techniques and the application of stable hydrogen and oxygen isotopes in water enabling to trace water movement and plant responses at much finer spatial and temporal scales. Such approaches provide detailed insights into soil water dynamics and plant adaptation to changing environmental conditions under climate change. This study aims to provide a comprehensive characterization of dynamics of distinct soil water pools – mobile versus tightly bound water – along an elevation gradient, while simultaneously assessing the impact of the absent snow accumulation in lowland areas on soil water distribution compared to higher elevations. In contrast to conventional bulk water sampling, a key innovation of this study lies in the experimental design across the elevation gradient combined with a novel extraction method that selectively isolates tightly bound soil water for isotopic analysis. The results indicate a prolonged residence time of winter-derived soil water in lowland sites, in contrast to a rapid turnover at the highest elevation, where the winter water signal dissipates shortly after snowmelt. Distinct isotopic compositions among water pools – mobile versus tightly bound water – were particularly evident in lowland areas at the edges of the growing season (up to 3 ‰ and 21 ‰ for δ¹⁸O and δ²H, respectively), while tightly bound and bulk soil water exhibited – on average – only minor or no isotopic differences. In the context of the projected continued decline in snow cover at higher elevations in Central Europe, these findings are critical for improving predictions of soil water storage and, consequently, plant water availability under ongoing climate change.

1 Introduction

Soil drought is becoming increasingly prevalent due to climate change, which alters air temperature, total amount of precipitation and its intra-annual distribution (Gebrechorkos et al., 2025; Samaniego et al., 2018). These shifts have contributed to a sustained decline in vegetation-accessible water over the past three decades (Jiao et al., 2021). In response, there has been growing interest in the role of snowpack water storage and runoff generation in snow-dominated catchments, which are essential for groundwater recharge and soil moisture replenishment (Jenicek and Ledvinka, 2020; Jenicek et al., 2021; Šípek et al., 2021; Musselman et al., 2017). Numerous studies project a continued decline in snow cover across mountainous regions as a consequence of rising air temperatures (Musselman et al., 2017; Marty et al., 2017; Jenicek et al., 2018; Willibald et al., 2020), accompanied by a shift from snowfall to rainfall during the winter season (Harpold et al., 2017; Safeeq et al., 2016). The impacts of these changes on snow storage for the annual water balance represent a critical and unresolved question in hydrological research. Equally important issues are the downstream consequences for plant-available water in lowland ecosystems, particularly during the late part of the growing season when drought stress is most acute (Büntgen et al., 2021; Qin et al., 2020; Mankin et al., 2019). A thorough understanding of the soil–plant–atmosphere continuum, a concept originally introduced by Gradmann (1928) and later formalized by van den Honert (1948), is therefore crucial for predicting vegetation dynamics and adaptive responses under increasingly frequent and severe drought periods.

The relationship between plant water use and local hydrology has been studied since the early 20th century (Bates, 1921). Pioneering studies on water transport through soils and plants were summarized in comprehensive reviews (e.g., Tinker, 1976; Weatherley, 1976; Molz, 1981). A major shift in perspective occurred when Dawson and Ehleringer (1991) demonstrated that some riparian trees primarily access deeper groundwater, rather than the more readily available stream water. A decade later, Bond et al. (2002) appeared to challenge this finding by demonstrating diel fluctuations in stream baseflow attributable to plant transpiration, demonstrating clear interactions between transpiration and streamflow. Despite this apparent contradiction, Brooks et al. (2010) showed that mobile and tightly bound soil water (represented by the stream and plant water, respectively) are isotopically distinct. Their results suggested that, especially during the dry season, mobile water traveling through macropores or pipes bypasses tightly bound soil water, which is instead more likely to be taken up by plants and not contribute to streamflow.

This conceptual breakthrough formed the basis for the ecohydrological separation framework, later termed the “Two Water Worlds” (TWW) hypothesis (McDonell, 2014). Since then, the TWW hypothesis has stimulated widespread debate, with numerous studies supporting (e.g., Goldsmith et al., 2012; Evaristo et al., 2015; Hervé-Fernandez et al., 2016) or challenging (e.g., Geris et al., 2015; Vargas et al., 2017; Dubbert et al., 2019) the existence of isotopically distinct water pools for vegetation use and runoff generation. These contrasting findings prove, that the hydrological connectivity between plant-accessible water and mobile water remains a central, unresolved question in ecohydrology and an issue of potential methodological and conceptual limitations in current approaches. Consequently, a new way forward has been proposed (Berry et al., 2017), emphasizing the need to investigate internal water cycling within the phloem and xylem, identifying potential sampling and methodological biases, and improve both the spatial and temporal resolution of sampling strategies.

However, an overwhelming majority of studies comparing soil water and xylem water (e.g., Zapater et al., 2011; Meunier et al., 2017; Vargas et al., 2017; Barbeta et al., 2019, 2020; Liu et al., 2021; Brighenti et al., 2024; Benettin et al., 2024) rely on mobile and so-called bulk soil water to be compared. Mobile water is typically extracted using suction lysimeters or other vacuum-based systems. The water obtained in this manner – usually under tension of −60 kPa (Brooks et al., 2010; Muñoz-Villers and McDonnell, 2012; Berry et al., 2017; Sprenger et al. 2018; Haagsma et al., 2024) – originates primarily from macropores and preferential flow paths. In contrast, bulk soil water encompasses the total soil water content, including both gravitational and capillary pore water. During dry periods, when macropores are emptied and suction lysimeters cannot collect water, bulk soil water reliably represents the tightly bound water fraction. However, during wet periods, when gravitational pores are partially or fully saturated, a significant portion of the bulk soil water may also consist of mobile water. In this case, bulk water may no longer represent the tightly bound fraction and may not be as suitable for direct comparison with xylem water, particularly under the assumption that plants preferentially utilize tightly bound soil water (McDonell, 2014). In our study, mobile water (MW) and tightly bound water (TBW) are defined consistently with previous work mentioned above: MW is extracted at −60 kPa, while TBW is defined as the soil water remaining after the removal of MW. Bulk water (BW) is considered to represent the total soil water content.

The aim of this study is therefore to determine a distinction between mobile soil water and experimentally extracted TBW, which is mostly examined only as a component of BW. We hypothesize that TBW differs in its isotopic composition from both MW and BW, and that this distinction persists for longer periods at higher elevations characterized by greater snow accumulation. To test this hypotheses, the experiments were conducted simultaneously at four sites spanning an elevation gradient exceeding 1000 meters. This elevation-focused design is particularly relevant in the rain-snow transition zone, which is highly sensitive to temperature-driven changes in snow storage. A site selection was considered to locations with similar soil texture, given its strong influence on the relative abundance of macropores and capillary pores. Soil water sampling was carried out at two-week intervals from February to November 2023, except at the highest-elevation site, which was accessible only from May onwards.

The primary objectives of the study are to:

1. determine whether tightly bound and mobile soil water differ significantly in their isotopic composition;

2. characterize the dynamics of the isotopic composition of the soil water along the elevation gradient in relation to declining snow cover and overall precipitation;

3. identify the dominant sources of tightly bound soil water; and

4. assess whether substituting bulk soil water with tightly bound soil water alters lead interpretations of soil water sources and their seasonal dynamics.

Table 1 Detailed characteristics of selected experimental areas. Climatic data (total annual precipitation and mean annual temperature) are for 2023. Snow data and climate classification (Köppen system) follow Tolasz et al. (2007). Further details for individual locations are provided in the cited references.

Figure 1 Location of selected experimental areas. Panels: top left: Zvěřínek; top right: Trhové Dušníky; bottom left: Liz; bottom right: Rokytka. Map data: Digital Vector Database of the Czech Republic ArcČR^® version 4.3 (ARCDATA PRAHA, s.r.o., 2024).

2 Study sites

The experiment was conducted in Czechia at four experimental plots (Table 1 and Fig. 1) strategically selected to span a considerable elevation gradient covering corresponding variations in temperature, precipitation, and snow cover extent and duration. The study sites ranged from the fertile agricultural lowlands of the Elbe River Basin – Zvěřínek (185 m a.s.l.) – through the highlands and foothills – Trhové Dušníky (430 m a.s.l.) and the Liz catchment (870 m a.s.l.) – to the upper montane zone of the Bohemian Forest – Rokytka catchment (1260 m a.s.l.). For simplicity, these locations are referred to by the abbreviations ZV, TD, LI, and RO, respectively.

Although the soil types vary due to the natural pedogenetic context of each site (Regosol, Gleyc Fluvisol, Cambisol, and Podzol), all plots exhibit similar soil texture (loamy sand or sandy loam). The minimal differences in clay content provide a more direct assessment of elevation-related influences on soil water behaviour. Due to considerable elevation differences, it was not possible to maintain identical vegetation cover across all sites. However, vegetation at each location represents typical plant communities at the corresponding elevation zones in Central Europe – ranging from agricultural land in the lowlands, to meadows, spruce forest, and beech forest at higher elevations. The agricultural land at the lowest elevation (ZV) lacked vegetation cover from February to March, with post-harvest stubble remaining from October onward.

Throughout the year, groundwater levels at most sites remain at least four meters below the soil surface. The TD site represents the only exception, for wich groundwater rises to approximately one meter below the surface during the spring months, potentially influencing the isotopic composition of the overlying soil profile. However, due to the sandy texture of the soil, capillary rise is limited, and this influence is therefore assumed to be confined to the lower soil layer sampled.

Each site was equipped with a meteorological station (Fiedler, Czechia) with rainfall gauges (Meteoservis, Czechia), Palmex precipitation collectors (Palmex Ltd., Zagreb, Croatia), and a tensiometer-regulated vacuum lysimeter system (VS-Pro, UMS, Germany) for MW sampling. Precipitation amounts, soil moisture, and air temperature were measured in 10 min intervals and calculated on daily basis.

3 Data and methods

3.1 Field sampling

At all study sites, precipitation was sampled at two-week intervals, except at the LI site, where samples were collected up to three times per week. Snow cover was sampled when present; however, the lower-elevation sites experienced sparse or no snow cover. For the comparison between precipitation and soil water, only precipitation samples from the respective hydrological year were used, as older water is not expected to persist at such shallow depths.

Furthermore, MW samples and soil cores for TBW extraction were collected at two-week intervals for the stable isotope analysis from February to November 2023. The exception was the RO location, accessible only from May 2023 onwards due to its remote location in the heart of the Bohemian Forest National Park, as well as heavy snow conditions during winter and early spring.

MW was extracted using a tensiometer-controlled vacuum system (VS-Pro) with a maximum applied tension of -60 kPa (Brooks et al., 2010; Berry et al., 2017; Sprenger et al., 2018). Extraction was conducted at two depths of 20 and 40 cm (hereafter referred to as D20 and D40, respectively); at the TD site, an additional depth of 60 cm was included. These depths often play a significant role in root water uptake (Hackmann et al., 2026). Samples collected in this manner represents a composite of water collected over the preceding two weeks.

For TBW extraction, undisturbed soil cores (100 cm³) were collected from the same depths, with five replicates per depth. The soil cores were wrapped in Parafilm^® and stored in a portable refrigerator during transport to the laboratory. In total, 805 soil cores, 320 MW samples, and 195 precipitation samples were collected during the sampling period (Table 2).

Table 2 Overview of samples collected throughout the sampling period.

The Type column indicates sampling depth (cm). For soil cores, the number of locations corresponds to the number of sampling campaigns (five samples per depth from a single soil profile per campaign). For mobile water, it represents the number of suction cups installed at the respective depths. For groundwater, it indicates the number of wells sampled.

3.2 Laboratory processing of soil samples

To obtain TBW, intact soil cores were placed in a pressure plate apparatus for the determination of the soil water retention curve (5 bar Pressure Plate Extractor, Soil Moisture Equipment Corp., CA, USA), following a modified protocol based on Orlowski and Breuer (2020). A pressure of 60 kPa (∼ pF 2.4) was applied using a 1 bar pressure plate cell for a two-week period to remove the MW fraction. Consequently, the top and bottom sections of each core were excised to eliminate portions potentially affected by contact with the ceramic plate and by evaporative losses. Only the central section of the core was retained for further extraction, yielding approximately 50 g of soil material per sample.

Figure 2 Sample processing procedure during the mass-balance mixing method.

For the subsequent extraction of TBW, the mass-balance mixing method was chosen due to its accessibility, simplicity and high throughput. A comprehensive description of the method and validation through spike experiments for the investigated soil type is presented in the Sects. S1 and S2 in the Supplement. Briefly, the experiments were carried out in the following steps:

• The soil samples (20–30 g) were transferred into 40 mL glass vials (ND24 (EPA), VERKON s. r. o., Czech Republic) equipped with plastic screw caps fitted with silicone/PTFE septa (VERKON s. r. o., Czech Republic).

• 20–25 mL of isotopically labelled water (distilled tap water) was added to each sample, and the remaining headspace was filled with glass beads (5 mm in diameter) to eliminate residual air and to enhance mechanical disaggregation of soil aggregates, thereby promoting efficient mixing (Fig. 2).

• The vials were mounted on a laboratory-constructed rotating device (Fig. S2 in the Supplement) and continuously rotated for 16 h at a constant speed of 15 rpm.

• The samples were refrigerated and left to settle to allow sedimentation.

• A 0.75 mL sample of the clear liquid phase was collected and filtered through a 0.45 µm mixed cellulose ester syringe filter. The remaining sample was oven-dried at 105 °C for 48 h and weighed to determine the soil water content of the soil sample.

3.3 Isotope analysis and calculations

Stable isotope analyses were performed at the Institute of Hydrology (Czech Academy of Sciences) with a L2140-i isotope analyser (Picarro Inc., Santa Clara, US). Standard mode (precision of ±0.03 ‰ and ±0.15 ‰ for δ¹⁸O and δ²H, respectively) was used with 6 injections per sample, out of them the first 3 injections were discarded. The isotope ratios are reported in per mil (‰ ) relative to Vienna Standard Mean Ocean Water (VSMOW) (δ²H or δ¹⁸O = $(R_{\mathrm{sample}}/R_{\mathrm{standard}}-\mathrm{1})\times \mathrm{1000}$ ‰, where R_sample is the isotope ratio of the sample and R_standard is the known reference value (i.e., VSMOW), see Craig (1961).

The isotopic composition of TBW was calculated according to the mass balance mixing model (Eqs. 1, 2).

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{TBW}}={\displaystyle \frac{m_{\mathrm{M}}}{m_{\mathrm{TBW}}}}\cdot {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{M}}-{\displaystyle \frac{m_{\mathrm{L}}}{m_{\mathrm{TBW}}}}\cdot {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{L}}\text{(2)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{TBW}}={\displaystyle \frac{m_{\mathrm{M}}}{m_{\mathrm{TBW}}}}\cdot {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{M}}-{\displaystyle \frac{m_{\mathrm{L}}}{m_{\mathrm{TBW}}}}\cdot {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{L}}\end{array}$$

where the sub-indices M, TBW and L denote the mixture, tightly bound soil water, and isotopically labelled water, respectively; m represents the mass of water in each pool (determined gravimetrically using a pressure plate apparatus), with the total mass of the mixture given by $m_{\mathrm{M}}=m_{\mathrm{TBW}}+m_{\mathrm{L}}$. The terms δ¹⁸O and δ²H denote the stable isotopic composition of the respective water pools.

Following the estimation of TBW, the equivalent mass-balance mixing model was applied to calculate the potential stable isotopic composition of BW. In this approach, BW was represented as a mixture of MW, obtained from suction lysimeters, and TBW, as derived in the previous step. The relative proportions of these two components were determined based on measurements from the pressure plate apparatus.

To establish the local meteoric water lines (LMWLs) and regression lines for soil water isotopic compositions, a reduced major axis (RMA) regression was applied (Harper, 2016) instead of the conventional linear regression. A necessity of the application of the RMA regression is caused by significant uncertainties in both δ²H and δ¹⁸O measurements, moreover in the case of soil water, substantially increased by the mixing-extraction method, systematically introducing a shift. In contrast, isotopic measurements of precipitation were associated primarily with analytical precision only. Since the classical least squares regression a priori assumes that all errors are related to the y axis only and entirely ignores uncertainties along the x axis, the RMA regression was considered appropriate for both soil water and precipitation, providing a more accurate representation of the δ²H–δ¹⁸O relationship by minimizing the orthogonal distances between data points and the fitted approximation.

The slope β_RMA (Eq. 3) of the RMA regression was calculated as the ratio of the standard deviations of δ²H and δ¹⁸O, multiply by the sign of their Pearson correlation coefficient, and the intercept αRMA (Eq. 4) was determined as the difference between the mean δ²H and the product of the slope and mean δ¹⁸O:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{RMA}}=\text{sign}\left(r\right)\cdot {\displaystyle \frac{{\mathit{\sigma }}_{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}}{{\mathit{\sigma }}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}}}\text{(4)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{RMA}}=\stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}-{\mathit{\beta }}_{\mathrm{RMA}}\cdot \stackrel{\mathrm{‾}}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}\end{array}$$

where r is the Pearson correlation coefficient between δ²H and δ¹⁸O, and ${\mathit{\sigma }}_{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}$ and ${\mathit{\sigma }}_{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}$ are their standard deviations. This method minimizes the orthogonal distances between data points and the fitted line, making it suitable for hydrological isotope data where both variables are subject to analytical and natural variability.

To quantify the uncertainty of the RMA regression predictions, a bootstrap procedure with 10 000 resamples was applied. For each bootstrap iteration, a dataset of equal size was sampled with replacement, the RMA regression coefficients were recalculated, and residual variation was added to generate a distribution of predicted δ values. 95 % confidence intervals were derived at each δ¹⁸O value from the 2.5th and 97.5th percentiles of these predicted values, providing prediction intervals for the regression lines.

The slopes of the regression lines were compared among different data sets (precipitation, MW, TBW) using their bootstrap distributions. For each comparison, differences between bootstrap slopes were calculated, and 95 % confidence intervals of these differences were obtained. A slope difference was considered statistically significant if the 95 % interval did not include zero. Additionally, pseudo-R² values were calculated as the squared Pearson correlation coefficient for each dataset within the aim to indicate the goodness-of-fit.

Furthermore, the line-condition excess (lc-excess; Landwehr and Coplen, 2006) was determined (Eq. 5) to identify and exclude data contaminated during the water extraction process.

$$\begin{array}{}\text{(5)}& \text{lc-excess}={\mathit{\delta }}^{\mathrm{2}}\mathrm{H}-a\cdot {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}-b\end{array}$$

where a and b are the coefficients of the LMWL from the individual experimental plots. This contamination was manifested by abnormally high lc-excess values relative to the rest of the dataset. An lc-excess value of 10 was taken as a pragmatic screening threshold; however, this criterion was not applied rigidly. Values slightly exceeding 10 were retained in the final dataset when the overall distribution of the respective sample group remained below or close to this threshold. In contrast, apparently contaminated samples, typically exhibiting lc-excess values between 20 and 50, were removed entirely. In such cases, all associated measurements were excluded, as well as occasional values marginally below 10, when they represented isolated initial observations followed by consistently elevated lc-excess values within the same sample set.

The methodological nature of the error was further supported by the observation that these anomalous values appeared randomly across different sites and sampling dates, with the only consistent factor being that the affected samples were processed together within the same run of the overpressure apparatus. In total, 6 out of 32 extraction runs of TBW were discarded due to this methodological error and significant sample contamination.

3.4 Seasonal Origin Index (SOI)

To characterize whether the extracted soil water originated from winter or summer precipitation, the Seasonal Origin Index (SOI) was calculated (Eq. 6; Allen et al., 2019) for individual sampling dates, to provide a more detailed representation of gradual changes in water origin resulting from mixing with newly infiltrating water.

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

where δ_x denotes the δ¹⁸O isotopic composition of soil water, and δ_annP, δ_winterP, and δ_summerP represent the δ¹⁸O values of volume-weighted annual precipitation, the characteristic winter (δ_annP− fitted amplitude) and summer (δ_annP+ fitted amplitude) precipitation, respectively. These precipitation components were derived from sinusoidal fitting to the time series. Prior to analysis, samples with negative lc-excess were corrected according to Benettin et al. (2018) to account for potential biases in isotope measurements.

Specifically, a sine function was fitted to the isotope data using an iteratively reweighted least squares (IRLS) regression approach with externally supplied weights for precipitation (following the methodology of Kirchner, 2016). This approach is widely applied to quantify the seasonal isotope cycle in precipitation, streamflow, soil water, and groundwater. The fitted amplitude was subsequently obtained from Eq. (7) and calculated according to Eq. (8).

$$\begin{array}{}\text{(7)}& {\displaystyle }c\left(t\right)=a\cdot \cos \left(\mathrm{2}\mathit{\pi }ft\right)+b\cdot \sin \left(\mathrm{2}\mathit{\pi }ft\right)+k\text{(8)}& {\displaystyle }A=\sqrt{a^{\mathrm{2}}+b^{\mathrm{2}}}\end{array}$$

where t is the time, c(t) represents the isotopic time series of the dataset, a and b are the cosine and sine coefficients determined by the IRLS regression, f is the frequency of annual isotopic fluctuation (f=1 yr⁻¹ for a seasonal cycle), k is the vertical shift of the sine wave, and A is the amplitude of the fitted sine wave.

Figure 3 Dual-isotope plots of all water samples collected in this study, with corresponding regression lines. Panels: top left: ZV (lowland agricultural field); top right: TD (mid-elevation meadow); bottom left: LI (submontane spruce forest); bottom right: RO (montane beech forest). Precipitation is shown in blue, mobile soil water in green, and tightly bound soil water in brown. Blue shading shows the 95 % prediction interval for precipitation.

4 Results

4.1 Precipitation and soil water data

The collected data closely follow the LMWLs of the respective locations, however, a clear difference is observed for TBW, which follows the LMWL with greater dispersion. In comparison with MW, several samples lie outside the 95 % prediction interval for precipitation (Fig. 3), indicating the influence of isotopic fractionation associated with evaporation and subsequent condensation and internal mixing processes within the soil matrix.

Table 3 RMA regression parameters (slope and intercept) along with their corresponding standard deviations (SD), estimated using bootstrap analysis.

Due to the pronounced elevation gradient between the study sites, a systematic shift in both the slope and intercept of the LMWLs was observed. The ZV lowland agricultural site (on average 150 m a.s.l.) is characterized by a reduced slope (7.6) and a low intercept (3.7), indicative of enhanced sub-cloud evaporation and stronger kinetic fractionation under warmer and drier boundary-layer conditions (Table 3). Both parameters increase with elevation, reaching values of 8.2 (slope) and 13.6 (intercept) at the highest-elevation site, RO (on average 1100 m a.s.l.). These values suggest limited secondary evaporation and a greater contribution of equilibrium-controlled condensation associated with orographic uplift. These observed differences are therefore primarily attributed to the elevation-dependant climatic conditions.

Mobile soil water (MW) closely tracked the LMWLs, with slopes ranging from 7.4 to 8.2 and intercepts from 2.0 to 15.8, suggesting minimal isotopic fractionation associated with evaporation or subsequent condensation processes. However, its isotopic range was narrower and on average more depleted than that of precipitation. Similar to precipitation, an elevation trend was also apparent in the isotopic composition of MW.

In contrast, TBW exhibited a reversed elevation trend both in slope (6.6–7.8) and in intercept (0.3–6.5) values decreasing with elevation. The differences became most pronounced at higher elevations. The differences between slopes were assessed using bootstrap resampling applied to the RMA regressions. At the high-elevation sites, most slope differences were significant, whereas the slopes of precipitation and MW at RO, and those of MW and TBW at LI, were not significantly different at the 95 % confidence level. In contrast, at lowland sites, no significant differences among the slopes of the regression lines were observed.

Figure 4 Seasonal comparison of the stable isotopic compositions of mobile and tightly bound soil water at two depths at selected study sites (the remaining two sites are shown in the Sect. S1, Fig. S1). Left panels: ZV (lowland agricultural field); right panels: LI (submontane spruce forest). Blue rectangles denote the daily precipitation (mm), and the light-blue dashed sine curve represents the weighted fit of the annual cycle of precipitation isotopic composition.

4.2 Comparison of mobile and tightly bound soil water

The differences in stable isotopic composition between MW and TBW were observed at all experimental sites. Among all components, D20 MW exhibited the greatest annual isotopic variability. At the ZV site, this variability reached 7.6 ‰ and 57.2 ‰ for δ¹⁸O and δ²H, respectively, but decreased with increasing elevation to only 5.6 ‰ and 43.2 ‰ at the LI site (Fig. 4). The largest contrast between D20 MW and TBW occurred in spring and autumn, with the maximum difference recorded on 7 March 2023 at the ZV site (3.3 ‰ and 20.7 ‰ for δ¹⁸O and δ²H, respectively).

A distinct phase shift between D20 MW and TBW was observed between February and May at all sites except RO (Fig. S1). The largest lag occurred at the ZV site, where the response of soil water to precipitation exceeded three months, representing the slowest response across the gradient. The lag decreased with elevation, shortening to approximately six weeks at LI site and becoming negligible at the highest-elevation site. The rapid response at the highest site likely reflects high annual precipitation (∼1400 mm), which frequently flushed the saturated soil profile. D40 MW and TBW exhibited broadly similar temporal dynamics and phase relationships to those observed for D20 TBW across all sites.

With the onset of the dry season and severe drought in July, soil profiles nearly dried out and isotopic differences among soil water pools diminished. With the decline of autumn period, isotopic amplitude peaks became synchronized, resembling the pattern of D20 MW but with attenuated signals. This attenuation likely resulted from mixing between newly infiltrated precipitation and residual water stored in the profile. Despite this homogenization, a temporal offset between precipitation and MW persisted at all sites except RO.

The influence of precipitation on soil water isotopic composition increased with elevation and precipitation input. At lower elevations (e.g., ZV) new precipitation gradually diluted pre-existing soil water, resulting in progressive isotopic enrichment. In contrast, at higher elevations, large precipitation inputs often replaced older soil water, resulting in soil water isotopic composition close to that of precipitation. Consequently, soil water in both horizons – particularly at the RO site (Fig. S1) – showed high variability and limited predictability.

Figure 5 Seasonal origin index (SOI) values across sites, depths, and soil water pools throughout the year 2023. Panels represent four study sites: ZV (lowland agricultural field), TD (mid-elevation meadow), LI (submontane spruce forest), and RO (montane beech forest). Soil water was sampled at two depths (D20 = 20 cm; D40 = 40 cm). MW and TBW denote mobile water and tightly bound soil water, respectively. SOI values near −1 indicate a dominant contribution from winter precipitation, whereas values approaching +1 reflect a dominant contribution from summer precipitation. Point colours represent SOI magnitude. The dashed black line illustrates the seasonal transition from winter- to summer-dominated precipitation inputs.

4.3 Origin of the soil water

Seasonal Origin Index (SOI; Fig. 5) values exhibited a consistent transition from winter-dominated signals in spring to summer-dominated signals later in the growing season across all sites. Site-specific differences reflected elevation-related gradients in precipitation and hydrological dynamics: higher-elevation forest sites (LI, RO) showed a stronger dominance of summer precipitation, whereas lower-elevation sites (ZV, TD) retained a more mixed seasonal signal. Moreover, the seasonal increase in SOI became steeper in mountainous areas, indicating a more rapid shift from winter- to summer-dominated precipitation and a shorter transitional period between these seasonal regimes. For the RO site, earlier seasonal data are not available (as explained above); however, the relatively high SOI values observed soon after snowmelt, together with the highest precipitation totals among the study sites, suggest that the seasonal turnover may be even more pronounced than at the LI site.

The differences between water pools were most evident at the lower-elevation sites (ZV and TD), where the transition from winter- to summer-dominated precipitation was delayed both between water pools at the same soil depth (e.g., July for MW vs. approximately August for TBW at D20 at the ZV site) and with increasing soil depth within the same water pool (August for D20 TBW vs. September for D40 TBW). Furthermore, D20 MW followed the annual course of precipitation most closely, whereas all other water pools exhibited a delayed response, particularly during the first half of the year, and retained an SOI signal from precipitation of the previous year. However, following a severe drought in July–August, the soil profiles were replenished by current precipitation. Consequently, the delay relative to precipitation was reduced, and differences among individual water pools largely disappeared, except for D20 MW.

4.4 Bulk soil water

Although BW was not directly sampled in this study, we present a conceptual illustration of its potential isotopic composition based on the obtained data. The results showed, that the stable isotopic composition of BW, estimated indirectly using a mass-balance mixing model (Eqs. 1 and 2), may vary statistically from TBW for both δ¹⁸O and δ²H. The unpaired t test (P<0.05) revealed this difference on at least one sampling date at both lowland study site (ZV, TD), with no difference observed at higher elevations. During the summer drought period, however, when soil desiccation removed almost all MW from the soil profile, BW effectively represented solely TBW.

Since the isotopic signature of BW depends both on the relative proportions of MW and TBW, and the isotopic contrast between them, the greatest deviations were observed during the spring and autumn seasons. During these periods, both the isotopic differences among water pools and the proportion of MW in the soil reached their annual maxima. At the ZV site, the discrepancy between BW and TBW was primarily driven by the pronounced isotopic contrast between mobile and tightly bound fractions, despite the low proportion of MW in the profile. In contrast, at the TD site, the difference was mainly attributed to a higher proportion of MW, while the isotopic contrast among the components was less pronounced.

For different BW and TBW, the average isotopic offset was 0.4 ‰ for δ¹⁸O and 2.3 ‰ for δ²H, with maximum differences reaching 1.7 ‰ and 5.2 ‰, respectively. When BW values were used to calculate the SOI, only minimal or no differences between BW and TBW were observed from a broader perspective. Notable deviations, however, occurred during individual sampling campaigns in April and May. In these instances, the SOI differences between BW and TBW reached up to 0.3.

5 Discussion

5.1 Isotopic changes due to soil properties and precipitation amount

All soil samples obtained in this study fell close to or directly on the LMWL. In consistency with previous studies (e.g., Goldsmith et al., 2012; Hervé-Fernández et al., 2016; Oerter and Bowen, 2017; Sprenger et al., 2018), MW was closely aligned with the LMWL, whereas TBW exhibited a lower slope and greater variability with some samples even outside the 95 % precipitation prediction interval. This reflects its longer residence time in the soil profile and indicates the influence of isotopic fractionation associated with evaporation and subsequent condensation and internal mixing processes within the soil matrix (Goldsmith et al., 2012; Sprenger et al., 2016).

In agreement with previous studies (e.g., Goldsmith et al., 2012; Geris et al., 2015; Hervé-Fernández et al., 2016; Sprenger et al., 2018), we observed distinct isotopic differences between MW and TBW, particularly in the upper part of the soil profile. These differences decreased with increasing depth, most likely reflecting longer residence times and enhanced mixing within deeper soil layers. Although measurements of cation exchange capacity (CEC) were not available, the very low clay content at all four study sites suggests that isotope fractionation associated with high clay-related CEC (Araguás-Araguás et al., 1995; Meißner et al., 2014; Oerter et al., 2014) was negligible. The observed differences in isotopic composition are therefore more likely attributable to differences in water retention and transport processes between macropores and capillary pore domains.

Despite the occurrence of extreme drought during the sampling year, which should leave an isotopically enriched signal in soil water, no such enrichment was observed at any of the study sites, regardless of precipitation regime or land cover. The reason can be most likely attributed to the sampling depth (20 and 40 cm), as isotopic enrichment from evaporation typically occurs at shallower depths, between 5 and 15 cm (Barnes and Allison, 1988; Sprenger et al., 2017; Oerter and Bowen, 2019; Dubbert et al., 2019). However, while Floriancic et al. (2024) reported no evaporative effect even at 10 cm depth, while other studies (Brooks et al., 2010; Sprenger et al., 2016) observed significant evaporative enrichment down to 30 cm. This discrepancy may be caused by the differences in soil texture or extraction methodologies and their associated, often unquantified, errors, particularly under low soil moisture conditions during drought periods (Sprenger et al, 2015; Orlowski et al., 2018).

In agreement with Kleine et al. (2020), we observed a greater phase shift of individual isotopic data in non-forested areas. This phase shift also increased with soil depth, particularly for MW. In contrast, the phase shift observed in TBW remained similar between shallow and deep layers. The greater lag observed in non-forested areas is likely driven by two main factors:

• Precipitation amount (Hervé-Fernández et al., 2016) for which higher rainfall can enhance leaching, thereby diminishing the isotopic distinction between MW and TBW.

• Vegetation cover as both soil texture and vegetation significantly influence the velocity of the wetting front (Xue et al., 2024).

Preferential flow pathways promote deeper and more rapid infiltration in forested areas, whereas under bare soil or grass, water infiltration proceeds more slowly and diffusively.

The unexpectedly rapid turnover in isotopic composition at our highest-elevation site (RO) contrasts with the results from other studies. For example, Floriancic et al. (2024) reported significant differences in soil water isotopic composition within the top 40 cm, even at forested sites with vegetation cover and precipitation amounts similar to those at our highest-elevation site. These discrepancies may be explained by differences in elevation, mean annual temperature, and soil texture. Lower elevations combined with higher temperatures contribute to a prolongation of the vegetation growing season, thereby increasing interception and evapotranspiration and reducing the infiltration of precipitation into the soil profile. In addition, the slightly higher silt content at their site likely enhances capillary water retention. Such capillary pores can hold water more effectively and may be bypassed by preferential flow paths, in contrast to the coarse sandy soils at our study sites. This comparison suggests that vegetation cover and soil properties may exert a stronger influence on soil water dynamics than precipitation amount alone.

5.2 Bulk soil water

Bulk soil water is commonly used as a proxy for the immobile fraction of soil water and is frequently compared with xylem water (e.g. Oliveira et al., 2025; Floriancic et al., 2024; Brighenti et al., 2024; Benettin et al., 2024; Barbeta et al., 2019, 2020; Goldsmith et al., 2019; Dubbert et al., 2019; Sprenger et al., 2016). However, the present results demonstrate a possible substantial difference between BW and TBW during certain periods of the year. Assuming that vegetation may access a less mobile or more tightly bound soil water pool partially disconnected from the mobile water contributing to groundwater recharge and streamflow (Brooks et al., 2010; Goldsmith et al., 2012; McDonell, 2014; Evaristo et al., 2015, 2019), then such differences may introduce uncertainties when identifying the sources of plant xylem water or may lead to an apparent mismatch between BW and xylem water isotopic composition (Oerter and Bowen, 2017; Bowling et al., 2017; Vargas et al., 2017; Barbeta et al., 2019; Lehmann et al., 2025).

Figure 6 Seasonal origin index (SOI) values of mobile soil water at the TD site at three available depths during 2022–2023. SOI values near −1 indicate a dominant contribution from winter precipitation, whereas values approaching +1 reflect a dominant contribution from summer precipitation. Point colours represent SOI magnitude. The horizontal dashed gray line illustrates the seasonal transition from winter- to summer-dominated precipitation inputs. Vertical dashed lines indicate the timing of maximum dominance of summer (red) and winter (blue) precipitation and their transition (black), highlighting the temporal lag of these events with increasing soil depth.

Although the differences observed in our study were most pronounced at the beginning of the growing season, historical data on MW indicate that these differences may vary over time. Some sites exhibit consistent behavior year by year – for example, the site at the highest elevation (RO), where no differences were observed between 2021 and 2023. Here, high total precipitation resulted in an isotopic composition of soil water that closely mirrored that of precipitation each year. In contrast, at lowland sites (e.g., TD; Fig. 6), MW exhibited a clear phase shift relative to precipitation throughout most of 2022, a year without extreme drought conditions. During the following year, when sampling including TBW was conducted, an extreme summer drought occurred, attenuating this phase shift. These observations suggest that differences between MW and TBW – and consequently between BW and TBW – may occur not only at the beginning of the year but also at other time intervals in dependence on dry or precipitation-rich conditions.

Although substituting TBW for BW may result in differences in SOI values of up to approximately 0.3, these differences do not appear to substantially affect the interpretation of long-term isotopic trends across years or study locations. However, this substitution should be considered with caution in short-term or single-event experiments (Lehmann et al., 2025; Muhic et al., 2024; Goldsmith et al., 2012, 2019), particularly under soil and environmental conditions indicating a higher proportion of MW. The high variability of MW can significantly dilute the BW signal and thereby mask the TBW-derived isotopic response.

5.3 Tightly bound water extraction

To obtain TBW, first the mobile fraction had to be removed. There are several studies proposing the methods for soil water extraction held in the soil matrix at different tensions (Geris et al., 2015; Bowers et al., 2020; Orlowski and Breuer, 2020). The results of these studies show different isotopic compositions of individual water pools, both with laboratory-prepared (Orlowski and Breuer, 2020; Bowers et al., 2020) and real soil samples (Geris et al., 2015). In this study, we use the pressure plate apparatus, similar to Orlowski and Breuer (2020), but using a different procedure. In their study, a spike experiment was performed, after which a pressure of 15 bar was applied, and the outflowing water was collected for isotopic analysis. Although labelled water was recovered during a specific time window, the initial and final stages of the experiment yielded water with isotopic signatures differing from the input. This method exhibits two basic limitations hindering its applicability to natural soil samples:

• The true isotopic composition of soil water is typically unknown, making it difficult to determine whether the observed isotopic composition already corresponds to soil water. This ambiguity arises from mixing between the soil water and the water used to saturate the ceramic plates within the apparatus, making the collected outflow likely a composite of both sources.

• To validate the method, every ceramic plate would have to be conditioned exclusively with samples from a single location and soil depth, to prevent internal mixing of different water sources. This requirement significantly reduces the practicality and scalability of the approach.

In our study, the pressure plate apparatus was employed in a modified configuration. Collected soil samples were subjected to a pressure of 0.6 bar, corresponding to the operational threshold for mobile water typically targeted by field-based suction lysimeters (Brooks et al., 2010; Muñoz-Villers and McDonnell, 2012; Berry et al., 2017; Sprenger et al., 2018; Haagsma et al., 2024). Unlike previous approach that rely on collecting the outflow water (Orlowski and Breuer, 2020), the presented method involves a subsequent extraction from pre-dried soil samples. This modification enables the simultaneous processing of up to 24 samples from various depths and locations within a single run. The subsequent extraction step is conceptually based on isotope mass-balance principles commonly used in hydrology (Haig et al., 2020; Zhao and Wang, 2021; Qiu et al., 2025) or isotopic modelling in general (e.g., Haagsma et al., 2024) and is applied here as an integral part of the extraction procedure itself to reconstruct the isotopic composition of the targeted soil-water fraction.

However, to obtain TBW, we recommend not to use the procedure described in the present study. Instead, the following approach can be recommended:

• Extract MW using standard suction lysimeters,

• Extract BW using each laboratory's standardized procedures (e.g., CVE, DVE-LS),

• Apply the mixing equation from this study only for the calculation of TBW, based on the measured isotopic composition of MW and BW and their relative absolute proportions determined gravimetrically (i.e., using a pressure plate apparatus).

This approach allows for the calculation of TBW while removing one procedural step used in this study, thereby reducing error propagation and the uncertainty in the final isotopic composition of TBW.

5.4 Data correction

Numerous studies have attempted to compare soil water (including both BW and MW) with xylem water (e.g., Zapater et al., 2011; Meunier et al., 2017; Vargas et al., 2017; Barbeta et al., 2019, 2020; Liu et al., 2021; Lehmann et al., 2025). To enable a meaningful comparison between soil water and xylem water, it is essential to employ an extraction technique minimizing isotopic alteration of the sample. However, it is well known that no currently available extraction method can extract soil water from all soil types and moisture contents without introducing some degree of isotopic bias (Sprenger et al., 2015; Orlowski et al., 2018; Kocum et al., 2025).

For this reason, various corrections are often applied to the measured data, although not universally. These include, for instance, adjustments to account for the presence of organic compounds (Martín-Gómez et al., 2015), or corrections based on Rayleigh-type fractionation models (Araguás-Araguás et al., 1995).

Another important but frequently overlooked limitation relates to the interpretation of method-validation experiments. Newly developed extraction techniques ate typically tested using soils of different textures and moisture contents (Dalton, 1988; Revesz and Woods, 1990; Leaney et al., 1993; Scrimgeour, 1995; Wassenaar et al., 2008; Kocum et al., 2025), and the results commonly reveal method-specific offsets expressed as shifts (± SD) relative to the isotopic composition of labelled reference water. Although these deviations are useful for method comparison, they are rarely incorporated as corrections into subsequently measured samples. To our knowledge, only Yang et al. (2023) explicitly addressed this issue by attempting to account for such method-specific offsets in the interpretation of extracted soil water isotope data.

The correction applied in this study relies on conducting spike experiments using the chosen extraction method with the site-specific soil types and varying water contents, in our case different dilution rates. This approach allows for the assessment of method performance for specific soil types collected at our study sites. Since soil texture and moisture content significantly affect extraction efficiency (Hendry et al., 2015; Orlowski et al., 2016), such validation experiments should be conducted separately for each soil type, ideally across a range of moisture conditions in each study. The resulting isotopic deviations (relative to labelled water) should then be incorporated into data correction procedures for actual samples. This condition seems to be essential for enabling meaningful comparison between individual water samples and between laboratories as well, as each CVE setup can differ and yield variable results (Orlowski et al., 2018).

Ultimately, the use of different extraction methods across laboratories, each associated with varying degrees of systematic errors, does not necessarily constitute a critical limitation. The results with low standard deviations, even in the presence of systematic offsets, can be quantified and subsequently corrected. Such calibration should allow meaningful comparisons across studies and research groups.

6 Conclusions

This study demonstrated that soil water stable isotope dynamics vary systematically along an elevational gradient. Soil water at lower-elevation sites with sparse or no snow cover exhibited longer residence times, whereas high-elevation sites with substantial winter snow accumulation showed more rapid isotopic turnover. All soil water samples plotted close to their respective local meteoric water lines, indicated minimal isotopic bias introduced by the applied extraction methods. Mobile soil water most closely mirrored the isotopic composition of precipitation, while tightly bound water, extracted using the newly presented mixing method, reflected its longer residence time in the soil profile and indicated the influence of isotopic fractionation associated with evaporation, subsequent condensation, and internal mixing processes within the soil matrix. The largest contrasts between mobile and tightly bound soil water pools occurred during spring and autumn, with maximum differences of 3.3 ‰ and 20.7 ‰ for δ¹⁸O and δ²H, respectively. In such cases, the difference between tightly bound and theoretically calculated bulk soil water reached up to 1.7 ‰ and 5.2 ‰, respectively. These findings highlight the importance of accounting for such variability, especially in short-term or single-time-point studies comparing soil and xylem water for plant source attribution. For such cases, a procedure was proposed to obtain tightly bound soil water. Future research should further explore how these dynamics interact with vegetation type, rooting depth, and changes in precipitation regimes under ongoing climate change.