Modelling temporal variability of in situ soil water and vegetation isotopes reveals ecohydrological couplings in a riparian willow plot

The partitioning of water fluxes in the critical zone is of great interest due to the implications for understanding water cycling and quantifying water availability for various ecosystem services. We used the tracer-aided ecohydrological model EcH₂O-iso to use stable water isotopes to help evaluate water, energy, and biomass dynamics at an intensively monitored study plot under two willow trees, a riparian species, in Berlin, Germany. Importantly, we assessed the value of in situ soil and plant water isotope data in helping to quantify xylem water sources and transit times, with coupled estimates of the temporal dynamics and ages of soil and root uptake water. The willows showed high water use through evapotranspiration, with limited percolation of summer precipitation to deeper soil layers due to the dominance of shallow root uptake (>80 % in the upper 10 cm, 70 %–78 % transpiration/evapotranspiration). Lower evapotranspiration under grass (52 %–55 % transpiration/evapotranspiration) resulted in higher soil moisture storage, greater soil evaporation, and more percolation of soil water. Biomass allocation was predominantly foliage growth (57 % in grass and 78 % in willow). Shallow soil water age under grass was estimated to be similar to under willows (15–17 d). Considering potential xylem transit times showed a substantial improvement in the model's capability to simulate xylem isotopic composition and water ages and demonstrates the potential value of using in situ data to aid ecohydrological modelling. Root water uptake was predominately derived from summer precipitation events (56 %) and had an average age of 35 d, with xylem transport times taking at least 6.2–8.1 d. By evaluating isotope mass balances along with water partitioning, energy budgets, and biomass allocation, the EcH₂O-iso model proved a useful tool for assessing water cycling within the critical zone at high temporal resolution, particularly xylem water sources and transport, which are all necessary for short- and long-term assessment of water availability for plant growth.

1 Introduction

Understanding how water is partitioned in the critical zone (CZ), the near-surface zone from the top of the vegetation canopy to groundwater (Grant and Dietrich, 2017), is essential for improving knowledge of landscape functionality while providing an evidence base for sustainable water management strategies. The partitioning of CZ water is strongly dependent on evapotranspiration (ET) (which accounts for > 60 % of terrestrial precipitation) (Oki and Kanae, 2006; Zhang et al., 2016), with vegetation water use globally accounting for 65 %–70 % of evapotranspiration (Good et al., 2015; Schlesinger and Jasechko, 2014). While measurements of ET and transpiration fluxes help to quantify water partitioning, these measurements usually do not constrain the dynamics of how water is taken from different water sources (i.e. different soil depths), which may greatly change with wetness conditions (Rothfuss and Javaux, 2017), the climate zone (Amin et al., 2020), and seasonally (Barbeta and Peñuelas, 2017). In regions where ET dominates precipitation water partitioning (> 90 %, Zink et al., 2017), seasonal variations in ET greatly reduce water availability during the growing season. Additionally, the importance of evaluating the water footprint of biomass production in water-limited regions where ET is dominant highlights the importance of partitioning ET into its components of interception and soil evaporation and transpiration (Kool et al., 2014; Xiao et al., 2018). However, partitioning ET to quantify the transpiration component throughout the growing season is complicated by multiple factors, including atmospheric demand, vegetation and root conductance, rooting distribution, stomatal resistance, and water potential throughout the soil profile (Dubbert and Werner, 2019; Jones and Tardieu, 1998; Sperry and Love, 2015). While total transpiration water usage may be constrained by direct measurement at plot or stand scales (e.g. via sap flow and eddy covariance, Kool et al., 2014) the source of transpired water from the rooting zone throughout the growing season is not easily measured and remains highly uncertain (Brantley et al., 2017; Dubbert and Werner, 2019).

Continued efforts to close the knowledge gap on root uptake sources are essential, as increasing climatic variability and accentuated extremes will likely affect future crop and timber production (Lobell and Gourdji, 2012; Lobell et al., 2011) as well as groundwater and stream water availability (Gudmundsson et al., 2019; Taylor et al., 2013). Tracers in soils and xylem, specifically conservative water stable isotopes of deuterium (²H) and oxygen-18 (¹⁸O), have previously been shown to be effective tools to help constrain root water uptake sources using various approaches (Rothfuss and Javaux, 2017). These approaches usually use mixing relationships, including linear mixing models of water derived from different pools (e.g. Barbeta and Peñuelas, 2017), the similarity of potential source waters and xylem water (e.g. ellipsoid method) (Amin et al., 2020; Tetzlaff et al., 2021), Bayesian mixing frameworks (cf. von Freyberg et al., 2020), and physically based modelling approaches (Knighton et al., 2020; Ogle et al., 2014; Sutanto et al., 2012). However, many of these previous studies have utilized data obtained from destructive sampling, which limits the number of samples and may not account for the transport lag of soil to xylem sample locations (von Freyberg et al., 2020); diurnal variability (De Deurwaerder et al., 2020); or stem internal storage and exchange with xylem (Steppe et al., 2006). Moreover, increasing uncertainty surrounds whether common destructive methods, which involve cryogenic extraction, actually represent xylem water (Chen et al., 2020). Recent developments of in situ tracer measurements of soil and xylem (Marshall et al., 2020; Oerter and Bowen, 2017) have enhanced the possibilities for spatiotemporal evaluation of root uptake distributions, in particular when integrated with other data in physically based ecohydrological models. Higher-resolution sampling can be beneficial, as it can provide a wider range of temporal conditions (e.g. event rewetting), which helps to test model performance and in the assessment of model structure.

Physically based ecohydrological modelling approaches have a wide range of applicability, including the estimation of water and energy fluxes and storages, flux partitioning, and biomass accumulation (Asbjornsen et al., 2011; Fatichi et al., 2012; Maneta and Silverman, 2013). Recent developments linking tracers into ecohydrological models (e.g. Kuppel et al., 2018a) have further expanded this potential to track water flow paths and associated ages, which can aid in understanding water cycling and mixing within the CZ (Geris et al., 2017; Penna et al., 2018; Sprenger et al., 2019). Such modelling approaches can help overcome the limitations of using isotopic data alone regarding temporal frequency and spatial heterogeneity, which may affect estimates of water partitioning and source identification (Goldsmith et al., 2018; Rothfuss and Javaux, 2017; Sprenger and Allen, 2020). However, relatively few studies have utilized physically based models to estimate storage–flux–age dynamics while considering the mixing of water after root uptake (Knighton et al., 2020), particularly at high resolution, to account for sub-daily variability (e.g. De Deurwaerder et al., 2020) or with consideration of root length properties (Gessler et al., 2021; Seeger and Weiler, 2021).

Here, we utilized soil and xylem water isotope data from in situ monitoring over the 2020 growing season in Berlin, Germany, to help calibrate a tracer-aided ecohydrological model and estimate water flux and uptake dynamics. The in situ location, which, crucially, also monitored complementary soil moisture and vegetation growth dynamics, comprises a small stand of riparian willow trees and surrounding grassland in a situation typical of north-eastern Germany. The model, EcH₂O-iso, is a distributed physically based model that couples water and energy fluxes, with vegetation carbon allocation across the soil–plant–atmosphere continuum (SPAC). The primary goal of this study was to utilize the EcH₂O-iso model as a learning framework to evaluate the temporal linkages of energy, water, and the ecology of the plot site throughout the growing season. The primary research goal was evaluated through the following research questions. (1) How well can the tracer-aided ecohydrological model reproduce in situ measurements of fluxes, storages, and stable water isotopes of two juxtaposed vegetation types (willows and grass)? (2) Does distance-based isotopic mixing of root water uptake better approximate xylem water isotopes than uniform instant root uptake mixing? (3) What are the transit times of vegetation water, and what are the implications for water usage? We used multicriteria calibration of water, energy, biomass, and isotopic datasets as equally informative data within the EcH₂O-iso model, with hydrological, isotopic, and ecological model outputs utilized as tools to evaluate these research questions. Such an overall assessment of ecohydrological partitioning, in fluxes as well as root uptake, was intended to help to improve the conceptualization of water cycling in the CZ at high-temporal resolutions and contribute to an evidence base for management strategies in water-sensitive areas.

2 Study site and background

2.1 Study site

