Mechanisms of consistently disjunct soil water pools over (pore) space and time

The storage and release of water in soils is critical for sustaining plant transpiration and groundwater recharge. However, how much subsurface mixing of water occurs, and how much of the water is available for plants or otherwise percolates to streams and the groundwater is not yet understood. Based on stable isotope (²H and ¹⁸O) data, some studies have found that water infiltrating into soils can bypass older pore water. However, the mechanisms leading to the separation of water routed to the streams and water held tightly in smaller pores are still unclear. Here, we address the current limitations of the understanding of subsurface mixing and their consequences regarding the application of stable isotopes in ecohydrological studies. We present an extensive data set, for which we sampled the isotopic composition of mobile and bulk soil water in parallel with groundwater at a fortnightly temporal resolution and stream water and rainfall at a much higher resolution in a Mediterranean long-term research catchment, in Vallcebre, Spain. The data reveal that the mobile and tightly bound water of a silty loam soil in a Scots pine forest do not mix well; however, they constitute two disjunct subsurface water pools with little exchange, despite intense rainfall events leading to high soil wetness. We show that the isotopic compartmentalization results from the rewetting of small soil pores by isotopically depleted winter/spring rain. Thus, stable isotopes, and, in turn, water residence times, do not only vary across soil depth, but also across soil pores. Our findings have important implications for stable isotope applications in ecohydrological studies assessing the water uptake by plants or the process realism of hydrological models, as the observed processes are currently rarely implemented in the simulation of water partitioning into evapotranspiration and recharge in the critical zone.

1 Introduction

Rainfall infiltrating into soils often does not replace already stored water in the soil pores as described by translatory flow (Hewlett and Hibbert, 1967; Fig. S1 in the Supplement of this paper), and it can bypass large amounts of water held in pores prior to infiltration (Beven and Germann, 1982). It is important to understand which of these two processes is dominating the subsurface water flow due to their impact on plant water availability, nutrient or contamination transport, and groundwater recharge. Bypass flow has been studied for decades using artificial tracers and hydrometric methods (Beven and Germann, 2013), but its relevance for the partitioning between transpiration and groundwater recharge from the plot scale to the global scale has also been of interest over the last 10 years based on stable isotope (²H and ¹⁸O) data (Brooks et al., 2010; Good et al., 2015). Several studies have found different isotopic compositions between mobile soil water, which predominantly contributes to groundwater recharge and stream runoff, and water that is more tightly bound to the soil matrix (Brooks et al., 2010; Goldsmith et al., 2012; Hervé-Fernández et al., 2016; Gierke et al., 2016; Sprenger et al., 2018b). Following Sprenger et al. (2018b), here we define mobile water as soil water sampled with suction lysimeters, bulk water as soil water sampled with cryogenic extraction, and tightly bound water as the difference between mobile and bulk soil water, as it is the water in the pore space fraction that is not accessible by suction lysimeter. It was hypothesized by Brooks et al. (2010) that rainwater refilling the dry soil at the end of a dry season would lead to a distinct signal between mobile and tightly bound water. However, the mechanisms leading to a separation between water routed to the streams and water held tightly in smaller pores are unclear and currently under debate (Sprenger et al., 2016b; Berry et al., 2017; Dubbert et al., 2019). Nevertheless, understanding the stable isotopic composition and its variation in space and time is essential in the application of stable isotopes of water for paleoclimatological tree ring analyses (Gessler et al., 2014), plant water uptake assessments (Dubbert and Werner, 2019), partitioning between soil evaporation and transpiration (Wang et al., 2010), testing atmosphere–land surface models (Haese et al., 2013), understanding runoff processes (Klaus and McDonnell, 2013), and water residence time estimates (Sprenger et al., 2019).

So far, the experimental setups of studies addressing the disconnect between water flow across different pore spaces have had limited sampling frequencies and sample numbers (maximum twice per year in Brooks et al., 2010; Goldsmith et al., 2012, and five times per year in Hervé-Fernández et al., 2016), which has impeded progress in understanding the mechanisms explaining the disjunct subsurface water pools. Consequently, here we address the following research questions:

1. How do the isotopic compositions of mobile and tightly bound water vary over time and different hydrometeorological conditions?

2. Can we identify the source of the water more tightly bound to the soil pores?

