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

The partitioning of water fluxes in the critical zone is of great interest due to the implications for understanding 10 water cycling and quantifying water availability for various ecosystem services. We used the tracer-aided ecohydrological model EcH2O-iso to evaluate water, energy, water stable isotope, 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 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 evapotranspiration water use, with limited percolation of 15 summer precipitation to deeper soil layers due to the dominance of shallow root-uptake (>80% in the upper 10cm). Lower evapotranspiration under grass 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 similar to under willows (15-17 days). Considering potential xylem transit times showed a large improvement in the model’s capability to estimate xylem isotopic composition and water age, and revealed the high value of 20 in-situ data within modelling. Root-uptake was predominately derived from summer precipitation events (56%) and had an average age of 35 days, with xylem transport times taking at least 6.2 – 8.1 days. By evaluating water partitioning, energy and isotope mass-balance, along with biomass allocation, the model revealed multifaceted capabilities for assessing water cycling within the critical zone at high temporal resolution, including xylem water sources and transport, which are all necessary for short and long-term assessment of water availability for plant growth. 25

Physically-based modelling approaches, such as ecohydrological models, have a wide range of applicability, including the estimation of fluxes and storages, flux partitioning, and biomass accumulation (Asbjornsen et al., 2011). Recent 65 developments linking tracers into ecohydrological models (e.g. Kuppel et al. (2018a)) have further expanded the 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 measurement limitations of isotopic data 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 70 few approaches have utilized physically-based models to estimate storage-flux-age dynamics while considering 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 calibrate a tracer-aided ecohydrological model and estimate flux and uptake dynamics. Importantly, we also 75 used destructive methods to validate the in-situ data (Landgraf et al., 2021). The in-situ location, which also monitored soil moisture and vegetation growth dynamics, comprises a small stand of willow trees in a situation typical of NE Germany. The model, EcH 2 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 evaluate the soil storage-flux-age and vegetation uptake dynamics during the growing season by calibration across multiple data sets. This 80 primary goal was achieved by: 1) assessing the partitioning of water usage and biomass dynamics and considering how well the tracer-aided ecohydrological model could reproduce in-situ measurements of fluxes, storages, and stable water isotopes; 2) exploring potential mixing dynamics of vegetation xylem water using sap flow and xylem isotopic composition; and 3) evaluating soil water ages and transit times of vegetation water. Such 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-85 temporal resolution, and contribute to an evidence base for management strategies in water sensitive areas.

Study site
The study site is in a peri-urban area, situated in the grounds of the Leibniz-Institute of Freshwater Ecology and Inland Fisheries (IGB) in south-east Berlin (Fig 1b). The site is situated near Lake Muggelsee (80m north of the lake edge), 90 encircled in the North and West by buildings (40m), and east and south (30 and 20m, respectively) by additional vegetation.
A stream draining nearby fish retention 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.2m below the ground surface with limited annual variability (<0.1m). A previous study at the site has excluded groundwater and the nearby stream water as likely sources of water to the trees (Landgraf et al., 2021).

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&d). Willows are a riparian species, which is well adapted to moderate and high moisture typical of riparian areas or with high groundwater levels (Isebrands and Richardson, 2014). Additional selection criteria for 105 the study site is provided in Landgraf et al. (2021). The willows have a similar age of ~15 years with willow in the North (Willow 1, Fig 1c&d) being slightly larger (~9m height and stem diameter: 398mm, 16.07.2020) and the southern willow (Willow 2, Fig 1c&d) being slightly smaller (~8m height and stem diameter: 353mm, 16.07.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 south of the willows (Fig. 1d). 110 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 https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License. a previous building was demolished and the site backfilled with sandy brown earth topsoil. Soil samples taken to 1m depth reveal a relatively uniform soil structure, with a slightly higher organic horizon in the near-surface soils, which is more developed below the open grassland area south of the willows (Site B, Fig. 1c&d). 115

Climate
The climate is continental with a maritime influence (Köppen Index) and experiences substantial interannual variability in precipitation ( Table 1). The range of annual precipitation from within the climate normals  is 80% of the long term average precipitation (526 mm/year, 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). 120 The growing season is less humid and is warmer than annual averages (Table 1), with large sub-daily variability 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).