The study site is in a peri-urban area on the grounds of the Leibniz-Institute of Freshwater Ecology and Inland Fisheries (IGB) in south-eastern Berlin (Fig. 1b). The site is situated near Lake Müggelsee (80 m north of the lake edge), encircled in the north and west by buildings (40 m) and east and south (30 and 20 m, respectively) by additional vegetation. A stream draining nearby ponds fringes the east of the study area. A lake water extraction facility (using bank filtration) is situated immediately north of the IGB site, and groundwater is ∼2.2 m below the ground surface with limited annual variability (<0.1 m). A previous isotope-based study at the site has excluded groundwater and the nearby stream water as likely sources of water to the trees (Landgraf et al., 2021).

Figure 1 Study site location and site description. (a) Location of the site within Germany. (b) Location of IGB within Berlin (2021). (c) Plan view of the study site including the location of the willow trees, soil water isotope and moisture, and automatic weather station (AWS). (d) Conceptual diagram of the study site, with three soil depths and two primary vegetation units (with relative prominence). Pink and blue columns in (d) indicate soil moisture and isotope measurement columns, and grey box indicates AWS location. Berlin maps © OpenStreetMap contributors, 2021. Distributed under the Open Data Commons Open Database License (ODbL) v1.0. The area surrounding the study site is grass covered, similar to Site B.

2.2 Soils and vegetation

The study site was chosen based on its open nature (not shaded) for the two willow trees (Salix alba) and an adjacent automatic weather station (Fig. 1c and d). Willows are a riparian species which is well adapted to moderate and high moisture typical of riparian areas or those with high groundwater levels (Isebrands and Richardson, 2014). Additional site selection criteria are provided in Landgraf et al. (2021). The willows have a similar age of ∼15 years with the willow in the north (Willow 1, Fig. 1c and d) being slightly larger (∼9 m height and stem diameter of 398 mm, 16 July 2020) and the southern willow (Willow 2, Fig. 1c and d) being slightly smaller (∼8 m height and stem diameter of 353 mm, 16 July 2020). Sparse grass is present below the willow trees, with some bare-soil patches (Fig. 1d). Grass coverage is greater with no clear bare-soil patches in the open area surrounding the willows (Fig. 1d).

The site is situated in the North European Plain, with geology and surficial soils in the surrounding areas of Berlin deposited during the Weichselian glaciation (Deutschen Nationalkomitee, 2016). The willows are situated on reclaimed ground where a previous building was demolished and the site was backfilled with sandy brown earth topsoil. Soil samples taken to 1 m depth reveal a relatively uniform soil structure, with a slightly higher organic horizon in the near-surface soils where rooting is densest, which is more developed below the open grassland area south of the willows (Site B, Fig. 1c and d).

2.3 Climate

The climate is continental with a maritime influence (Köppen index Cfb) and experiences substantial interannual variability in precipitation (Table 1). The range of annual precipitation from within the climate normals (1980–2010) is 80 % of the long-term average precipitation (526 mm yr⁻¹, Table 1) (DWD, 2021). While precipitation during the study period was similar to the climatic normal, the site experienced lower-than-normal humidity and higher wind speeds and air temperature (Table 1). The growing season is less humid and is warmer than annual averages (Table 1), with large sub-daily variability in hydroclimate. For the vegetation growing season in the surrounding region as well as Berlin, evapotranspiration (ET) is the dominant hydrologic flux, accounting for ∼90 % of total precipitation (Gillefalk et al., 2021; Smith et al., 2020a).

Table 1 Climate conditions at a nearby long-term weather station, Berlin Brandenburg Airport (DWD, 2021). Average (standard deviation in round brackets) climate conditions are shown for the climate normals (1980–2010) and the study year (2020). The growing season is May–September (inclusive). n/a is not applicable.

3 Materials and methods

3.1 In situ measurements

The study site was set up to continuously measure hydroclimate, soil moisture, and vegetation productivity. Hydrometeorological conditions at the site were monitored with a mobile eddy covariance system (LI-COR Biosciences, Lincoln, NE, USA), which included measurements of precipitation amounts (event) and a collection of precipitation samples for isotope analysis, humidity, air temperature, wind speed, short-wave radiation, surface temperature, and latent heat (Site B, Fig. 1c). Climate data were quality-checked against nearby weather stations that have been in operation longer. These include measurements from the roof of the nearby IGB building (2013), the station on Lake Müggelsee (2013), and surrounding DWD (Deutscher Wetterdienst) weather stations (DWD, 2021). Precipitation samples were collected at 4-hourly intervals. Each sample bottle was filled with 0.5 cm of paraffin oil to prevent evaporation. Precipitation samples were analysed with a Picarro L2130-i cavity ring-down laser spectrometer (Picarro, Inc., Santa Clara, CA, USA). Hourly pictures were taken from the soil surface towards the canopy to capture changes in leaf coverage. Images were translated to an estimated leaf area index (LAI) using the LAI package for R (Martin, 2015), which estimates the gap fraction of the histogram-based unimodal threshold method.

Soil moisture, isotopes (δ²H and δ¹⁸O), and temperature were measured at two sites in the study area, Site A near Willow 2 and Site B in the open grass area. Soil moisture and temperature were measured at three depths (10, 40, and 100 cm) with water content reflectometers (CS616, Campbell Scientific, Inc., Logan, UT, USA; ±2.5 % for VMC, volumetric water content) and thermistors (BetaTherm 100K6A1IA, Campbell Scientific, Inc., Logan, UT, USA; ±0.2 ^∘C). Visual inspection of rooting densities during installation suggested a higher proportion of roots near the surface and very few tree roots at 1 m.

Soil water isotope samples were collected using both destructive bulk soil sampling and with an in situ vapour isotope analyser. Installation of polypropylene membranes (7 cm long, 0.2 µm pores, Kübert et al., 2020) for soil water vapour extraction at three depths (10, 40, and 100 cm) was conducted at the end of May 2020 (20 May 2020). Extracted in situ soil vapour samples were analysed every 2 h using a Picarro L2130-i cavity ring-down laser spectrometer (Picarro, Inc., Santa Clara, CA, USA), calibrated using the approach presented for this site in Landgraf et al. (2021).

To complement the climate and soil measurements, detailed measurements indicative of vegetation water fluxes and biomass accumulation were conducted continuously throughout the study period. Measurements included sap flow (heat ratio method, SFM1 Sap Flow Meter, ICT International, Australia), dendrometers (DR, Radius Dendrometer, Ecomatik, Dachau, Germany), and xylem water vapour (Picarro, Inc., Santa Clara, CA, USA). Polypropylene membranes (as per the soil) were installed in two bore holes in both willow trees starting in June (Willow 2). COVID-19-related complications delayed the installation in Willow 1 (August). Xylem water vapour was extracted from two bore holes (inner diameter of 8–10 mm, described in Marshall et al., 2020) on each willow tree, at 30 and 170 cm, and analysed in the same analyser as for the soil vapour isotopes every other hour. To prevent contamination of atmospheric air within the xylem vapour samples, the probe was embedded within the stem and sealed using commercial silicon. Analysed xylem water isotopes were evaluated for wounding effects (in early periods after installation), and any were removed from the analysis. Wounding effects were identified by daily inspection of the trees and unrealistically enriched isotopic compositions which exponentially decreased in time.

Insulated heated cables were installed with the tubing for all of the vapour to avoid condensation effects and modulate vapour temperatures. Additionally, tubing was flushed with dry air for 10 min each morning if condensation was detected to remove any residual water. To test and correct for water concentration and temperature dependencies, a linear regression of vapour concentration and slopes of ²H and ¹⁸O, as well as a nonlinear polynomial regression of daily average xylem and soil water isotopes with temperature (air and soil, respectively), was conducted to determine the strength of the relationship for water concentration and temperature-dependent offset. Corrected in situ xylem and soil water isotopes for water concentration and temperature were consistent with bulk soil water isotope and twig cryogenically extracted isotope samples (see Supplement and Landgraf et al., 2021, for more details on isotopic measurement).

3.2 EcH₂O-iso model

EcH₂O-iso is a physically based, distributed, tracer-aided model coupling vegetation, soil energy and water balance, and carbon utilization (Maneta and Silverman, 2013). The isotope module (simulating δ²H, δ¹⁸O, and water age) was coupled with the water balance to track the movement of water throughout the model domain (Kuppel et al., 2018a). The following section presents a brief synopsis of the model conceptualization of energy, water, and isotope balances (conceptual diagram, Fig. S1 in the Supplement). Further details of the model are provided in Maneta and Silverman (2013) and Kuppel et al. (2018a).