3. What do the isotopic differences in subsurface water pools tell us about the hydrological processes, such as pore scale variability in water transport and age?

Our objective was to identify the dynamics of subsurface transport and flow by combining fortnightly isotope sampling of groundwater, and mobile and bulk (encompassing both mobile and tightly bound water) soil water during various hydrometeorological conditions (May to December 2015) with long-term rainfall and runoff isotope sampling and soil moisture and groundwater measurements (2011–2015).

2 Methods

2.1 Sampling site

Our study site is in the Can Vila catchment (0.56 km²) in the Vallcebre research area in the southeastern part of the Pyrenees (Spain; 42^∘11^′43^′′ N, 1^∘49^′13^′′ E) at 1200 m a.s.l. (Fig. 1). The climate is characterized as humid Mediterranean with an average precipitation of 880 mm yr⁻¹ that is distributed over about 90 d yr⁻¹ (Latron et al., 2009). The least precipitation occurs during winter and the most rain falls in May, October, and November, although the most intense rain events occur during summer. Snow plays a minor role with only <5 % of precipitation falling as snow. The average air temperature is 9.1 ^∘C and varies seasonally. The potential evaporation is about 823 mm yr⁻¹ with a strong seasonal dynamic ranging from about 20 mm per month in winter to up to 150 mm per month in summer (Llorens et al., 2018). The dry winters and the negative water balance during summer result in a succession of wet and dry periods with soils usually wetting up in early spring and early fall (Gallart et al., 2002). The soils studied are of silty loam texture (4 % sand, 72 % silt, and 24 % clay) with about 3 % gravel content and 6 % organic matter content (Molina et al., 2019). The soils' bulk density is 0.99 g cm⁻³, and the porosity is 62 % (Molina et al., 2019). The study area was terraced before and during the 19th century and abandoned during the second half of the 20th century. The slope at the study site is about 15^∘. The vegetation is Scots pine (Pinus sylvestris). The Scots pine trees are about 17 m high, the basal area is 45 m² ha⁻¹, and their canopy cover is 70 % (Molina et al., 2019).

Figure 1 Locations of the hydrometric measurements and the isotope sampling of soil water and groundwater (red circle), precipitation (red square), and stream water (red triangle) within the Can Vila catchment. The elevation is shown using contour lines, and the colors represent land use classes. The inset in the top left indicates the location within Spain (light grey) and Catalonia (dark grey).

2.2 Stable isotope data

We sampled the mobile soil water at depths of 20, 50, and 100 cm with suction lysimeters between May and December 2015 at approximately a fortnightly temporal resolution. The lysimeters consisted of 15 cm long porous cups (RSK ADAS Ltd., UK) with two tubes inserted that allowed for the creation of the vacuum in the lysimeter and for the sampling of soil water by injecting air into the lysimeter. We applied a suction of −700 hPa and sampled the water extracted within a few hours. In parallel, we took soil samples at 10, 20, 30, 50, 100 cm for bulk water cryogenic extraction (temperature: 110–120 ^∘C; duration: 120 min; vacuum: 10⁻² mbar; Martín-Gómez et al., 2015). We further sampled water from a piezometer that reached down to 263 cm and was screened along the bottom 60 cm. We took rainwater samples on an event basis from January 2011 to May 2013 and again from May 2015 until December 2015. We also sampled the stream at the outlet of the Can Vila catchment at least every 12 h, and at higher frequency during high flows, using a discharge-triggered autosampler throughout 2015 (Cayuela et al., 2018b). We retrieved monthly long-term rainfall stable isotope data from the Global Network of Isotopes in Precipitation (GNIP; IAEA/WMO, 2019) database for the Girona station (longitude: 41.91, latitude: 2.76; 129 m a.s.l., ca. 100 km southeast of Vallcebre). For the Vallcebre rainfall isotopes, we calculated the weighted monthly averages and weighted averages over 30 d prior to each day of soil sampling. We used the GNIP data to reveal long-term relationships between soil moisture and rainfall stable isotope compositions and to fill the gap in the Vallcebre data set (August 2013 to May 2015) via linear regression on monthly data: δ²H in Vallcebre = 1.1512 δ²H in Girona −5.1584 ‰  (r=0.79; p<0.01). We sampled throughfall and stemflow on an event basis between May 2015 and May 2016, but we do not present this information here. The isotopic composition of throughfall and stemflow was generally more enriched than that of rainfall, and all samples fell along the local meteoric water line (LMWL; Cayuela et al., 2018a). All water samples were analyzed for their stable isotopic composition (²H and ¹⁸O) using cavity ring-down spectroscopy (Picarro L2120-i, Picarro Inc., USA). The precision of the measurements is <0.1 for δ¹⁸O and <0.4 ‰ for δ²H. All isotope data are expressed in the δ-notation as parts per mill (‰) relative to Vienna Standard Mean Ocean Water (Coplen, 2011).