In-situ measurements
The study site was set up to continuously measure hydroclimate, soil moisture, and vegetation productivity. 130 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), collection of precipitation samples for isotope analysis, humidity, air temperature, wind speed, short wave radiation, and latent heat. Climate data from the weather station 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 Muggelsee (2013), and surrounding DWD weather 135 stations (DWD, 2021). Precipitation samples were collected at 4-hourly intervals. Each sample bottle was filled with 0.5cm 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). water content reflectometers (CS616, Campbell Scientific, Inc. Logan, UT, USA) and thermistors (BetaTherm 100K6A1IA, Campbell Scientific, Inc. Logan, UT, USA). Soil isotope samples were collected using both destructive bulk soil sampling and with an in-situ vapour analyser. Installation of polypropylene membranes (7cm long) for soil water vapour extraction at three depths (10, 40, and 100cm) was conducted at the end of May 2020 (May 20, 2020). Extracted in-situ soil vapour samples were analysed every two hours using a Picarro L2130-I cavity ring down laser spectrometer (Picarro, Inc., Santa 145 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 instrument, ICT International, Australia), dendrometers (DR, Radius Dendrometer, Ecomaticl, Dachau, Germany), and xylem water vapour (Picarro, Inc., Santa Clara, CA, USA). Polypropylene membranes were installed in two 150 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 on each willow tree, at 30cm and 170cm, and analysed within the same analyser as 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. Analyzed xylem water isotopes were evaluated for wounding effects (in early periods after installation) and any were removed from the analysis. 155 Heated cables were installed with the tubing for all vapour extraction to avoid condensation effects and modulate vapour temperatures. Additionally, tubing was flushed for 10 minutes 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 2 H and 18 O, and nonlinear regression of daily average xylem and soil isotopes with temperature (air and soil, respectively) was conducted to determine the strength of the relationship for water concentration and temperature-dependent offset. 160 Corrected in-situ xylem and soil isotopes for water concentration and temperature were consistent with bulk soil isotope and twig cryogenically extracted isotope samples (see Landgraf et al. (2021) for more details).

EcH 2 O-iso model
The EcH 2 O-iso model is a physically-based, distributed, tracer-aided model coupling vegetation, and soil energy and water balance, and carbon utilization (Maneta and Silverman, 2013). The isotope module (simulating δ 2 H, δ 18 O, and water age) 165 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 conceptualisation of energy, water, and isotope balances (conceptual diagram, Fig. S1). Further details of the model are provided in Maneta and Silverman (2013) and Kuppel et al. (2018a).

EcH 2 O-iso energy balance 170
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 175 using a Jarvis-type stomatal conductance model dependent on soil moisture, vapour pressure deficit, air temperature, and incoming radiation. Plant hydraulics were incorporated using the Soil Plant Atmosphere Continuum module (SPAC) which tracks leaf water potential and conductivity limitations due to cavitation. The SPAC module includes supply-demand functions in the rhizosphere (as a function of soil hydraulic conductivity, root area index, and pore-disconnectedness index), and stem and leaf (function of vegetation conductivity and leaf area index). The supply-demand further regulates 180 transpiration by soil and vapour pressure deficits. To account for potential root-uptake from outside model cells containing vegetation, estimated radial rooting lengths (Section 3.3.1) were used to estimate the proportion of water used from adjacent 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 185 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).