3.2.1 EcH₂O-iso energy balance

The energy balance is resolved with a top-down approach, conducting energy balance within the canopy before the energy balance at the surface. The canopy energy balance partitions incoming radiation into latent heat (interception evaporation and transpiration), sensible heat, and net radiation as a function of canopy temperature (Fig. S1). The canopy energy balance is very sensitive to the canopy stored water (maximum canopy storage parameter, CWS_max), where higher intercepted water storage decreases energy availability for transpiration. Transpiration is limited by environmental constraints, implemented using a Jarvis-type stomatal conductance model dependent on soil moisture, vapour pressure deficit, air temperature, and incoming radiation. EcH₂O-iso was modified to account for potential root uptake from outside the vegetation model cells. The radial rooting parameterization (Sect. 3.3.1) estimated the rooting proportion in each surrounding model cell which assumes a radial distance from the cell centre. Rooting proportions with depth (k_root, described in Kuppel et al., 2018) and laterally were utilized to account for total transpiration water availability and lateral proportions used to withdraw water from neighbouring cells.

The surface energy balance utilizes energy translated from the canopy energy balance to resolve latent heat (soil evaporation), sensible heat, net radiation, ground heat, and snow and melt heat (only if snow is present) using surface temperature. Ground heat is resolved using two thermal layers, where the depth of the first thermal layer is defined as half of the depth where the thermal wave is damped by 37 % (Arya, 2001). Interpolation of soil temperatures for each soil layer are estimated using surface and soil temperature at the bottom of the thermal layers and the linear damping equation (Arya, 2001) solved at the bottom of each soil layer. Soil temperature interpolation does not influence the energy balance.

3.2.2 EcH₂O-iso water balance

Similar to the energy balance, the water balance is estimated with a top-down approach including canopy, surface, and sub-surface storage (Fig. S1). Canopy storage is estimated using a linear bucket approach, with maximum storage limited by a storage parameter (CWS_max) and canopy interception driving drawdown of canopy storage. Net precipitation as throughfall accumulates on the soil surface (ponded water) and infiltrates into the shallow soil (layer 1) using the Green–Ampt model (parameters for Brooks–Corey, air entry pressure, and vertical hydraulic conductivity which are λ_BC, ψ_ae, and K_v, respectively). Ponded water at the end of each time step is directly routed to ponded water in the next downstream cell. Soil water redistribution is conducted using gravitational drainage when field capacity is exceeded. Redistribution is estimated for each model layer. The model assumes that there is no preferential flow in any soil layer. Very dry soils at the study site prompted the introduction of the sub-discretization of shallow soils (layer 1) to estimate the moisture at discrete depths in addition to the average moisture of the layer. Shallow soils were sub-discretized into 1 cm increments with incoming water (infiltration and return flow) entering from the layer boundaries and redistributed using gravitational drainage. Sub-discretization was implemented for informational purposes only, not for use in calibration. Gravitational drainage rates linearly increase from zero (at field capacity) to saturated vertical conductivity when the soil is fully saturated. Upward redistribution of water occurs if deeper soil storages are fully saturated. Vertical downward flux from the deepest soil layer (layer 3) can occur due to leakance out of the model domain. Lateral flow may occur in the deepest soil layer, with water above field capacity routed to the next downstream model cell using a linear kinematic wave model.

3.2.3 EcH₂O-iso isotope mixing and water ages

Mixing of isotopes (δ²H and δ¹⁸O) and water age is conducted using a complete mixing assumption (fully mixed at the end of each time step) in each soil storage (layers 1–3 independently) and in the canopy and ponded water. Storage mixing is conducted with amount-weighted averaging of isotopic compositions (or age) with incoming fluxes. Outfluxes have the same isotopic composition (and age) as the mixed water in storage. Evaporative fractionation is calculated for the shallow soils (layer 1) using estimated soil evaporation and the Craig–Gordon fractionation model (Craig and Gordon, 1965). Humidity in the soil is estimated using the method proposed by Lee and Pielke (1992), with the kinetic fractionation factor modified for use in soils (Braud et al., 2005). Isotopic composition and the age of transpired water are estimated by amount-weighting the contribution of root uptake from each soil layer (root uptake proportion dependent on water availability and rooting distribution). Amount-weighted root uptake assumes instant translocation, complete mixing, no fractionation, and uniform isotopic composition throughout the vegetation using soil water isotopes of the same time step. At the end of each time step the age of water in storage is advanced by one time step (e.g. 1 h).

To evaluate the distribution of water ages, parameter sets were re-run for each time step that included precipitation (925 h with precipitation over the study period). The input precipitation concentration was changed in a stepwise manner with a concentration of “1” on the precipitation event evaluated (Smith et al., 2020b). For example, for the 25th precipitation event (i.e. 25th hour with precipitation, here, 7 January at 11:00), the precipitation concentration was set to 1, while all other precipitation and initial storages had a concentration of 0. The concentration was tracked through the domain to provide a proportion of water in each storage originating from each precipitation event. The cumulative concentration of different precipitation events sums to 1 when all of the water in storage is younger than the duration of the simulation.

3.3 Tree water mixing and transit times

Knighton et al. (2020) showed that the improved simulation of xylem isotopes could result from considering potential vegetation storage and the mixing of soil water older than one time step within EcH₂O-iso. Similarly, Seeger and Weiler (2021) showed a promising conceptualization of isotope mixing within vegetation based on a convolution approach using in situ measurements of soil and xylem water. To evaluate the dynamics of transit times from the root to measurement height during the study period, a combination of modelled results of sap flow, soil water isotopes, and spatiotemporal root uptake proportions were used to estimate vegetation xylem water composition and age at the average measurement height (1 m) while accounting for spatial variability in soil water isotopes (laterally and with depth). Root uptake proportions were calibrated using parameterized exponential distributions for depth and lateral contributions in EcH₂O-iso. The tree water mixing routine (applied following EcH₂O-iso calibration) was developed with the assumption that isotopic mixing is dependent on the rooting distance (similar to Seeger and Weiler, 2021). Mixing of source water was conducted using the convolution equation and expanded on Seeger and Weiler (2021) to include lateral soil heterogeneity (e.g. outside model cell or different soil layers). Contrary to Seeger and Weiler (2021), rooting distributions are not parameterized through mixing and are calibrated through energy and water balance in EcH₂O-iso. This approach differs from tree mixing within EcH₂O-iso which instantly mixes uptake water throughout the tree and from previous tree mixing approaches in EcH₂O-iso (Knighton et al., 2020) by specifically exploring the capabilities of rooting distributions and uptake travel times to describe xylem water isotopes. Hereafter, the new xylem mixing will be referred to as distance-based mixing, and the method used by EcH₂O-iso is referred to as instantaneous mixing.

3.3.1 Vertical and radial rooting length

The average vertical and horizontal distances (and total distance) of roots in each layer to the measurement height were estimated and calibrated in EcH₂O-iso with a modification of the root depth and radial spread approach used in Sperry et al. (2016). The parameterization of this approach is in synchrony with the rooting distribution within EcH₂O-iso. Following EcH₂O-iso, the rooting distribution (k_root) was used to determine the proportion of total roots in each soil layer (Kuppel et al., 2018b). The vertical distance (v) from the base of the vegetation to the vertical centre of the biomass in each layer was estimated using a log function of rooting proportions modified from Sperry et al. (2016, 2022) to account for unequal rooting proportions with layers using cumulative rooting proportions:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}v\left(i\right)& =\mathrm{0.01}\cdot \ln (\mathrm{1}-\mathrm{0.995}\cdot (p\left(i\right)\\ & -\left(p\right(i)-p(i-\mathrm{1}\left)\right)/\mathrm{2}\left)\right)/\ln \left(\mathit{\beta }\right),\end{array}\end{array}$$

where i is the soil layer, p is the cumulative proportion of roots (from the surface), 0.995 is a coefficient indicating 99.5 % of root biomass, and β is a rooting distribution parameter. The rooting parameter in Sperry et al. (2016, 2022) was estimated using a linear correlation to the EcH₂O rooting distribution (k_root) to enable estimation in EcH₂O ($\mathit{\beta }=-\mathrm{0.0089}\cdot k_{\mathrm{root}}+\mathrm{0.9947}$). The radial distance from the stem of the vegetation is estimated as a function of the volume of rooting (Sperry et al., 2016, 2022):