As cryogenic extraction has been found to be reliable for clayey loams with soil moisture levels > 30 % (Orlowski et al., 2016), we can exclude methodological issues with the isotopic analysis for our silty loam soils that had moisture contents > 30 % on the sampling days (Fig. 2d). Furthermore, we do not see a fixed offset between the isotopic compositions of mobile and bulk water, nor does the difference between the mobile and bulk soil water correlate with soil moisture. Thus, there are no signs of a bias introduced by the methods applied.

Figure 2 Temporal dynamics of hydrological conditions and stable isotopes during the soil water sampling period in 2015. (a) Catchment runoff at outlet (cyan), groundwater (GW) depth (grey) and daily rainfall volumes (black) at the study site. Note the logarithmic scale for discharge. (b) Stable isotope composition (δ²H) of stream water, groundwater, and rainfall. Blue and red shaded areas show the range of mobile water (MW) and bulk soil water (BW). Note that the size of the black squares indicates the rainfall amount, and the black lines represent the time span over which rainwater was cumulatively sampled. (c) Stable isotope composition of mobile water and bulk soil water at different depths. (d) Volumetric soil moisture at three depth intervals.

2.3 Hydrometric data

We measured the volumetric soil moisture at depths of 0–30, 30–60, and 60–90 cm using 30 cm long TDR probes (CS605, Campbell Scientific) at 20 min intervals from January 2013 onward (Molina et al., 2019) and computed daily and monthly averages. We continuously gauged the depth to the groundwater level using a Mini-Diver (Van Essen Instruments B.V.) compensated for barometric pressure in a PVC tube reaching down to 263 cm and screened along the bottom 60 cm. Catchment runoff was monitored at a 90^∘ V-notch weir using a water pressure sensor and a rating curve to get discharge values. Rainfall was measured in an open field adjacent to the study site using a 0.2 mm tipping bucket rain gauge (AW-P, Institut Analitic, Spain).

2.4 Data analysis

We used a non-parametric Spearman rank correlation to describe relationships between two variables, because the isotope data were not normally distributed according to the Shapiro–Wilk test. We assessed if the isotopic compositions of water samples were significantly different from long-term weighted mean values for precipitation using the Wilcoxon signed-rank test. Comparisons among piezometer, bulk and mobile soil water isotopic compositions were carried out with the Kruskal–Wallis test and a subsequent post-hoc Dunn test. The significance level for all statistical tests was set to the 95 % confidence interval.

We calculated the lc-excess as defined by Landwehr and Coplen (2006) for each water sample as follows: lc-excess = δ²H − a⋅ δ¹⁸O − b, with a being the slope of the LMWL (a=7.96) and b being the intercept of the LWML (b=12.89 ‰) (Casellas et al., 2019). A negative lc-excess value indicates kinetic fractionation of the isotopic composition of a water sample (Craig, 1961), which is depicted by plotting to the right of the LMWL in a dual isotope plot.

We estimated the isotopic composition (both ²H and ¹⁸O) of tightly bound water (δ_TW) based on an isotope mass balance approach. First, we derived the water retention curve, θ(h), using a pedotransfer function based on the soil texture and bulk density information (Schaap et al., 2001) (see the black line in Fig. S2). Afterwards, based on θ(h) and the pressure limit of the suction lysimeter, h_MW ($=-\mathrm{700}$ hPa), we determined the maximum soil moisture of tightly bound water (θ_TW=29 %). The soil moisture of the mobile water fraction (θ_MW) is the difference between the measured (bulk) soil moisture θ_BW and θ_TW. (Note that θ_TW was always bigger than θ_MW.)

We derived δ_TW as a function of δ_BW and δ_MW, which were the bulk and mobile water stable isotopic compositions, respectively:

$$\begin{array}{}\text{(1)}& {\mathit{\delta }}_{\mathrm{TW}}={\displaystyle \frac{{\mathit{\theta }}_{\mathrm{BW}}\cdot {\mathit{\delta }}_{\mathrm{BW}}-{\mathit{\theta }}_{\mathrm{MW}}\cdot {\mathit{\delta }}_{\mathrm{MW}}}{{\mathit{\theta }}_{\mathrm{TW}}}}.\end{array}$$

We calculated δ_TW for the depths at which bulk and mobile water were sampled in parallel (20, 50, and 100 cm).

We used the abovementioned pedotransfer function and applied the concept of dual-permeability as suggested by Gerke and van Genuchten (1993) to derive θ(h) and the hydraulic conductivity curve, K(h), for the fast flow domain (soil fracture), the slow flow domain (soil matrix), and the bulk soil. We calculated the pressure head according to the measured soil moisture based on θ(h) for each soil sampling day. We further calculated the hydraulic conductivity for the long-term rainfall isotope data set based on observed monthly mean soil moisture (as both plotted in Fig. S2).

For each sampling time (t), we calculated the fraction of event water (f_e [–]), which we defined as the share of water that was newly infiltrated from rainfall since the last sampling time (t−1). These calculations were carried out for sampling time steps in which the total soil moisture volume (V_a [L] = θ [L³ L${}^{-\mathrm{3}}]\cdot z$, with z being the representative depth of soil moisture measurements; z=30 cm) increased ($V_t>V_{t-\mathrm{1}}$), and the soil water isotopic composition changed from δ_t−1 to δ_t into the direction of the rainfall event (δ_e). We estimated the share of event water (V_e) from all soil water (V_a) as follows:

$$\begin{array}{}\text{(2)}& f_{\mathrm{e}}={\displaystyle \frac{V_{\mathrm{e}}}{V_{\mathrm{a}}}}={\displaystyle \frac{{\mathit{\delta }}_t\cdot {\mathit{\theta }}_t-{\mathit{\delta }}_{t-\mathrm{1}}\cdot {\mathit{\theta }}_{t-\mathrm{1}}}{{\mathit{\delta }}_{\mathrm{e}}\cdot {\mathit{\theta }}_t}}.\end{array}$$

3 Results and discussion

3.1 Continuous isotopic separation between mobile and bulk soil water

Mobile and bulk soil waters were significantly distinct in their stable isotopic composition (²H and ¹⁸O) at all depths throughout the study period (Fig. 2c) (Kruskal–Wallis post-hoc Dunn test, p<0.02). We found that the mobile water was more enriched in heavy isotopes than the bulk soil water for each sampling campaign. At every sampled depth, bulk soil water δ²H was significantly more depleted than mobile water, and bulk soil water below 10 cm was also significantly more depleted than groundwater (Kruskal-Wallis post-hoc Dunn test, p<0.01). Mobile water δ²H was significantly more enriched than groundwater (Kruskal–Wallis post-hoc Dunn test, p<0.05) at all depths. With soil depth, we observed a depletion of heavy isotopes in soil water, a decrease in temporal variability (Table S1 in the Supplement), and less difference between mobile and bulk soil water (Fig. 3). Our groundwater samples had intermediate δ²H and δ¹⁸O values and generally little variation (Fig. 2b). Groundwater was not significantly different to the weighted mean δ²H values for spring (one-sample Wilcoxon signed-rank test, p value = 0.09377) and fall precipitation (p value = 0.7982) and to the weighted mean δ¹⁸O values of the annual precipitation (p value = 0.4144). Stream water base flow isotopic compositions were similar to groundwater in terms of values and variability during base flow conditions. However, during rain events, the runoff isotopic composition responded immediately according to the rainfall δ²H and δ¹⁸O input (Fig. 2b).

Figure 3 Depth profiles showing the variability of the stable isotopic composition across the soil profile over the sampling period (color coded) for bulk soil water (stars), mobile soil water (circles), and groundwater sampled from the piezometer (diamonds). Vertical lines represent seasonal and annual weighted averages of precipitation. Location of the piezometer samples according to the depth of the water table in the piezometer.