EcH 2 O-iso water balance
Similar to the energy balance, the water balance is estimated with a top-down approach including canopy, surface, and subsurface storage (Fig. S1). Canopy storage is estimated using a linear bucket approach, with maximum storage limited by a 190 storage parameter (CWS max ) and canopy interception driving draw-down 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 (Brooks-Corey, air-entry pressure, and vertical hydraulic conductivity parameters, λ 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. 195 Very dry soils at the study site prompted the introduction of 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 1cm 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 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.

EcH 2 O-iso isotope mixing and water ages 205
EcH 2 O-iso estimates isotopes (δ 2 H and δ 18 O) and water age for each storage using a complete mixing assumption (fully mixed at the end of each time-step). In storage, mixing is conducted with amount weighted-averaging of isotopic compositions (or age) with incoming fluxes. Out-fluxes have the same isotopic composition (and age) as the mixed water in storage. Evaporative fractionation is estimated 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 210 and Pielke (1992), with the kinetic fractionation factor modified for use in soils (Braud et al., 2005). Isotopic composition and water ages 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). The amount-weighting assumes instant mixing of soil water within the vegetation and encompasses no mixing of previous soil isotopes or water ages. At the end of each time-step the age of water in storage is advanced by one time-step. 215 To evaluate the distribution of water ages, calibrated simulations were re-run for each time-step that included a precipitation event (925 hours with precipitation over the study period). For each of the 100 'best' simulations, the input precipitation concentration was changed in a stepwise manner with a concentration of "1" on the precipitation event evaluated. For example, the 25 th precipitation event (i.e. 25 th hour with precipitation, here, January 7 th , 11:00am), the precipitation concentration was set to 1 while all other precipitation and initial storages had concentration of 0. The concentration was 220 tracked through the domain to provide a proportion of water in each storage originating from the 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. Knighton et al. (2020) showed that improved estimation of xylem isotopes could result when considering potential 225 vegetation storage, and mixing of soil water older than one time-step within EcH 2 O-iso. Similarly, Seeger and Weiler (2021) showed promising conceptualization of isotope mixing within vegetation based on a convolution approach using in-situ measurements. To evaluate the dynamics of transit times in xylem during the study period, a combination of modelled results (sap flow and soil isotopes) and assumed vegetation root distributions were used to estimate vegetation xylem water composition and age at the average measurement height (1m) while accounting for spatial variability of soil water isotopes 230 (laterally and with depth). The tree water mixing routine (applied post-EcH 2 O-iso calibration) was developed with the assumption that isotopic mixing is dependent on the rooting distance (cf. Seeger and Weiler (2021)) and using the convolution equation to mix water from different pools (e.g. outside model cell or different soil layers). This differs from the https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License.

Tree water mixing and transit times
traditional approach currently utilized in EcH 2 O-iso which instantly mixes uptake-water throughout the tree. Hereafter, the new xylem mixing will be referred to as distance-based mixing, and the method used internally by EcH 2 O-iso is referred to 235 as instantaneous mixing.

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 using 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 and SPAC module within EcH 2 O-iso. 240 Following EcH 2 O-iso, the rooting distribution (k root ) was used to determine the proportion of total roots in each soil layer (Kuppel et al., 2018b). Vertical distances (v) from the base of the vegetation to the vertical centre of the biomass in each layer were estimated using a log function of rooting proportions (Sperry et al., 2016): where i is the soil layer, p is the cumulative proportion of roots (from the surface) and β is a rooting distribution parameter 245 ( = −0.0089 * + 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): where Vol is the volume of roots in layer 1, d 1 is the depth of layer 1, D is the total depth of the soil, and a is an aspect parameter, controlling radial distance. The radial distance is then estimated as: 250 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.