$$\begin{array}{}\text{(2)}& \mathrm{Vol}=d_{\mathrm{1}}\cdot \mathrm{\Pi }\cdot (D\cdot a{)}^{\mathrm{2}},\end{array}$$

where Vol is the volume of roots in layer 1; d₁ is the depth of layer 1; D is the total depth of the soil (equivalent to total root depth); and a is an aspect parameter, controlling radial distance. The radial distance is then estimated using the proportion of roots in the soil layers (modified from Sperry et al., 2016, 2022):

$$\begin{array}{}\text{(3)}& r\left(i\right)=(\mathrm{3}\cdot \mathrm{Vol}\cdot p(i)/(d\left(i\right)\cdot \mathrm{\Pi }){)}^{\mathrm{0.5}},\end{array}$$

where r(i) is the average radial distance in layer i, p is the proportion of roots in layer i estimated from k_root, and d is the depth of layer i (Fig. 2a). The total average distance (D(i)) for roots in each layer is the sum of v(i) and r(i). The radial extent of roots in each layer is then used to estimate the proportion of roots within the cell and in adjacent cells.

3.3.2 Xylem transit time and calibration

Since much of the rooting system below trees is not well characterized, including limited information on the proportion or total lengths of fine roots to transport roots being available, and due to the translation of xylem velocity to fine root velocity, we simplified the conceptualization of root transport for each soil layer into convolution integrals (Rothfuss and Javaux, 2017). We conceptualize the rooting system in each soil layer as being similar to a large river (primary root transport) network with many tributaries (fine roots) with precipitation input (root uptake volume). In this way, root water velocity is implicitly accounted for (averaged) across different root diameters. Root uptake volume-weighted (combined in cell and off cell as estimated from EcH₂O-iso, Fig. 2c) convolution equations (gamma distribution) were used to estimate xylem isotope values from each soil layer for each hourly time step (Fig. 2d). Fractionation during uptake and along the flow path to measurement height was assumed to be negligible. The approach introduces a scale parameter (β, hours) for each soil layer and one shape parameter (α) as a function of total root length (four total parameters) (Fig. 2b). This approach assumes that there are no additional effects of cavitation apart from sap velocity decreases (measured or simulated) within the xylem or roots of the plant.

Figure 2 The estimation of rooting density (a) vertically and (b) the transit time distribution describing root velocity and distance to measurement height weighted by (c) sap flow to get the (d) cumulative contribution of soil water from entry time (red line).

To evaluate and test the information content of sampling frequency, the convolution integral was calibrated against xylem isotopes at different time step intervals: using 6, 12, and 24 h intervals. Measured and simulated xylem isotopes were averaged over the respective intervals prior to evaluation with the Kling–Gupta efficiency (KGE, Kling et al., 2012). While not used in calibration, the KGE of standardized hourly xylem was evaluated to examine how well sub-daily dynamics of soil water isotopes propagated through the roots. Lastly, to evaluate the effectiveness of distance-based and instantaneous mixing, the transit time and xylem isotopes were calibrated twice, (1) using modelled soil isotopic compositions and sap flow and (2) using measured soil water isotopes and sap flow. The use of measured soil water isotopes and sap flow tests the maximum potential for how each model performs and is not limited to the performance of EcH₂O-iso for sap flow or soil water isotopes. Further, the twin calibration approach shows how the error in the simulation of soil water isotopes potentially propagates through to the estimation of xylem isotope and transit time mixing. Both xylem isotope model calibrations used the root water uptake proportions estimated from the multicriteria EcH₂O-iso calibration and are inclusive of off-cell root water contributions (Sect. 3.4). The Akaike information criteria (AIC, Akaike, 1998) was used to assess the significance of the additional parameters used by the distance-based mixing (four additional parameters). The AIC is estimated using the log likelihood of each model fit and number of parameters, where a lower AIC indicates better model performance.

3.4 EcH₂O-iso and tree mixing model calibration and set-up

3.4.1 Model set-up

To best leverage measurements in both willows and the soil pits, the model was set up with three square grids (6 m). Two grids contained one willow each, and the third contained the open grass area (Fig. 1c). The model was set up to run with hourly time steps between 1 January and 31 October 2020, encompassing a spin-up and the primary growing season. Soil layer depths were fixed at the depth of the soil moisture measurements (10, 40, and 100 cm) for all model cells. Further, previous analysis and observation during the soil moisture probe installation suggest that primary water sources of the willow trees were within the upper 100 cm of the soil (Landgraf et al., 2021), which is consistent with shallow rooting of willow trees observed elsewhere (e.g. Phillips et al., 2014). Additionally, there is little isotopic evidence of trees using water with distinct groundwater or stream water signatures, despite their proximity (Landgraf et al., 2021). Initial testing revealed inadequate water supply for willows within a single cell (from all soil layers), while measurements revealed no notable decrease in sap flow during the growing season. This suggested potentially substantial water use from outside the measurement location (i.e. model cell), and therefore parameterization of rooting radius (aspect parameter) was set to extend beyond cell boundaries (radius > 3m) and permits water use from adjacent cells (assumed to be grass in adjacent cells outside the model domain, Fig. 1c). Increased rooting radius increases the proportion of water used from outside the model cell. Additionally, root distributions with depth were parameterized to be in synchrony with common willow rooting distributions, with a greater proportion of roots near the surface (k_root > 5). Due to the notable differences in shallow soil moisture between Site A and Site B (soil organic matter content), soil parameters were different below the willows and the open grass (see calibrated parameters for soil and vegetation in Table S1 in the Supplement). Initial soil moisture was set to field capacity, which was shown by model testing to have negligible influence on soil moisture simulations during the growing season. Soil water isotopes were initialized using the average measured isotopic composition of soil water. Model calibration of the shallowest soil depths was not sensitive to initial soil water isotope conditions, as simulated isotopes stabilized prior to the beginning of measurement (end of May 2020). The basal diameter of the willows was set to the measured diameter at the start of monitoring.

3.4.2 Model forcing data, calibration, and evaluation

Hourly model forcing data (Table 2) primarily consisted of data measured with the mobile weather station at the study site (Sect. 3.1) corrected and gap-filled (where necessary) using the surrounding weather stations (Fig. S2). As long-wave radiation was not measured, estimated long-wave radiation was used from the ERA5 reanalysis dataset (Hersbach et al., 2020) (Table 2). With no direct measurements of LAI in the grass site, MODIS Terra (MOD15A2, Myneni et al., 2015) dynamics were used to calibrate the LAI (calibration step 2). Due to the coarser resolution of MODIS (500 m), the LAI range for the grass was scaled to have a maximum of 2 (cf. Smith et al., 2021). In the willow calibration, the MODIS dataset was used in conjunction with canopy images (Sect. 3.1) in willow LAI calibration.

To best constrain the model performance and evaluate how well the model estimates the dynamics of water, isotope, and energy fluxes during the growing season, all available data (except δ¹⁸O and soil temperature) were used to constrain model parameterization, with δ¹⁸O retained for further validation (Fig. S3). To simplify the presentation of results, we present δ²H, as the temporal dynamics of deviation from measurements were not greatly different from δ¹⁸O (Supplement). As the full duration of the single growing season was necessary to constrain fluxes, energy, and biomass (e.g. biomass growth and decay) and the model was not used for predictive purposes, a temporal split of the datasets for calibration and validation was not practical. Model results were informally validated using retained datasets (δ¹⁸O, sensible heat, and soil temperature) and statistical analyses (root uptake distribution). To simultaneously evaluate the variability, bias, and correlation of the measured and simulated datasets, the KGE was used for multicriteria calibration of all datasets. Calibration was conducted in two steps, (1) calibration of the isotope, energy, and water balance and (2) calibration of biomass indices (leaf area index and basal area). To minimize the effects of biomass calibration on the energy and water balance, the temporal dynamics of LAI from MODIS Terra were used as input time series. The second step of calibration resampled parameter sets from calibration step 1 for soil and non-biomass vegetation parameters and coupled them with Latin hypercube sampling of biomass vegetation parameters. Each calibration utilized 100 000 parameter sets generated using Latin hypercube sampling, with parameter combination feasibility evaluated prior to simulations. Infeasible parameter combinations identified using standard ranges of soil properties (e.g. field capacity) were rejected and resampled from parameter space. An empirical cumulative distribution function (eCDF) was created for each output. For each parameter set, the minimum eCDF value for all variables was used to rank the simulations, retaining the 100 “best” simulations for analysis (Ala-aho et al., 2017). Simulation ranking checked for soil parameter uniqueness to ensure no repeated sampling of soil and non-biomass vegetation parameters occurred within the 100 best simulations.