The general pattern of disjunct isotopic compositions in mobile and bulk soil water was persistent, despite variable hydrometeorological conditions during our sampling period: there was a dry period between mid-June and mid-July that caused catchment runoff to cease (cyan line in Fig. 2a) in addition to a strong depletion of the soil moisture (Fig. 2d). Within our sampling period, we also covered four intense rain events between the end of July and the end of October that resulted in short-term increases of the water table in the piezometer up to the ground surface (grey line in Fig. 2a). We further observed continuously high soil moisture volumes between August and October with little variation (Fig. 2d).

δ²H values of mobile soil water at a depth of 20 cm were significantly correlated with 30-day weighted averages in δ²H of the antecedent rainfall of each sampling day (ρ=0.67, p=0.02). Such a relationship did not exist for the 50 and 100 cm mobile water sampling depths, nor for any of the sampling depths of bulk soil water. Averages over shorter time spans (e.g., 7 or 14 d) for the weighted average input did not correlate with soil water stable isotope compositions. Notably, we observed a significant trend in the isotopic enrichment of bulk soil water during the growing season (May to September) at depths of 10, 20, and 30 cm ($\mathit{\rho }=\mathrm{0.68},$ ρ=0.94, and ρ=0.85, respectively; p<0.03). This relationship indicates a slow refill of small pore spaces with isotopically enriched summer precipitation at depths of 20 and 30 cm; however, at a depth of 10 cm, the enrichment is caused by evaporative fractionation, as lc-excess decreased over the growing season (ρ=0.72, p=0.02). We also saw an increasing trend in δ²H values for mobile water at 50 cm (ρ=0.65, p=0.04), but no trends for lc-excess or samples taken at 20 and 100 cm. Mobile water δ²H values at depths of 20 and 50 cm were similar to the long-term weighted average for summer rainfalls (Fig. 4). Towards a depth of 100 cm, the mobile water approached the isotopic compositions of the piezometer samples, although, they were significantly different with respect to their δ²H values (Wilcoxon rank-sum test, p<0.001). In contrast, bulk soil water isotope values ranged mainly between average winter and annual rainfall isotope values. We observed kinetic fractionation (samples plot to the right of the local meteoric water line, LMWL, in Fig. 4) in the bulk water samples at a soil depth of 10 cm, which is indicative of soil water evaporation. This deviation from the LMWL, described as the lc-excess, was significantly more intense for the bulk water at 10 cm compared with other bulk water samples at deeper depths (Kruskal–Wallis post-hoc Dunn test, p<0.05), and all mobile water, piezometer, and stream water samples, which all plotted along the LMWL.

Figure 4 Water isotopes (δ²H and δ¹⁸O) of mobile water (MW) and bulk soil water (BW) as well as piezometer water. We further show weighted seasonal and annual averages of rainfall stable isotopes (squares), the local meteoric water line (LMWL; δ²H = 7.96 δ¹⁸O + 12.89 ‰), the global meteoric water line (GMWL; δ²H = 8 δ¹⁸O + 10 ‰), and the mass-balance-derived tightly bound water (TW). The box and whisker plots represent the mean (black line), the interquartile (box), the range (whiskers), and the outliers (black dots) of all MW (blue), piezometer (grey), BW (red), and TW (green) water samples.

We conclude that there is little mixing of tightly bound water (here defined as water that cannot be accessed via suction lysimeter) with recently infiltrated rainfall or with mobile soil water. While it was recently shown that mobile and tightly bound water would be in exchange during wet conditions in a silty clay soil (Hervé-Fernández et al., 2016), our observation for a silty loam soil does not support this, as the isotopic distinction persisted throughout extreme wet conditions with soil pores close to saturation. Our data do not indicate any periodic well-mixed conditions across soil depths (as frequently assumed in hydrological catchment models; McDonnell, 2014) or across the pore space at individual depths (as frequently assumed in soil hydraulic models; Sprenger et al., 2018b). We only recognized a replacement of older water with newly infiltrating water for the shallow mobile water as suggested in the translatory flow concept (Hewlett and Hibbert, 1967). Furthermore, the limited variation in bulk water stable isotopes and the similarity of the isotopic composition to that of isotopically depleted winter rainfall, indicates that tightly bound water is composed of relatively old water (see also discussion in Sect. 3.3). To underline this, we calculated the isotopic composition of the tightly bound water based on a mass balance approach (Eq. 1). This water, which cannot be sampled with a suction lysimeter, had isotopic compositions ranging between the weighted averages of spring/fall and winter rainfall (green points in Fig. 4). Therefore, our data set reveals that a large part of the subsurface water storage does not – or at least to a very limited degree and thus on long timescales – take part in the water flux towards groundwater and streams.