Xylem transit time, and calibration 255
Since much of the root systems below trees are not well characterized, including limited information on the proportion or total lengths of fine roots to transport roots are available, and 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 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 indirectly accounted 260 (averaged) across different root diameters. Root-uptake volume weighted (estimated from EcH 2 O-iso, Fig. 2c) convolution equations (gamma distribution) were used to estimate xylem δ 2 H from each soil layer for each hourly time-step (Fig. 2d).
The approach introduces a scale parameter (β, hours) for each soil layer and one shape parameter (α) as a function of total https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License. root length (four total parameters) (Fig. 2b). This approach assumes that there is no cavitation (no drought conditions) within the xylem or roots of the plant. 265 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-hour intervals. Measured and simulated xylem isotopes were 270 averaged over the respective intervals prior to evaluation with the Kling-Gupta efficiency (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 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 isotopes and sap flow. The use of measured soil isotopes and sap flow tests the maximum potential for 275 how each model performs and is not limited to the performance of EcH 2 O-iso for sap flow or soil isotope. Further, the twin calibration approach shows how the error in simulation soil 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 multi-criteria EcH 2 O-iso calibration (Section 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 (4 additional parameters). 280 The AIC is estimated using the log-likelihood of each model fit and number of parameters, where lower AIC indicates better model performance.

Model set-up
To best account for measurements in both willows and the soil moisture measurement locations, the model was set up with 285 three square grids (6m). 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 January 1, 2020, and October 31, 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 100cm) for all model cells. Further, previous analysis suggests that primary water sources of the willow trees were within the upper 100cm of the soil, with little isotopic evidence of trees using water with distinct groundwater or stream water signatures, 290 despite their proximity (Landgraf et al., 2021). The low local measured soil moisture and high transpiration in the willows suggests notable water use from outside the measurement location (model cell), and therefore parameterization of rooting radius (aspect parameter) was set to result in water being available from outside the cell. Due to the notable differences in shallow soil moisture between Site A and B (organic content), soil parameters were different below the willows and the open grass (see calibrated parameters, Table S1). Initial soil moisture was set to field capacity which was revealed by model 295 testing to have negligible influence on soil moisture simulations during the growing season. Soil isotopes were initialized using the average measured isotopes. Initial soil isotopes for the shallowest soils were not sensitive 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.

Model forcing data, calibration and evaluation 300
Hourly model forcing data primarily consisted of data measured with the mobile weather station at the study site (Section 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).
To best constrain the model performance, and provide an evaluation of how well the model estimates the dynamics of water, tracer, and energy fluxes during the growing season, all available data (except δ 18 O) were used to constrain model 305 parameterization. To simplify the presentation of results, we present δ 2 H as the results were not greatly different from δ 18 O.
Furthermore, the evaluation of a single growing season at the site limits the feasibility of a split calibration approach. To simultaneously evaluate the variability, bias, and correlation of the measured and simulated datasets, the KGE was used for multicriteria calibration of all datasets.
Step-wise calibration was utilized to calibrate tracer, energy, and water balance, Infeasible parameter combinations 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. 315

EcH 2 O-iso soil water, isotope and energy balance evaluation
Model calibration provided adequate simulations of each calibration variable. Similar to measured data, simulations of soil 320 moisture and soil isotopes showed decreasing dynamics with increasing depth from the surface (Fig. 3). Despite the close proximity, Site A and B had largely different soil moisture and isotopic response in the upper soils ( Fig. 3a-d). Simulations were able to capture the drier soil moisture in Site A relative to Site B (Fig. 3a&c); however, the simulated dynamics of average moisture in layer 1 were greater than the measured moisture at 10cm likely due to greater moisture retention in the organic material near the surface preventing infiltration reaching 10cm. Sub-discretization of the shallow layer 1 soil using 325 calibrated parameters revealed that while average moisture in the upper 10cm is high (blue lines, Fig. 3a), percolation of infiltrated water down to 10cm during the summer months is limited (Fig. 3a, red line). Late season soil moistures (average and at 10cm) were over-estimated relative to measurements, coinciding with the underestimation of sap flux during the same period (Fig. 4).
The general dynamics of the soil isotopes at 10cm in Site A was reasonably represented, however the larger variability 330 during a period of depletion in August was not captured. Soil moisture and isotope dynamics in the deeper soils (40cm) were much more damped than 10cm, with a slight under-estimation bias in 40cm at Site B despite appropriate dynamics (Fig. 3d).
The much lower variability of isotopes and soil moisture at 100cm was similarly captured by the model calibration. higher depletion simulated at Site A relative to measured was consistent with the deeper (100cm) soil isotopes measured at Site B. Outliers of in-situ soil isotopes towards the end of the growing season (end of August and September, Fig. 3c-f) were 335 a result of rapid temperature changes which were too marked for the system (i.e. heated cables) to control and caused temporary condensation effects in the tubing. These data were few and did not have notable effect on the soil isotope investigation. Estimated energy balance (latent heat, sensible heat, transpiration flux, and temperature) were additionally shown to match in-situ measurements quite well (Fig. 4). Despite only utilizing latent heat and surface temperature at Site B in the calibration, simulated sensible heat and soil temperature at multiple depths were additionally reasonably captured. Latent during the early simulation period with increasing soil depth is likely due to propagation of a slight over-estimation of early-350 growing season surface temperature at each site (Fig. 4) potentially due to lower estimated leaf area index at the beginning of the growing season (Fig. 5).

Growing season dynamicsfluxes and biomass accumulation
Simulated biomass for foliage and tree diameter produced reasonable results for the two willow trees throughout the simulation period, with only minor deviations from both measured and remote sensed (MODIS) datasets (Fig. 5a&b). Average simulations (solid lines) under-estimated LAI at the beginning of the growing season; however, these difference 360 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 (Fig. 5c). The dynamics and total net stem growth was adequately captured by the model for each willow; however, there was a slight underestimation of average net growth in Willow 2 (Fig. 5c&d). Modelled partitioning of water fluxes and biomass allocation in the willow and grass sites showed marked differences despite the relative proximity. Evapotranspiration (ET) dominated water partitioning throughout the summer, with the greatest ET between June 15 and September 15. Total ET was dominated by transpiration (T r ) particularly for the willows 370 (>70%). This coincided with a decreased interception evaporation (E i ) and soil evaporation (E s ) proportions 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 difference in infiltration between the willow and grass sites. The high root uptake rate of the willow resulted in negligible percolation of water to deeper soil layers during the growing season. Continuation of drier conditions 375 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 slight shift of root growth in the pre-growing season to foliage growth during 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 380 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 https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License.

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-water mixing was conducted in separate calibrations using both simulated and measured soil isotopes. Instantaneous mixing of root-uptake water throughout the tree (instant uniformity in all xylem at each time-step, as simulated by the EcH 2 O-iso model structure) was able to capture the general dynamics of xylem water in both willows (Fig.  395 6a&b), though simulated day-to-day variability was limited. Unsurprisingly, the instantaneous mixing approach showed improved performance when measured soil isotopes were used rather than simulated soil isotopes, increasing seasonal and day-to-day variability (Table 4; Fig. 6c&d). Distance-based mixing simulations of xylem water at 1m were significantly better than the instant mixing simulations (lower AIC for distance-based mixing). Similar to the instant mixing simulations, the use of measured soil isotopes within the distance-based mixing approach showed improved results (AIC) compared to 400 the use of simulated soil isotopes with the exception of the 24-hour time-step in Willow 2 ( Table 4). The decrease in KGE in Willow 1 from simulated to measured data (Table 4) was an artefact of the shorter measured 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).