Table 2 Forcing and calibration data of EcH₂O-iso calibration. Round brackets in the calibration datasets indicate the spatial locations of all available measured data.

4 Results

4.1 EcH₂O-iso soil water, isotope, and energy balance evaluation

Model calibration provided reasonable approximations of the measured variables (Table S2), with the average goodness of fit of model simulations within measurement uncertainties (soil moisture MAE – mean absolute error – range of 0.01–0.06 m³ m⁻³, Table S2). Similar to measured data, simulations of soil moisture and soil water isotopes showed decreasing dynamics with increasing depth (Fig. 3). Despite the close proximity, Site A and Site B had largely different soil moisture and isotopic response in the upper soil layers (Fig. 3a–d). Simulations were able to capture the drier soil moisture in Site A relative to Site B (Fig. 3a and c); however, the simulated dynamics of average moisture in layer 1 were greater than the measured moisture at 10 cm. The sub-discretization of the shallow layer 1 soil (1 cm increments using EcH₂O soil water redistribution) using calibrated parameters revealed that while average moisture in the upper 10 cm is high (blue lines, Fig. 3a), percolation of infiltrated water down to 10 cm during the summer months is limited (Fig. 3a, red line). Late-season soil moisture values (average and at 10 cm) were overestimated relative to measurements, coinciding with the underestimation of sap flux during the same period (Fig. 4).

The general dynamics of the soil water isotopes at 10 cm in Site A were reasonably represented, with the primary deviation of simulations from measurements due to large sub-daily measured isotope variability (δ²H MAE = 8.6 ‰, Table S2). Additionally, the period of soil water isotopic depletion in August was not captured. Soil moisture and isotope dynamics in the deeper soils (40 cm) were much more damped, with a slight underestimation bias at 40 cm at Site B despite appropriate dynamics (Fig. 3d). The much lower variability in isotopes and soil moisture at 100 cm was similarly captured by the model calibration. The higher depletion simulated at Site A relative to the measured values was consistent with the soil water isotopic values at greater depth (100 cm) measured at Site B. Outliers of in situ soil water isotopes towards the end of the growing season (end of August and September, Fig. 3c–f) were a result of rapid temperature changes which were too strong for the system (i.e. heated cables) to control and caused temporary condensation effects in the tubing. These data were few and did not affect the soil water isotope modelling, and uncalibrated soil water isotopic simulations (δ¹⁸O, Fig. S3) were broadly similar.

Figure 3 Simulation mean and upper and lower bounds of soil moisture and δ²H at Site A at (a) 10 cm, (c) 40 cm, and (e) 100 cm and at Site B at (b) 10 cm, (d) 40 cm, and (f) 100 cm. The top row shows the hourly precipitation for both sites. Isotopic values (circles) are average daily measurements, with bars showing the daily range in isotopic measurements. δ¹⁸O simulations are shown in Fig. S3.

Estimated energy balance (latent heat, sensible heat, transpiration flux, and temperature) was shown to match in situ measurements within the limits of measurement uncertainty (Fig. 4, Table S2). Despite only utilizing latent heat and surface temperature at Site B in the calibration, simulated sensible heat and soil temperature at multiple depths were reasonably captured (MAE of 4.5 W m⁻² and 1.2–2.8 ^∘C for sensible heat and soil temperature range, respectively, Table S2). Latent heat and sensible heat at Site A (Willow 2) were noticeably more variable than at Site B, driven by the greater transpiration demands of the willow relative to the grasses. Sap flux in the willow was adequately captured by the model (MAE = 0.02 m³ d⁻¹), with the primary deviation – underestimation – from measurements during the late growing season (September–October). Greater variability in surface and soil (10 cm, Fig. S4) temperature was additionally simulated at Site A; however, estimated differences in ground heat storage between Site A and Site B resulted in the increased damping of soil temperature with depth at Site A. Progressive overestimation of soil temperature during the early simulation period with increasing soil depth is likely due to the propagation of a slight overestimation of early-growing-season surface temperature at each site (Fig. 4) potentially due to a lower estimated leaf area index at the beginning of the growing season (Fig. 5).

Figure 4 Simulation mean and upper and lower bounds of energy balance components of (a, b) latent and (c, d) sensible heat, (e, f) volumetric vegetation water use, and (g, h) surface and (i, j) soil temperature (40, 10, and 100 cm in Fig. S4). Black lines show the mean of measured data with grey-shaded areas showing uncertainty and sub-hourly variability. Volumetric vegetation water at Site A shows sap flow within the willow, and that at Site B shows total volumetric water use of the grass. Note that soil temperature was not calibrated, and transpiration rates are shown in Fig. S4; sap flow for Willow 1 is identical to Willow 2 (Fig. 4e).

4.2 Growing season dynamics – fluxes and biomass accumulation

Simulated biomass production for foliage and tree diameter produced reasonable results for the two willow trees throughout the simulation period, with only minor deviations (willow LAI MAE = 0.5 m² m⁻², Table S2) from both measured and remotely sensed (MODIS) datasets (Fig. 5a and b). Average simulations (solid lines) underestimated LAI at the beginning of the growing season; however, these differences were relatively small, and uncertainty bounds of simulated LAI were within measurement uncertainty. Calibrated grass leaf area index dynamics were consistent with MODIS, with only a small deviation in leaf area index at the beginning of the growing season (MAE = 0.1 m² m⁻²; Table S2; Fig. 5c). The dynamics and total net stem growth was adequately captured by the model for each willow (MAE = 1.5–2.2 µm; Table S2); however, there was a slight underestimation of average net growth in Willow 2 (Fig. 5c and d).

Figure 5 Simulated and measured leaf area index of (a) Willow 1, (b) Willow 2, and (c) grass and net growth of willow diameter of (d) Willow 1 and (e) Willow 2.

Modelled partitioning of water fluxes and biomass allocation in the willow and grass sites showed marked differences despite the relative proximity (Table 3). Evapotranspiration (ET) dominated water partitioning throughout the summer, with the greatest ET between 15 June and 15 September. Total ET was dominated by transpiration (T_r) particularly for the willows (> 70 %, Table 3). This coincided with a decreased proportion of interception evaporation (E_i), and soil evaporation (E_s) from drier soils. Pre-growing-season infiltration was higher below the willows relative to the grass, with moderately lower infiltration below the willows during the growing season. The lower density of grass coverage (and increased bare soil) below the willow reduced growing season differences in infiltration between the willow and grass sites. The high root uptake rate of the willow resulted in a negligible percolation of water to deeper soil layers during the growing season. Continuation of drier conditions below the willow resulted in a gradual shift toward near-surface root uptake compared to a higher proportion (24 %) from deeper soils at the beginning of the growing season (Table 3). The willows showed little change in biomass allocation, with only a slight shift in root growth in the pre-growing-season period to foliage growth during the peak growing season (Table 3). Allocation to stem growth showed negligible change throughout the growing season. The grasses had much higher percolation to deeper soils relative to the willows as a result of low ET. Wetter soils below the grasses resulted in greater soil evaporation proportions. Higher soil moisture throughout the soil profile resulted in more stable root uptake proportions from layers 1 and 2, with an approximately equal mixture of near-surface and mid-depth (40 cm) soil water. Biomass allocation in grasses showed a greater allocation of carbon to foliage growth than roots (note that there is no allocation to stem growth).

Table 3 Mean (and standard deviation) of water partitioning before (January–May) and during the early growing season (1 May–15 June), mid growing season (16 June–1 August), and late growing season (1 August–31 October). Also shown are the root uptake (RU) proportions from each soil layer (L1, L2, and L3 are layers 1, 2, and 3, respectively) and biomass allocations. Depth of the L1, L2, and L3 are 0–10, 10–40, and 40–100 cm, respectively. Note that infiltration and percolation are cell averages; ET is specific to vegetation (proportion of willow in cell < 100 %). ET is evapotranspiration; E_s is soil evaporation; T_r is transpiration; and E_i is interception evaporation.