3.2 Filling of small pores during rewetting

To understand how the distinct isotopic composition in mobile and bulk soil water samples evolves, we combined our 8-month field data with multi-year monitoring data. The long-term water balance (2011–2017) is positive (rainfall > potential evapotranspiration) for the months between October and May. However, this is not the period of highest soil wetness, because there is little precipitation from December to March (<40 mm per month on average) which results in a dry-down of soils (Fig. 5b). Losses via evapotranspiration clearly surpass rainfall inputs (>10 mm per month) between June and August, but rainfall in these months often occurs as intense precipitation events (Fig. 2a).

Figure 5 Long-term hydrological and stable isotope dynamics. (a) Daily rainfall (bars) and potential evapotranspiration (PET; red line). (b) Rainwater δ²H values and soil moisture dynamics. Note that the size of the black and grey squares indicates the rainfall amount, and the line represents the time span of the water sample. Black squares show rainwater sampled at the study site, and grey squares show monthly values derived from GNIP data. The blue and red shaded area shows the range of mobile (MW) and bulk (BW) soil water, respectively. The green shaded area indicates the average value  ± standard deviation estimated for the tightly bound water (TW).

We connected these hydrometeorological and soil moisture data with rainfall stable isotopic compositions measured on an event basis at the study site and monthly GNIP data. We found that, in the long run, the rainfall happens to be more depleted in heavy isotopes when soil moisture is low (Fig. 6a). This significant relationship between rainfall isotopic composition and soil moisture (ρ>0.41, p<0.05) results from the combination of dry soil during winter, which is rewetted during early spring, and the temperature dependency of the rainfall isotopic composition that leads to its seasonality, with isotopically depleted rainwater during winter and spring (Fig. 5b). Rainfall was generally more depleted in heavy isotopes when monthly rainfall exceeded evapotranspiration ($\mathit{\rho }=-\mathrm{0.5}$, p<0.05, Fig. 6b). Thus, isotopically depleted rainfall from winter and spring preferably infiltrates into deeper soil layers and is not recycled into the atmosphere, as is the case for isotopically enriched summer rainfall. As a result of isotopically depleted rainfall infiltrating into dry soils when atmospheric demands are low, smaller pores are preferentially filled with the isotopically depleted waters.

Figure 6 Relationships between wetness and rainfall stable isotopes and refilling of soil pores. (a) Monthly rainfall stable isotopes (GNIP data) are more depleted when soils are drier. (b) Monthly rainfall at the study site is more depleted when monthly rainfall volumes are higher than potential evapotranspiration (positive water balance). (c) Fraction of newly infiltrated water (derived from mass balance) is higher when soil moisture was low before the infiltration event. All relationships (ρ shows Spearman rank correlation coefficients) are significant (p<0.05) with the exception of Fit BW (dashed red line), for which p=0.1.

The fraction of newly infiltrated event water in soil pores based on an isotope mass balance (Eq. 2) showed a strong negative relationship with the soil moisture before an event for mobile water ($\mathit{\rho }=-\mathrm{0.86}$, p<0.05) and a weaker relationship with bulk soil water ($\mathit{\rho }=-\mathrm{0.62}$, p<0.1). Mobile water was almost fully replaced by event water at low soil moisture and was less affected when rain infiltrated into wetter soil (Fig. 6c). The fraction of newly infiltrated event water in soil pores was generally lower for the bulk water than for mobile water underlining the high impact of tightly bound water that did not mix well with infiltrating rainwater.

Thus, our findings highlight that the ecohydrological separation of mobile and tightly bound water, as defined by the isotope study of Brooks et al. (2010), will occur in all settings where a soil of relatively fine texture (high potential for tightly bound water) dries out occasionally and is then rewetted by an isotopic composition that is distinct from the rainfall isotopic composition during wet periods. However, if there is little or no isotopic variability in the precipitation, tightly bound and mobile water may not be hydraulically connected, although we would not be able to detect that separation via stable isotopes.

3.3 Tightly bound water is old