Figure 6: Simulated xylem deuterium using simulated soil deuterium with distance-based and instant mixing in (a) Willow 1, and (b) Willow 2. Simulated xylem deuterium using measured soil deuterium with distance-based and instant mixing in (c) Willow 1, and (d) Willow 2. For each simulation, standardized sub-daily δ 2 H (S(δ 2 H)) shows the measured and simulated soil and xylem average hourly variability.
Xylem simulations with instant-mixing and modelled soil isotopic composition were unable to capture the observed 410 standardized hourly variability of the xylem isotopic composition (Fig. 6a&b). This is consistent with the limited standardized hourly variability of simulated soil isotopes relative to measured soil isotopes (Fig. 6a&b). Distance-based simulations with modelled soil isotopes could reasonably capture the timing of sub-daily dynamics, due to the sub-daily variability of sap flow volumes. Xylem water mixing estimations using measured soil isotopes revealed a reasonable standardized hourly variability of xylem isotopes (Fig. 6c&d). Instant mixing unsurprisingly showed dynamics similar to soil 415 water isotopes, which peaked before the standardized measured xylem water. Distance-based mixing with measured soil 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 isotopes at 10cm 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). The increase in model performance is consistent with the reduced variability of measured xylem with 420 https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License.
increasing time-steps, increasing in similarity to the variability observed in soil isotopes. Likewise, the KGE increased most rapidly with time-step on average when simulated soil isotopes were used in mixing as they lower variability than measured soil isotopes (Fig. 3).

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&b). Average layer 1 water ages below the willow (15 ± 1 days) were similar to the grass 430 (17 ±1 days), with larger proportions of 7-day (36±2%), 14-day (59±2%), and 30-day (84±1%) relative to the below grass (31±1, 52±1, and 67±1%, respectively). The higher proportions of young water below the willow were primarily the result of lower moisture, 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 waters. Percolation throughout the summer under the grasses resulted in a stabilization of water age (average: 73 ± 2 days), with 81±1% of water older than 30 days. Average root 435 water age in the root-tip was a mixture of water ages for all soil layers, with the high root uptake from shallow soils willow resulting in a younger average uptake water age in the willow roots (28 ± 6 days) compared to the more equal mixture of shallow and mid-depth soils uptaken by the grass roots resulting in older water uptake (51 ± 3 days) (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 440 10cm) resulted in low percentages of late winter and spring water (Pre-April 1) in the early growing season (13 and 18% below willow and grass, respectively, Fig. S3, Table S2&S3), which were negligible at the end of the growing season (0 % below willow and grass, Table S2&S3). 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 445 distribution and the temporal dynamics of infiltration resulted in differences of spring and winter water usage for the grass https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License. and willow throughout the growing season, with 33, 16, 9% for early, mid-, and late growing season willow (57, 35, 24%, respectively for grass) root-uptake of spring/winter water.