4.3 Evaluation of xylem water mixing

To avoid influencing the evaluation of tree water mixing by the results of simulated soil isotopes (lower variability than measured, Fig. 3a), tree xylem-water mixing was conducted in separate calibrations using both simulated and measured soil water isotopes. The instantaneous mixing of root uptake water throughout the tree (instant uniformity in all xylem at each time step, as simulated by the EcH₂O-iso model structure) was able to capture the general seasonal dynamics of xylem water in both willows (Fig. 6a and c), though simulated diurnal variability could not fully reproduce the measured variability. Unsurprisingly, the instantaneous mixing approach showed improved performance when measured soil water isotopes were used rather than simulated soil water isotopes, increasing seasonal and day-to-day variability (Table 4; Fig. 6b and d). The lower AIC for the distance-based mixing simulations of xylem water at 1 m suggested that the additional parameterization aided in the simulation performance of xylem water isotopes relative to the instant mixing simulations. The decrease in KGE in Willow 1 from simulated to measured data (Table 4) was an artefact of the shorter measured soil water isotope time series (stopping before the end of October), which did not encompass the large depletion at the end of October (time series length incorporated in AIC). Differences in the AIC and KGE results using measured datasets relative to simulated datasets (Table 4) were due to differences in the method of goodness of fit, whereby squared differences (in AIC) were minimized by the distance-based approach but were accompanied by a penalization of the variance coefficient (in KGE) due to more damped dynamics than the instant mixing approach.

Importantly, xylem simulations with instant mixing and modelled soil isotopic composition were unable to capture the observed standardized hourly variability in the xylem isotopic composition (right panels of Fig. 6a and c). This is consistent with the limited standardized hourly variability in simulated soil water isotopes relative to measured soil water isotopes (Figs. 3a and 6a and c). Distance-based simulations with modelled soil water isotopes could reasonably capture the timing of sub-daily dynamics due to the sub-daily variability in sap flow volumes. Xylem water mixing estimations using measured soil water isotopes revealed a reasonable standardized hourly variability in xylem isotopes (Fig. 6b and d). Instant mixing unsurprisingly showed dynamics similar to soil water isotopes, which peaked before the standardized measured xylem water. Distance-based mixing with measured soil water isotopes was able to capture the sub-daily variability and timing of the peak; however, the rapid decrease (and late peak) in sub-daily standardized soil water isotopes at 10 cm slightly degrades the simulated dynamics late in the day.

Model performance improved when the mixing model was calibrated to daily compared to sub-daily (6-hourly or 12-hourly) data (higher KGE). While simulations show deviations from the temporally averaged xylem water isotopes (e.g. 12 h average, Fig. 6), simulations were generally within the minimum and maximum measured range over the time step. The increase in model performance is consistent with the reduced variability in measured xylem with increasing time steps, increasing in similarity to the variability observed in soil water isotopes. Likewise, the KGE increased most rapidly with time steps on average when simulated soil water isotopes were used in mixing, as they lower variability more than measured soil water isotopes (Fig. 3).

Figure 6 Simulated 12 h average xylem deuterium using (a) simulated and (b) measured soil deuterium with distance-based and instant mixing in Willow 1. Simulated 12 h average xylem deuterium using (c) simulated and (d) measured soil deuterium with distance-based and instant mixing in Willow 2. Measurement bars represent the maximum and minimum measured xylem water isotopes during the 12 h period. For each simulation, standardized sub-daily δ²H (S(δ²H)) shows the measured and simulated soil (cell-contribution-weighted average) and xylem average hourly variability. Note that standardized simulated soil water isotopes overlap for all soil layers.

Table 4 Akaike information criteria (AIC) and average KGE (± standard deviation) for model performance for each willow evaluated at different time steps. Values in round brackets for KGE show the KGE of sub-daily variability. Mixing models used either output (isotope and sap flow) from EcH₂O-iso or measured isotope and sap flow. All models used estimated uptake proportions from EcH₂O-iso calibrations.

4.4 Evaluation of fluxes, uptake, and water ages through the growing season

Differences in infiltration and percolation quantities between the willow and grass (Table 3) further propagated to water ages below each vegetation type (Fig. 7a and b). Average layer 1 water ages below the willow (15 ± 1 d) were similar to the grass (17 ± 1 d), with larger proportions of 7 d (36 ± 2 %), 14 d (59 ± 2 %), and 30 d (84 ± 1 %) relative to those below grass (31 ± 1 %, 52 ± 1 %, and 67 ± 1 %, respectively). The higher proportions of young water (more recent precipitation) below the willow were primarily the result of lower moisture storage, causing a greater effect of infiltration on water ages. The limited percolation to layer 2 below the willow resulted in a continuously increasing water age in the deeper soil water. Percolation throughout the summer under the grasses resulted in a stabilization of water age (average of 73 ± 2 d), with 81 ± 1 % of water older than 30 d. Average root water age in the root tip was a mixture of water ages for all soil layers, with the high root uptake from shallow soils by the willow resulting in a younger average uptake water age in the willow roots (28 ± 6 d) compared to the more equal mixture of shallow and mid-depth soils taken up by the grass roots resulting in uptake of older water (51 ± 3 d) (Table 3).

Water ages in storage and uptake were further discretized into contributions of water from each monthly precipitation amount to characterize the seasonal origins of water stored and utilized. Fast turnover of water in the shallow soils (i.e. upper 10 cm) resulted in low percentages of late winter and spring water (before 1 April) in the early growing season (13 % and 18 % below the willow and grass, respectively, Fig. 7 and Tables S3 and S4), which were negligible at the end of the growing season (0 % below the willow and grass, Tables S3 and S4). Limited percolation from shallow to deeper soils under the willow resulted in large proportions of deeper water originating from spring or winter precipitation (77 %), while under grass the spring and winter precipitation only accounted for 37 % of deeper soil water by the end of the growing season. Differences in rooting distribution and the temporal dynamics of infiltration resulted in differences in spring and winter water usage for the grass and willow throughout the growing season, with 33 %, 16 %, and 9 % for the willow in the early, mid, and late growing season (57 %, 35 %, 24 %, respectively, for grass) root uptake of spring/winter water.

Figure 7 Average water age (and range) and fractional monthly water contribution in layer 1 (10 cm) in (a) Willow 1 and (b) grass and in layer 2 (40 cm) in (c) Willow 1 and (d) grass; the average (and range) of water entering roots in (e) Willow 1 and (f) grass; and (g) the average (and range) of time between root uptake and 1 m using distance-based mixing using simulated or measured soil water isotopes. Fractional water ages are shown in Fig. S5, and transit times from different soil layers are in Table S5.

Since instantaneous mixing equates water age through all vegetation storage and pathways (i.e. roots and xylem), xylem ages are only shown for the distance-based mixing approach (using simulated and measured soil water isotopes). Using simulated soil water isotopes, distance-based mixing revealed an average transit time (all time steps) of 149 ± 30 h from the roots to 1 m height (Fig. 7g). Distance-based xylem mixing using measured soil water isotopes revealed a mean transit time of 237 ± 97 h (Fig. 7g). Smaller uncertainty in transit time estimated from simulated soil and sap flow was a direct result of smaller temporal variability (lower degrees of freedom) of simulated soil water isotopes compared to measured soil water isotopes.

5 Discussion

5.1 Assessing growing season flux–storage–tracer relationships, water partitioning, and biomass dynamics

The utility of ecohydrological modelling as a tool to evaluate the linkages between water partitioning and biological dynamics is directly related to the capability of such models to accurately reproduce a wide range of variables to indicate the acceptable representation of ecophysiological processes. At fine spatial scales, such as plot sites, reproducing flux–storage–biomass dynamics and feedback mechanisms is important for the transference of their nonlinearities to larger scales (Asbjornsen et al., 2011). Advancements of complex (e.g. physically based) models and fusion with multiple, high-temporal-resolution data streams can lead to the convergence of model performance for all variables (Asadzadeh et al., 2014) and higher confidence in the ability of the model to perform across a wide range of hydroclimatic conditions. The application of a tracer-aided ecohydrological model on sub-daily time steps using multiple data streams from an intensive in situ monitoring site provided a means to directly evaluate how well the model could simultaneously reproduce flux–storage–tracer–biomass dynamics through multiple criteria throughout the high variability in a growing season.