The combination of the abovementioned mechanisms of (i) mobile and bulk water being continuously disjunct and (ii) tightly bound water being refilled by isotopically depleted winter rainfall leads to a pronounced variability in water ages across the soil pore space. Water held in smaller soil pores of the clayey fraction appear to be considerably older than water in pores of the sandy soil texture within the same soil profile. Thus, the isotope data presented suggest that the subsurface flow at the study site cannot be conceptualized as a uniform flow with water being exchanged across the entire continuum from coarse to fine pores (e.g., as assumed in the Buckingham–Richards equation (van Genuchten, 1980), Fig. S1d). Instead, we conclude that concepts based on, e.g., the dual-porosity model as described by Gerke and van Genuchten (1993) (Fig. 7d) for soil physical modeling or the use of two conceptual parallel systems representing matrix and preferential flow paths (e.g., Stumpp et al., 2007) will provide a better representation of the pronounced variability in flow velocities across the pore space.

Figure 7 Conceptualization of pore space water stable isotope variability. (a) Seasonal variation of the isotopic composition in rainfall. Arrows indicate soil sampling days. (b) Isotope depth profile of the fast flow domain (mobile water) during winter (grey) and summer (orange). The arrows below the profile represent the hydraulic conductivity. (c) Same as (b), but for the slow flow domain (tightly bound water). (d) Water retention curve based on the dual-porosity model as proposed by Gerke and van Genuchten (1993) for the fast (blue) and slow (red) flow domain with h_MW indicating the threshold for suction of lysimeter sampling. (e) Isotopic compositions across the pore space at a specific depth as shown in (b) and (c) with blue and red shaded areas representing mobile and tightly bound water pore space, respectively. Note that (d) and (e) share the same x axis, representing the pressure head and pore diameter based on Schjonning (1992). Vertical arrows in (b), (c), and (e) represent hydraulic conductivity in the fast (blue) and slow (red) flow domain. (f) Isotopic composition of the discharge from the soil subdivided into fast (blue) and slow (red) flow contributions and the resulting total isotopic composition (black dashed line).

Figure 7 summarizes our field observations and their potential implications. Due to the nature of two disjunct pore spaces, a seasonally varying isotopic composition in the rainfall is not preserved across the soil profile (as in Sprenger et al., 2016a, and Fig. S1 in this paper), but the fast and slow flow domains contain waters of different isotopic compositions with higher seasonal variation in the mobile water (Fig. 7b) than in the tightly bound water (Fig. 7c). The different water retention characteristics (as derived in Sect. 2.4) for the fast and slow flow domains (e.g., Gerke and van Genuchten, 1993) are directly linked to the soil pore diameters (Schjonning, 1992). The different flow pattern and the lack of exchange between the mobile and tightly bound water result in distinct isotopic compositions across the pore space at the same soil depth (red and blue shaded areas in Fig. 7e).