(40cm) in (c) willow 1 and (d) grass, 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 1m using distance-based mixing using simulated or measured soil isotopes.
Since instantaneous mixing has an increase of xylem water age from the roots of zero days (reaches 1m instantly), xylem ages are only shown for the distance-based mixing approach (using simulated and measured soil isotopes). Using simulated soil isotopes, distance-based mixing revealed an average transit time (all time-steps) of 149 ± 30 hours from the roots to 1m 455 height (Fig. 7g). Distance-based xylem mixing using measured soil isotopes revealed a mean transit time of 237 ± 97 hours ( Fig. 7g). Smaller uncertainty of transit time estimated from simulated soil and sap flow was a direct result of smaller temporal variability (lower degrees of freedom) of simulated soil isotopes compared to measured soil isotopes. 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 acceptable reproduction of ecophysiological processes. At fine spatial scales, such as plot sites, reproducing flux-storagebiomass dynamics and feedback mechanisms are important as they can accentuate large-scale nonlinearities (Asbjornsen et al., 2011). Continued advancements of complex (e.g. physically-based) models and fusion with temporal high resolution, 465 multiple 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 in 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 intensive in-situ monitoring site provided a means to directly evaluate how well the model could simultaneously reproduce flux-storagetracer-biomass dynamics through multiple criteria throughout high variability of a growing season. 470 Overall, the model performed well against all measured data in both the willow and neighboring grass sites, including soil water storage, water and energy fluxes, tracer 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;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 (Coenders-Gerrits et al., 475 2014;Dubbert et al., 2014;Schlesinger and Jasechko, 2014;Zhou et al., 2016), including lower Tr/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 10cm, transpiration in both willows and grass were maintained throughout the study period (Fig. 4), with only minor under-estimation of the transpiration in the willows toward the end of the growing season. This was primarily possible within EcH 2 O-iso due to the capability of the willows 480 using water from adjacent model cells. As observed in other nearby studies, soil evaporation contribution to ET under both the willow and grass was damped throughout the growing season due to leaf shading (LAI) causing decreased energy availability at the surface Smith et al., 2021). Although willows are generally a riparian species accessing groundwater or stream and lake water, the willows at this site did not reveal notable contribution from the deeper (>2m) groundwater. The predominantly shallow root-uptake contribution (driven by higher near-surface water availability 485 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 the expected shallow rooting zone. The deeper rooting distribution in the grass is likely due to high uptake competition in the near surface soils with the nearby high water-490 use willows (Moustakas et al., 2013). Competition potentially drove deeper rooting of the grass, causing the more dynamic https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License. soil moisture at 40cm as measured below the grass. The slight decrease in root biomass production in the willows during the drier summer months suggested that vegetation was not under water stress in the study period (Eziz et al., 2017). Average hourly simulated transpiration dynamics were comparable to sap flow measurements (Fig. S4), including small, but notable uptake during the night, similar to observations in other woody plants (Dawson et al., 2007). 495 Within this study, the relatively low soil moisture did not propagate to a notable decrease in sap flow despite a high proportion of the rooting zones in the shallow soils. While the low moisture contrasts high water usage and shallow roots of the willows, in-situ isotope dynamics were essential in confirming measurements of shallow moisture were not underrepresenting infiltration events. Multicriteria calibration using high-temporal (and highly dynamic) resolution soil isotopes revealed a strong positive correlation of simulated shallow soil moisture KGE and simulated shallow soil isotope 500 KGE (i.e. model performance of soil isotopes only improved with performance of soil moisture). Bulk soil isotope samples, with a coarser temporal resolution, would likely have been insufficient to adequately capture the necessary isotope dynamics.
EcH 2 O-iso was only able to adequately estimate soil moisture at 10cm through sub-discretization of the soils above 10cm, which suggested a high retention of water in soils above the measurement depth (e.g. 0-5cm), up-taken by roots or soil evaporation before there is sufficient water to percolate to 10cm. These small-scale spatial variations, as well as the large 505 spatial differences from 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 structure can impact throughfall (Dalsgaard, 2007;Gerrits et al., 2010;Li and Liang, 2019), and if such data were available, it could likely further improve soil moisture simulations. In this regard it is possible that zones of higher soil moisture are 510 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.

Evaluation of mixing dynamics of root-uptake and implications of the rooting zone
We used output from EcH 2 O-iso (i.e. the root-uptake sources) and a simple, parsimonious approach to further explain the time-delay, and variability of xylem water from soil water. The approach showed significant improvements over the use of 515 the instantaneous mixing currently applied within most modelling approaches. The inclusion and calibration of rooting radial distance into EcH 2 O-iso showed comparable results to other willows (0.67 v. 0.5 tree height/diameter of root spread ;Phillips et al. (2014)) which increases confidence in 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 isotope composition underlines the potential influence of spatio-temporal variability of source waters on xylem 520 isotopes, as suggested by De Deurwaerder et al. (2020) and Seeger and Weiler (2021).
As with partitioning soil moisture and vegetation water usage, high-resolution in-situ isotopes (xylem) were indispensable for constraining the mixing dynamics of root-uptake. Coarser destructive xylem sampling could potentially over-estimate xylem variability due to the large sub-daily variability, damping seasonal patterns (e.g. enrichment in September, Fig. 6 influencing estimated root-uptake profiles. Further, while the use of tracer-aided ecohydrological model output for 525 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 (Fig. 6), seasonal magnitudes of xylem isotope dynamics were predominantly due to differences in simulated v. measured soil isotopes in the shallow soils (Fig. 3 & RU in Table 3). This additionally suggests the great 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 530 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 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). Uncertainty in isotopic measurements, complex physiological processes within the trees which may cause fractionation, limitations of the model structure and interactions between these were all potentially limiting factors for reducing the xylem mixing model 535 performance. In-situ measurement uncertainty, particularly in xylem, may be large for δ 2 H (5-10 ‰; Beyer et al. (2020)), greater than the standard deviation of the daily mean δ 2 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 δ 2 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 540 von Freyberg et al., 2020). In this study, condensation uncertainties were minimized by heating the tubes and were largely restricted to the autumn (September and October) when night time temperatures could be much cooler (Landgraf et al., 2021). Further, while fractionation during uptake or within vegetation has frequently been discussed and reexamined recently (Brinkmann et al., 2019;Poca et al., 2019;Vargas et al., 2017), the vegetation isotopes reflect the soil isotope sources ( Fig. S5; Landgraf et al. (2021)) which suggests that vegetation specific fractionation is not likely a major 545 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 (Sperry et al., 2016)) 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., 550 2014). While such a release of water may be present during peak transpiration hours, as suggested by the hourly variability of basal diameter (Fig. S6), the high water usage of willows potentially limits the proportional contribution of cell storage to xylem water.