Overall, the model performed well against all measured data in both the willow and neighbouring grass sites, including soil water storage, water and energy fluxes, isotope variations, and biomass dynamics throughout the growing season. Total green water usage (total ET) was notably higher for the willows compared to the grass, similar to other nearby studies (Douinot et al., 2019; Kleine et al., 2021; Smith et al., 2020b). Further, the proportions of ET for interception evaporation (as a proxy for total interception) and transpiration were well within expected ranges for both vegetation units (0.4–0.8 = $T_{\mathrm{r}}/\mathrm{ET}$) (Coenders-Gerrits et al., 2014; Dubbert et al., 2014; Schlesinger and Jasechko, 2014; Zhou et al., 2016), including lower $T_{\mathrm{r}}/\mathrm{ET}$ for the grass compared to the willows. This suggests adequate reproduction of water partitioning within the model. Despite the relatively dry characteristics of the soil at and below the soil moisture measurement at 10 cm, transpiration in both willows and grass was maintained throughout the study period (Fig. 4), albeit with an underestimation of the transpiration in the willows toward the end of the growing season. The underestimation was accompanied by an overestimation of soil moisture but an appropriate depletion in the shallow soils below the willow (Fig. 3a). While the late-season deviation of estimation of sap flow and soil moisture is directly linked, the model structure does not permit temporal variability in lateral root water usage (i.e. outside the model cell); therefore increased late-season transpiration could have detrimental effects on shallow soil moisture below the grass. Furthermore, an increase in late-season transpiration would only be the result of vegetation parametric changes outside feasible ranges. The high goodness of fit for the primary growing season for the willow trees was mainly possible within EcH₂O-iso due to the capability of the willows to use water from adjacent model cells. As observed in other nearby studies, the soil evaporation contribution to ET under both the willow and grass was damped throughout the growing season due to leaf shading (LAI) reducing energy availability at the surface (Kleine et al., 2021; Smith et al., 2021). Shading effects of the willows with higher interception capacities could be the cause of the late-season stability in isotopic measurements below the grass due to reduced precipitation input compared to the simulated large depletion of measured and simulated soil water isotopes below the willows. Although willows are generally a riparian species able to access groundwater or surface water, the willows at this site did not reveal contributions from the groundwater (deeper than in situ soil water isotope measurements, > 2m). The predominantly shallow root uptake contribution (driven by higher near-surface water availability and root distribution) was consistent with root distributions observed in willows elsewhere (Cunniff et al., 2015; Phillips et al., 2014) and only slightly lower than independent Bayesian mixing estimations conducted using soil and xylem isotopes for the same trees (Landgraf et al., 2021). The near-equal contribution of shallow and mid-depth uptake in the grass (and subsequent older water age in mid-depth soils) was more prominent than expected given the shallow rooting zone. The deeper rooting distribution in the grass is likely due to a combination of vegetation functional types (i.e. willow vs. grass) and high uptake competition in the near-surface soils with the nearby high-water-use willows (Moustakas et al., 2013). Competition could potentially drive deeper rooting of the grass, contributing to the more dynamic soil moisture at 40 cm as measured below the grass. The slight decrease in simulated root biomass production in the willows during the summer months suggested that vegetation was not under water stress in the study period, as allocation increased for the leaf (water usage) rather than the root (water sourcing) (Eziz et al., 2017). Average hourly simulated transpiration dynamics were comparable to sap flow measurements (Fig. S6), including small but notable uptake during the night, similar to observations in other woody plants (Dawson et al., 2007).

Within this study, the relatively low soil moisture did not propagate to a notable decrease in sap flow. While the low shallow soil moisture and the high proportion of near-surface roots of the willows seemed inconsistent with the measured high vegetation water usage, in situ isotope dynamics were essential in confirming that measurements of shallow moisture were not underrepresenting infiltration events (i.e. via spatial heterogeneity). Multicriteria calibration using high temporal resolution and highly dynamic soil water isotopes revealed a strong positive correlation of the KGE of simulated shallow soil moisture and the KGE of simulated shallow soil water isotope (i.e. model performance of soil water isotopes only improved with the performance of soil moisture). Bulk soil water isotope samples, with a coarser temporal resolution, would have been unlikely to adequately capture the necessary isotope dynamics for calibration. EcH₂O-iso was only able to adequately estimate soil moisture at 10 cm through sub-discretization of the soils above 10 cm, which suggested high retention of water in soils above the measurement depth (e.g. 0–5 cm), taken up by roots or soil evaporation before there is sufficient water to percolate to 10 cm. These small-scale spatial variations, as well as the large spatial differences between the soils below the willows and below the grass and between different soil layers (Fig. 3), reveal the significant heterogeneity of the site despite relatively immature soils and the local spatial scales. The heterogeneity of moisture and isotopes is likely further exacerbated by spatial heterogeneity under the canopy, as differing branch and leaf structures can impact detailed throughfall distribution (Dalsgaard, 2007; Gerrits et al., 2010; Li and Liang, 2019). If such data were available, it could likely further improve soil moisture simulations, as zones of higher soil moisture may be present at locations below the canopy where throughfall is concentrated (Gerrits et al., 2010), which could in turn influence and increase dynamics of near-surface soil water isotopic compositions. These heterogeneities would likely help improve the model performance of near-surface soils by accounting for additional variability in measurements. However, given the relatively small quantity of water encompassed in the shallow soils relative to total fluxes, further heterogeneity of measurements in soils would not likely influence the overall estimation of evapotranspiration fluxes (or partitions) or root water sources.

5.2 Evaluation of mixing dynamics of root uptake and implications of the rooting zone

We used output from EcH₂O-iso (i.e. the root uptake sources) and a simple, parsimonious approach to further explain the time delay and variability in xylem water isotopic composition from soil water isotopic composition. The approach showed substantial improvements over the use of the instantaneous mixing currently applied within most modelling approaches. While the instantaneous mixing had occurrences of higher KGE than the distance-based approach, this was the product of greater simulated variability rather than improved dynamics, as indicated by the AIC. Further, the additional capabilities of the distance-based model to capture the lagged response of diurnal soil water isotopes in xylem further underlines the advantages of the approach. The inclusion and calibration of rooting length into EcH₂O-iso showed comparable results to other willows (1.3–1.4 vs. 1.3 tree height/root length; Table S5; Phillips et al., 2014), which supports the initial estimation of spatial uptake sources. With around half of the uptake (by root length and water availability) estimated to occur outside of the willow cells, this consideration of spatial variability in soil water isotope composition underlines the potential influence of spatiotemporal variability in source waters on xylem isotopes, as suggested by De Deurwaerder et al. (2020) and Seeger and Weiler (2021).

As with partitioning soil moisture and vegetation water usage, the high-resolution in situ isotopic composition of xylem water was indispensable for constraining the mixing dynamics of root uptake. Coarser temporal destructive xylem sampling could potentially overestimate xylem variability due to the large sub-daily variability, damping seasonal patterns (e.g. enrichment in September, Fig. 6) and influencing estimated root uptake profiles. Further, while the use of tracer-aided ecohydrological model output for estimating xylem water composition with refined rooting distributions was able to reconcile general seasonal xylem dynamics and some sub-daily variability due to sap flow (7 d moving average and sub-daily variability, Fig. 6), deviations in seasonal magnitudes of simulated to measured xylem water isotope dynamics were primarily due to differences in simulated vs. measured soil water isotopes (e.g. Fig. 6a vs. Fig. 6b) in the shallow soils (Fig. 3 and RU in Table 3). This additionally highlights the importance of in situ isotopic soil and vegetation measurements to constrain model estimates. The use of measured data within the same framework further showed the capabilities of the approach to estimate xylem dynamics, including more direct translation of soil sub-daily variability in the xylem (Fig. 6) as suggested by von Freyberg et al. (2020). However, even the use of in situ measured soil water isotopes was unable to fully resolve the xylem dynamics with much larger day-to-day variability than could be estimated by the modelling approach (Fig. 6, Table 4).