The extent to which the isotopologues (i.e., ²H¹H¹⁶O, ${}^{\mathrm{1}}{\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$) of mobile and tightly bound water undergo gas-phase or dispersion exchanges is unclear (Oerter and Bowen, 2017), but for the soil studied in the Can Vila catchment, exchange appears to be very limited. One approach to this question could be testing the different hypotheses of mixing using a recently presented soil hydrologic isotope model that allows for the conceptualization of various intensities of isotopic exchange between a fast and a slow subsurface flow domain (i.e., mobile and tightly bound soil water) (Sprenger et al., 2018b). Thus, based on our field data, we currently cannot assess the actual age (months or years?) of the tightly bound water, but experimental approaches with deuterated (enriched in ²H) water could help assess how much exchange can occur between the mobile and tightly bound water (Evaristo et al., 2019).

The maximum water volume stored in the tightly bound water pool is about 290 mm if we consider the upper 100 cm studied (max θ_TW=29 %, Fig. S2). This represents about one-third of the annual average rainfall, but it cannot be considered to be in inactive storage and completely immobile, as it would be partly available for evaporation and transpiration (permanent wilting point is often assumed to be about −15 000 hPa) and percolates according to subsurface pressure differences (likely at lower conductivities as shown in Fig. S2).

While there is increasing acknowledgement that subsurface water is not well mixed, and the consequences of this with respect to hydrological modeling are being tested (Fenicia et al., 2010; McMillan et al., 2012; Hrachowitz et al., 2013; Knighton et al., 2017; Cain et al., 2019), we propose extending this notion from the variability of isotopic tracers (and thus residence times) over depth (Fig. 7b, c) to the variability of tracers (and residence times) over the pore space continuum (Fig. 7e). Due to the highly non-linear relationship between soil moisture, pressure heads and hydraulic conductivity, the tightly bound soil water will percolate several magnitudes slower than the mobile water (Fig. S2). The mechanisms outlined result in isotopically enriched rainwater infiltrating at higher soil hydraulic conductivities, which leads to a relatively rapid transport of young water to the groundwater and streams (Fig. 7f). When isotopically depleted rain infiltrates during dry conditions, the hydraulic conductivity is generally lower; furthermore, the differences in hydraulic conductivity between mobile and tightly bound water are also lower than during wet conditions, resulting in lowered percolation rates of isotopically depleted rainfall (squares in Fig. S2). Experiments with weighable lysimeters (e.g., Stumpp and Maloszewski, 2010; Benettin et al., 2019), where the stable isotopes of the outflow were sampled, reported similar results to our field observations for the mobile water: rapid percolation via preferential flow paths. However, the isotopic composition of the stored bulk soil water is unknown in such experiments, as destructive soil sampling is not possible in weighing lysimeter studies. Thus, to study interaction (or the lack of it) between tightly bound and mobile water, lysimeter experiments would need to be extended by in situ soil water vapor isotope measurements.

Our observations of the pore-scale flow velocity variability of disjunct subsurface water pools has consequences with respect to the contact time with mineral surfaces and will therefore affect nutrient concentrations (Asano et al., 2006), which will then be higher for the water in the slow flow domain. As hypothesized by Brantley et al. (2017), such nutrient concentration differences might be one reason why plants have been shown to take up tightly bound water even though highly mobile water is available (Brooks et al., 2010). The pore-scale water age differences outlined above could further explain why plant stem water stable isotopic compositions indicate that trees preferentially use winter precipitation (relatively old water) during summer, as recently shown by Allen et al. (2019). Moreover, the different origin (timing of infiltration and associated isotopic composition) of mobile and tightly bound water (as shown in Fig. 7e) can explain differences observed among different methods of soil water isotope analysis (Sprenger et al., 2015; Geris et al., 2015; Orlowski et al., 2016; Tsuruta et al., 2019). The distinct subsurface flow paths of a younger more mobile water domain and a domain of low flow velocities (consequently older soil water) has further implications for the interpretation of recharge estimates based on depth profiles of pore water stable isotopes (Koeniger et al., 2016) or tritium (Li et al., 2019). Moreover, residence time estimates based on isotope data sampled by suction lysimeter will overestimate the soil water turnover rates (Sprenger et al., 2018a).

4 Conclusions

The mechanism by which different isotopic compositions for the mobile and tightly bound soil water develop at the plot scale have not previously been presented. Our extensive field data show that mobile and tightly bound water of a silty loam soil in a Scots pine forest were continuously separated during various environmental conditions over our 8-month study period. These differences resulted in the preferential refilling of small soil pores by isotopically depleted rainfall during low soil moisture conditions. This water, which replenishes the smaller pores, is held more tightly and contributes to a lesser extent to the soil water flux; therefore, it is also older. Thus, the variability of the stable isotopic composition and soil water ages across the soil pore space continuum is contrary to the common assumptions made in stable isotope applications in both the hydrological modeling of water fluxes and ecohydrological studies that assess the root water uptake depth of trees.

In both applications of stable isotopes, it is usually presumed that the soil water is well mixed across the pore space, which means, for example, that the water within the sandy soil fraction has the same isotopic composition as the clayey soil fraction at the same depth in a soil profile, as the water exchanges over the entire pore space. Here we showed for our Mediterranean study site that this is not the case, and we explained the pore-scale isotopic variability using the synchronized water balance and rainfall stable isotope dynamics.

While the differences in stable isotopes between mobile and tightly bound soil water have been observed in other studies, the explanation for these differences have been limited. We suggest that the long-known soil physical processes in heterogeneously structured soils explain the variability of stable isotopes across the pore space continuum and that this will need to be acknowledged in ecohydrological field and modeling studies based on stable isotope data.