Implications of water ages in soil water and xylem fluxes
Quantifying the age of water has significant implications for understanding how water is cycled through the critical zone, 555 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 https://doi.org/10.5194/bg-2021-278 Preprint. Discussion started: 4 November 2021 c Author(s) 2021. CC BY 4.0 License.
how 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 . 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 560 (volume-weighted). This differs from other studies (e.g. Allen et al. (2019)) which have shown a higher usage of winter precipitation in leafy vegetation (beech and oak) but is consistent with estimated high contribution of recent precipitation during the growing season (Miguez-Macho and Fan, 2021). However, differences in hydroclimate, as well as variations in rooting depths between species and 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 565 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))implies that for the conditions present at the study site, inclusion of preferential flow or variable pore size would 570 not influence water ages for the uptake of either in 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 estimated in this study were of a similar order of magnitude to estimates for other woody vegetation; averaging between 6.2 and 8.1 days with a maximum of 10.3 days at our site compared to between two days and one month in other studies (Brandes et al., 2007;Meinzer et al., 2006). While our estimates were 575 towards the lower transit times compared to most other studies, the willows at our site were younger, had smaller diameters, and had higher xylem sap flow velocity, consistent with high water usage and low stomatal control (xylem cavitation control, Wikberg and Ögren (2007)). High-resolution in-situ measurements, capturing soil and xylem isotope dynamics were essential 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 580 release can have on xylem isotopes (Section 5.2) (Brandes et al., 2007;De Deurwaerder et al., 2020); there are likely implications for water ages and 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 utilize 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 585 water. Due to the relatively young uptake water ages of willows in this study, the inclusion of any cell water storage may significantly change the mean age of xylem water, and would increase the time period for mixing in modelling approaches, such as the one implemented here. Given the limited studies describing the ages of water and tracers in transpiring trees (Sprenger et al., 2019), more studies quantifying the transport through the xylem, as conducted here, and inclusion of additional mixing within vegetation (e.g. Steppe et al. (2006)) would be beneficial to further constrain plant water use 590 estimations.