Uncertainty in isotopic measurements, complex physiological processes within the trees which may cause fractionation, limitations of the model structure, and interactions between all of these were potentially limiting factors for reducing the xylem mixing model performance. In situ measurement uncertainty, particularly in xylem, may be large for δ²H (5 ‰–10 ‰; Beyer et al., 2020), greater than the standard deviation of the daily mean δ²H in xylem (3.7 ‰). Variation in the daily mean is additionally lower than the sub-daily standard deviation (6.8 ‰), which would suggest sub-daily variations are more significant than daily variations. These daily variations in xylem δ²H could be influenced by natural short-term climatic variability and monitoring errors, with condensation, mixing, and diffusion potentially causing instabilities in measurements and greater uncertainties (Beyer et al., 2020; von Freyberg et al., 2020). In this study, condensation uncertainties were minimized by heating the tubing and were largely restricted to the autumn (September and October) when nighttime temperatures could be much cooler (Landgraf et al., 2021). Further, while fractionation during uptake or within vegetation has frequently been discussed and re-examined recently (Brinkmann et al., 2019; Poca et al., 2019; Vargas et al., 2017), the vegetation isotopes reflect the fractionated soil water isotope sources (Landgraf et al., 2021), which suggests that vegetation-specific fractionation is not likely a major cause of isotopic variation on either hourly or daily time steps at our site. Finally, the model structure may play a role in the xylem mixing estimation, in particular with the assumed “static” transit time (consistent with the SPAC supply–demand models; e.g. Sperry et al., 2016; Simeone et al., 2019), limiting effects of dynamic tree cell water storage exchange with xylem. These water pools may provide a different source of contrasting isotopic composition to sap flow when transpiration is lower (De Deurwaerder et al., 2020; Secchi et al., 2017) and may make a significant contribution to transpiration fluxes at certain times (Urban et al., 2014). The linkage of the model structure of flow path mixing as shown here, as well as tree cell water storage, inclusive of simple parametric storage mixing approaches (e.g. Knighton et al., 2020), presents potential further model development to better estimate sub-daily and inter-daily xylem variability and interpret water use dynamics in vegetation. While such a release of water may be present during peak transpiration hours, as suggested by the hourly variability in measured basal diameter (Fig. S7), the high water usage of willows potentially limits the proportional contribution of cell storage to xylem water.

5.3 Implications of water ages in soil water and xylem fluxes

Quantifying the age of water has significant implications for improved understanding of how water is cycled through the critical zone, particularly for understanding the temporal changes and response times of water quantity and quality in different environments (Sprenger et al., 2019). Further discretization of average water ages into age distributions helps to illuminate how the evolution of water ages are influenced by different processes (Rodriguez et al., 2020), as well as the influence of water fluxes from different temporal periods or extremes (Allen et al., 2019). The modelled contributions of different water ages in the willow root uptake suggested that summer precipitation was the dominant source of water, accounting for 56 % of total uptake (volume-weighted). This differs from isotope-based studies in other locations (e.g. Allen et al., 2019, across Switzerland) which have shown a higher usage of winter precipitation in leafy vegetation (beech and oak), but it is consistent with the estimated high contribution of recent precipitation during the growing season in the study region (Miguez-Macho and Fan, 2021). However, differences in hydroclimate, as well as variations in rooting depths between species and the maturity of forests stands, likely give deeper rooting depths for beech and oak compared to shallow roots for the younger willows (in this study), which could drive differences in seasonal precipitation use by vegetation. Alternatively, the effect of separation of storage with fast and slow turnover times or the effects of preferential flow (Sprenger and Allen, 2020) could contribute to an increase in water age if younger summer water preferentially moves through the soil or if tightly bound water is used by vegetation. However, the good temporal agreements of soil and xylem isotope samples relative to the likely dichotomy between faster and slower moving isotopes in the soil (e.g. Sprenger et al., 2019) imply that for the conditions present at the study site, inclusion of preferential flow or variable pore size would not influence water ages for the uptake from either of the soils. This is also consistent with the immature nature of the soils and their relatively uniform sandy nature.

Age estimates for water being transported in xylem derived in this study were of a similar order of magnitude to estimates for other woody vegetation, averaging between 6.2 and 8.1 d with a maximum of 10.3 d at our site compared to between 2 d and 1 month in other studies (Brandes et al., 2007; Meinzer et al., 2006). While our estimates were towards lower transit times compared to most other studies, the willows at our site were younger, had smaller diameters, were measured at 1 m basal height compared to in the canopy, and had higher xylem sap flow velocity, consistent with high water usage and low stomatal control (xylem cavitation control; Wikberg and Ögren, 2007). Additionally, differences in xylem water ages relative to previous root uptake mixing studies (e.g. Seeger and Weiler, 2021) are likely due to stand size, vegetation age, and measurement heights. High-resolution in situ measurements, capturing soil and xylem isotope dynamics, were indispensable for improving the confidence in the estimated xylem transport age. Coarser xylem isotope samples would likely be insufficient to adequately constrain transport age uncertainty. Similar to the impact that stored tree water release can have on xylem isotopes (Sect. 5.2) (Brandes et al., 2007; De Deurwaerder et al., 2020), there are likely implications for water ages and the identification of sources due to water released from cell storage. The effect of stored water on xylem transport ages could be significant in vegetation that utilizes large fractions of previously stored water (e.g. which can be as much as 20 % of transpiration; Čermák et al., 2007). In such vegetation, total transpiration could be maintained from stored water for a week (Čermák et al., 2007), drawing water from storage cells which may be notably older than xylem water. Effects of such long-term storage components have previously been implicitly incorporated into post-root-uptake lumped vegetation mixing approaches (e.g. Knighton et al., 2020), which have previously shown order-of-magnitude changes in tree stored water age with changes in wetness conditions. Due to the relatively young uptake water ages of willows in this study, the inclusion of any cell water storage (∼ 16 L d⁻¹ if 20 % of transpiration, ∼600 L total storage to maintain transpiration for 1 week) may significantly increase the mean age of xylem water and would increase the total vegetation mixing time. Given the limited studies describing the ages of water and tracers in transpiring trees (Sprenger et al., 2019; Mennekes et al., 2021; Benettin et al., 2021), more studies quantifying the transport through the xylem, as conducted here, and the inclusion of additional mixing within vegetation (e.g. Steppe et al., 2006) would be beneficial to further constrain plant water use estimations.

6 Conclusions

In order to increase understanding of water cycling mechanisms in the critical zone, it is essential to evaluate how water is partitioned by vegetation and the dominant processes controlling water usage and movement. Furthermore, quantifying flux and storage dynamics, through methods such as ecohydrological modelling, are necessary to understand the linkages of water in vegetation and soil at high spatial and temporal resolution. We used a large in situ dataset under grass and willow trees, including soil moisture, energy balance, water stable isotopes, and biomass accumulation to test the capability of using a tracer-aided ecohydrological model at high temporal resolution to constrain water, energy, and isotope mass balance throughout a growing season. The model captured event and seasonal dynamics of soil moisture (MAE = 0.02–0.03 m³ m⁻³), soil and xylem water δ²H values (MAE = 4.2 ‰–8.6 ‰), latent and sensible heat (MAE = 1.8–20 W m⁻²), and biomass (MAE of stem growth and leaf growth of 1.5–2.2 µm and 0.5 m² m⁻², respectively). In situ soil and vegetation isotopes were indispensable in calibration for simulating water storage, sources and fluxes of vegetation uptake, mixing processes, and water age estimates; and the importance of such datasets in informing modelling approaches was demonstrated. Modelling sap flow of the willows revealed significant water usage from neighbouring cells and the heterogeneity of root water and xylem water isotope sources. Distance-based isotopic mixing of root uptake revealed an improved estimation of xylem water stable isotopes and showed the capability of the model to reproduce the diurnal isotopic variability in spatially heterogeneous soils as a lag in the xylem water isotopes. The results additionally pointed towards further model development needs for modelling vegetation mixing in natural environments. Such numerical modelling approaches, with a physical basis and the capability of accurate simulation of multiple, inter-related variables estimation, have a high potential for the further exploration of critical zone water cycling and improved understanding of spatiotemporal changes in water availability due to vegetation–soil interactions. Continuation of integrated modelling approaches using the leverage provided by in situ data will aid future ecohydrological investigations in constraining and informing modelling while providing high-spatiotemporal-resolution insights into ecohydrological processes.