Conclusion
In order to increase understanding of the mechanisms of water cycling within 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 595 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 temporal high resolution to constrain water, energy, and isotope mass balance throughout a summer growing season. The model captured event and seasonal dynamics of soil moisture, soil and xylem water deuterium, latent and sensible heat, and biomass (stem growth). In-situ soil and vegetation isotopes were 600 indispensable for simulating water storage, vegetation water use and source and mixing and water age estimates. Water usage at the willow plot sites showed higher evapotranspiration than grass, rapidly drying the soils down to 1m with most root-uptake from upper 10cm. The grasses showed a greater proportion of soil evaporation in evapotranspiration compared to the willows, and root-uptake equally mixed from the upper 40cm. Both vegetation units showed most biomass allocated to foliage growth throughout the growing season. Water age and isotopic mixing of root-uptake showed an improved 605 estimation of xylem water stable isotopes, and the capability of the model to estimate the isotopic lag from soils in xylem.
Estimated willow root-uptake showed that summer precipitation was the dominant source (56%), with an average uptake age of 28 days (near-surface sourced water), and xylem transport time from root to 1m above the surface between 6.2 -8.1 days.
Such numerical modelling approaches, with a physical basis and the capability of accurate multifaceted variable estimation, have a high potential for further exploration of critical zone water cycling and improved understanding of spatio-temporal 610 changes in water availability due to vegetation-soil interactions. Continuation of integrated modelling approaches using the great value provided by in-situ data will aid future ecohydrological investigations in constraining and informing modelling, while providing high spatio-temporal resolution information of ecohydrological processes.