3.1.1 Soil profiles

Figure 2a and b show a very stable soil water content profile and a more variable δ_soil profile from DaS 166 at 15:45 LT to DaS 168 at 05:00 LT. Soil was dry at the surface (0.058 m³ m⁻³ < θ < 0.092 m³ m⁻³ for the layer from 0.015 to 0.040 m), whereas it was closer to saturation at a depth of −1.30 m (θ=0.34 m³ m⁻³ ±0.012 m³ m⁻³, and the estimated θ_sat=0.40 m³ m⁻³; see Appendix A). According to the measured soil matric potentials (Fig. 2c), soil water was virtually unavailable ($\le -\mathrm{1.5}$ MPa) above a depth of −0.5 m. Soil moisture remained unchanged in the top 25 cm during the sampling period ($\mathit{\theta }=\mathrm{0.08}\pm \mathrm{0.00}$ m³ m⁻³) as well as at −1.30 m from DaS 166 at 15:45 LT to DaS 168 at 05:00 LT ($\mathit{\theta }=\mathrm{0.33}\pm \mathrm{0.01}$ m³ m⁻³), showing that roots were predominantly extracting water from deep soil layers.

Figure 2(a) Measured soil volumetric water content (θ), (b) oxygen isotopic composition (δ_soil), and (c) calculated soil matric potential (ψ_soil) profiles during the sampling period.

Download

Water in the top soil layers (−0.040 m < z < −0.015 m) was isotopically enriched (−3.2 ‰ < δ_soil < 0.3 ‰) in contrast to the deepest layer (${\mathit{\delta }}_{\mathrm{soil}}=-\mathrm{7.34}$ ‰ ±0.30 ‰ at −1.30 m). Following the labeling of the reservoir water on DaS 166 at 17:00 LT, δ_soil reached a value of 36.9 ‰ at −1.50 m on DaS 167 at 17:00 LT. The development of the vegetation on DaS 166–168 (LAI =5.6) and the observed surface θ values lead us to assume that the rhizotron water losses were solely due to transpiration flux (i.e., evapotranspiration equals transpiration). The soil water oxygen isotopic exponential-shaped profiles were due to fractionating evaporation flux and, to a great extent, the fact that the soil was bare or the tall fescue cover was not fully developed in the early stages of the experiment. Therefore, the differences in the soil water oxygen isotopic profile observed on the four different sampling dates were either due to lateral heterogeneity (e.g., upper soil layers), the soil capillary rise of labeled water from the reservoir (deep soil layers), or the hydraulic redistribution of water through roots (if the isotopic composition of the redistributed water differs from that of the soil water at the release location). We noted an isotopic enrichment of 1.0 ‰ of soil water observed on DaS 168 at 05:00 LT at a depth of −0.9 m with respect to the mean δ_soil value across previous sampling dates. This could partly be due to factors such as the upward preferential flow of labeled water from the bottom soil layers and could, therefore, be a sign of the lateral heterogeneity of the soil. Another reason for this would be the hydraulic redistribution of labeled water by the roots. However, it was not possible to evaluate the relative importance of these three processes (lateral heterogeneity, capillary rise/preferential flow, and hydraulic redistribution) in the setting of the soil water isotopic profile, as the physically based soil–root model presented in section 2.4 does not account for soil liquid and vapor flow. This was also not the primary intent of the present study.

The observed RLD profile (Fig. 1a) showed a typical exponential shape, i.e., a maximum at the surface (5.42±0.34 cm cm⁻³) down to a minimum at −1.10 m (0.540±0.35 cm cm⁻³), and it increased again from the latter depth up to a value of 1.660 cm cm⁻³ at −1.30 m. This significant trend was most probably a direct consequence of the high soil water content value in this deeper layer.

3.1.2 Plant water and isotopic temporal dynamics

The temporal variation in δ_tiller (Fig. 3a) was found to be either (i) moderate during day and night, i.e., from DaS 167 at 06:00 to 11:00 LT (${\mathit{\delta }}_{\mathrm{tiller}}=-\mathrm{2.6}\pm \mathrm{1.4}$ ‰) and from DaS 167 at 21:30 LT to DaS 168 at 00:00 LT (${\mathit{\delta }}_{\mathrm{tiller}}=-\mathrm{2.7}\pm \mathrm{0.4}$ ‰), (ii) strong during the day, i.e., from DaS 167 at 11:00 to 18:00 LT (maximum value of 20.9 ‰ on DaS 167 at 12:40 LT), or (iii) strong during the night, i.e., from DaS 167 at 04:00 to 06:00 LT (maximum value of 36.4 ‰ on DaS 167 at 05:15 LT) and from DaS 168 at 00:00 to 06:00 LT (maximum value of 14.6 ‰ on DaS 168 at 04:00 LT). Note that transpiration (Fig. 3b) also occurred at night during the sampling period, due to the relatively high temperature in the glasshouse leading to a value of atmospheric relative humidity smaller than 85 % (Fig. 3b). From 12:00 to 14:00 LT and from 16:00 to 17:00 LT on DaS 167 (case ii), high values of leaf transpiration corresponded to high values of δ_tiller.

Figure 3(a) Time series of tiller and leaf water oxygen isotopic compositions (δ_tiller and δ_leaf, respectively, ‰). (b) Transpiration flux (T, in m d⁻¹), relative humidity (HR, %), and leaf water potential (ψ_leaf, in MPa) from days after seeding (DaS) 167 at 04:00 LT to DaS 168 at 11:00 LT. Labeling was carried out on DaS 166 at 17:00 LT.

Download

3.1.3 Partial decorrelation between water and isotopic state variables

Figure 4 shows that variables describing plant water status, i.e., T and RH (Fig. 4a) and T and ψ_leaf (Fig. 4b), were well correlated: the coefficient of determination (R²) was equal to 0.78 and 0.70 for the entire experimental duration, respectively. However, linear relationships between water status and isotopic variables were either nonexistent, e.g., between T and δ_tiller (R²=0.01, Fig. 4c) and between ψ_leaf and δ_tiller (R²=0.00, Fig. 4h), or characterized by a low R² and high p value (e.g., between T and δ_leaf, R²=0.43, p > 0.05, Fig. 4d). The partial temporal disconnect between δ_leaf and T could not be attributed to problems with the isotopic methodology during processes such as the vacuum distillation of the water from the plant tillers and leaves: the water recovery rate was always greater than 99 % and Rayleigh distillation corrections (Dansgaard, 1964; Galewsky et al., 2016) were applied to standardize the observed oxygen isotopic composition values to a 100 % water recovery (based on the comparison of the sample weight loss during distillation and the mass of collected distilled water). The evolution of δ_leaf was strongly correlated with that of δ_tiller during the day (R²=0.90), whereas it was not correlated during the night (R²=0.00, Fig. 4j). These observed correlations are in agreement with the Craig and Gordon (1965) model revisited by Dongmann et al. (1974) and later by Farquhar and Cernusak (2005) and Farquhar et al. (2007). The model, which is extensively used in the current literature (e.g., Dubbert et al., 2017), states that, at isotopic steady-state, δ_leaf is a function of the input water oxygen isotopic composition (δ_tiller) among other variables, i.e., leaf temperature (not measured during the experiment), stomatal and boundary layer conductances, oxygen isotopic composition of atmospheric water vapor, and relative humidity.

Figure 4Correlations between measured variables: oxygen isotopic compositions of xylem and leaf waters (δ_tiller and δ_leaf, respectively, in ‰), transpiration rate (T, in m d⁻¹), relative humidity (RH, %), and leaf water potential (ψ_leaf, in MPa). The coefficients of determination (R² values) are reported for all data as well as separately for “day” (gray symbols) and “night” data (black symbols; see Appendix C for definition of “day” and “night” experimental periods). Regression lines are drawn for linear models with a p value < 0.01

Download

It is generally difficult to observe a statistically significant δ_leaf−δ_tiller (Fig. 4j) relationship at this temporal scale under natural abundance conditions in the field, as the soil water isotopic weak gradient translates into weaker δ_tiller temporal dynamics. The quality of the linear fit between δ_leaf and δ_tiller data collected during the day (R²=0.90) was made possible in this specific experiment by the artificial isotopic labeling pulse that enhanced the soil water isotopic gradient, which in turn increased the range of variation in δ_tiller, ultimately highlighting the δ_leaf−δ_tiller temporal correlation. Air relative humidity is a driving variable of δ_leaf in the model of Dongmann et al. (1974) via the competing terms of (1− RH) ⋅δ_tiller and RH ⋅δ_atm, where δ_atm is the atmospheric water vapor isotopic composition inside the glasshouse. An overall significant linear correlation was observed between RH and δ_leaf during the experiment (R²=0.57, Fig. 4g). During the two night periods (i.e., from 04:00 to 06:00 LT and from 20:30 to 07:00 LT), as relative humidity increased in the glasshouse (51 % < RH < 85 %, Fig. 3b), the influence of the isotopic labeling of the tiller water (due to the labeling of deep soil water) via the (1− RH) ⋅δ_tiller term decreased to the benefit of the RH ⋅δ_atm term (with δ_atm values ranging from −15.9 to −10.7 ‰ and a mean of $-\mathrm{13.1}\pm \mathrm{1.6}$ ‰, data not shown). This was especially visible between 04:50 and 06:00 LT on DaS 167 and between 01:00 and 06:00 LT on DaS 168, when δ_tiller reached greater values than δ_leaf.

From a different perspective, as three plant water samples were pooled to reach a workable volume for the isotopic analysis at each observation time without replicates, the isotopic signal fluctuations may reflect both its temporal dynamics and its variability within the plant population.

3.2.1 Rooting depth and transpiration rate control δ_tiller and ψ_leaf fluctuations, respectively

Despite the use of a global optimizer and 4 degrees of freedom (L_pr, k_axial, k_sat, and λ, see optimal values in Table 1) that specifically aimed at matching the simulated and observed temporal dynamics of δ_tiller, none of the 60 root system groups or the average population could reproduce the measured fluctuations in time (R²=0.00, Fig. 5a), regardless of the weight attributed to this criterion in the objective function. The predicted versus the observed δ_tiller distributions including all plant groups and observation times differed noticeably but not significantly (6.6±8.4 ‰ and 3.7±8.4 ‰, respectively) when pooling three simulated δ_tiller randomly at each observation time (P > 0.01 in 92 cases out of 100 repeated drawings), as in the measurements. In addition, the simulated ψ_leaf fitted the observations well (R²=0.67, with overall distributions of $-\mathrm{0.175}\pm \mathrm{0.053}$ MPa and $-\mathrm{0.177}\pm \mathrm{0.053}$ MPa, respectively; Fig. 5c). When analyzing the distributions of ψ_leaf and δ_tiller per maximum root system depth (Fig. 5b, d), it appears that the ψ_leaf signal is not sensitive to the rooting depth (Fig. 5d), whereas δ_tiller is more sensitive to rooting depth than to the temporal evolution of the plant environment (Fig. 5b).

Figure 5Variation in δ_tiller and ψ_leaf in time and across the 60 groups of simulated root systems. (a) The measured (thick red line) and simulated (thin gray lines, one line per root system group, following a “swarm” pattern) temporal dynamics of δ_tiller. (b) Boxplot of simulated δ_tiller values for each root system maximum depth, in 1 cm increments. (c) The measured (thick green line) and simulated (thin gray lines, one line per root system group, following a “roller-coaster” pattern) temporal dynamics of ψ_leaf. (d) Boxplot of simulated ψ_leaf values for each root system maximum depth, in 1 cm increments.

Download

Table 1Optimum and limits of the four-dimensional parametric space explored by the global optimization algorithm aiming at minimizing the difference between simulated and observed δ_tiller and ψ_leaf, as well as their standard deviation from average values during the full experiment.

Download Print Version | Download XLSX

This leaves us with two hypotheses. First, the “roller-coaster hypothesis” states that δ_tiller rapidly goes up and down with all individuals onboard the same car (i.e., little variability within the population, unlike the predictions in Fig. 5a but similar to the simulated ψ_leaf in Fig. 5c). If this holds true, the physical model lacks a process that would capture the observed temporal fluctuations of δ_tiller. Second, the “swarm pattern hypothesis” states that δ_tiller is rather stable in time, but its values within the plant population are dispersed like a flying swarm; thus, δ_tiller values sampled at different times fluctuate, not due to temporal dynamics but instead owing to the fact that different individuals are sampled (Fig. 5a).

The model suggests that the tall fescue population ψ_leaf follows a “roller-coaster” dynamics that is driven by transpiration rate, whereas the population δ_tiller follows a “swarm” pattern that is driven by the maximum rooting depth of the sampled plants. As no correlation could be expected between the drivers (the maximum rooting depth of the sample plants and the canopy transpiration rate), our analysis explains the absence of a correlation between δ_tiller and ψ_leaf or transpiration rate.

In future experiments and in the specific context of labeling pulses, sampling more plants at each observation time would help disentangle the spatial from the temporal sources of variability of ψ_leaf and δ_tiller. However, it would be at the cost of the temporal resolution of observations or would necessitate a larger setup with more plants in the case of controlled conditions experiments.

3.2.2 Independent observations support the validity of the hydraulic model predictions

In the last 12 h of the experiment (DaS 167 at 17:00 LT to DaS 168 at 05:00 LT), the measured soil water content increased by 0.029 m³ m⁻³ at a depth of −0.9 m, which could be a sign of nighttime hydraulic redistribution. During the same period, the physical model predicted a cumulative water exudation sufficient to increase soil water content by 0.003 m³ m⁻³, as the soil water potential was low enough to generate reverse flow but high enough not to disrupt the hydraulic continuity between the soil and roots (Carminati and Vetterlein, 2013; Meunier et al., 2017a). While this increase is smaller than the observed water content change, it is only a component in the soil water mass balance. Given the soil water potential vertical gradient, upward soil capillary water flow may have accounted for another part of the observed moisture change. Experimental observations also show that δ_soil increased by 1.0 ‰ at a depth of −0.9 m during that time (−6.2 ‰, which was a value significantly higher than −7.1 ‰ ±0.1 ‰ at earlier times based on an ANOVA analysis, P < 0.01), whereas our simulations of hydraulic redistribution generated a 0.34 ‰ increase of δ_soil. As soil capillary flow may not generate local maxima of δ_soil (no enrichment observed at surrounding depths, see Fig. 2b), and soil evaporation is assumed negligible at that depth, it is likely that the observed local enrichment was entirely due to hydraulic redistribution, which would then be underestimated by a factor of about 3 in our simulations. Increasing water exudation by a factor 3 would imply a simulated water content change due to exudation of 0.0090 m³ m⁻³ absolute water content, which remains compatible with the experimental observation. Between a depth of −1.1 m and −0.9, the nighttime water flow pattern transitioned from exudation to uptake in both measurements and predictions. At −1.1 m, the model predicted a cumulative water uptake sufficient to decrease the soil water content by 0.0101 m³ m⁻³, compared with the observed 0.0141 m³ m⁻³ total soil water content decrease. The remaining 0.004 m³ m⁻³ water content decrease may have contributed to the recharge of the soil layers above via capillary flow, which was not simulated. Therefore, all relevant measurements (local increase in soil water content and local enrichment of water isotopic composition) and simulation results (S < 0, i.e., local water release from roots) clearly converge to the conclusion that hydraulic lift occurred in the vicinity of a depth of −0.9 m in the early morning of DaS 168.

As far as fitted parameter values are concerned, L_pr ($\mathrm{2.3}\times {\mathrm{10}}^{-\mathrm{7}}$ m MPa⁻¹ s⁻¹) was in the range found by Martre et al. (2001) for tall fescue (${\mathrm{2.210}}^{-\mathrm{7}}\pm \mathrm{0.1}$ m MPa⁻¹ s⁻¹) and also falls in the range obtained by Meunier et al. (2017a) for another grass (Lolium multiflorum Lam., $\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{8}}$ to $\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{7}}$ m MPa⁻¹ s⁻¹). Our k_axial value cannot be compared to values of axial root conductance from the literature, as it transfers the water absorbed by roots in a single “big root” per group of five identical plants. The optimal value of k_sat was quite high (Table 1) but reportedly very correlated with λ (i.e., soil unsaturated hydraulic conductivity is proportional to k_sat but also to $S_e^{\mathit{\lambda }}$; van Genuchten, 1980), so that the low value of the latter compensated for the high value of the former; thus, they should be considered as effective rather than physical parameters.

3.2.3 Other sources of variability and observational error

Our treatment of the soil medium in this experiment (sieving and irrigation from the bottom) makes it laterally more homogeneous than natural soils. This method allowed us to specifically study the impact of the vertical gradients of δ_soil on δ_tiller. It also justified the use of a simplistic one-dimensional model adapted to the vertically resolved measurements. If lateral heterogeneity of soil water content remained and was accounted for, our predictions of root water uptake distribution, δ_tiller, and ψ_leaf would be altered. Observational errors in the gravimetric soil water content measurement (converted to soil water potential using the soil water retention curve) would also alter these predictions. In order to quantify the sensitivity of our simulated results to such heterogeneity or observational error, we varied the soil water content input by ±0.02 m³ m⁻³ at three critical depths (−0.9, −1.1, and −1.3 m, before interpolation) at the last observation time, during which measurements and simulations suggested that hydraulic lift occurred. Our results were mostly sensitive to soil water content alterations at −0.9 m, and they barely differed in response to alterations at −1.1 and −1.3 m, although the conclusions were not affected qualitatively. No statistically significant difference between the predicted and observed δ_tiller distributions for the overall dataset could be found when pooling three simulated δ_tiller randomly at each observation time (predicted and observed δ_tiller distributions were closest to differing when the soil water content was reduced by 0.02 m³ m⁻³ at a depth of 0.9 m; P > 0.01 in 76 cases out of 100 repeated drawings). Measured and simulated ψ_leaf remained very correlated in all cases (from an R² value of 0.69 to 0.74 when adding or removing 0.02 m³ m⁻³ at a depth of 0.9 m, respectively). Furthermore, when adding or removing 0.02 m³ m⁻³ at a depth of 0.9 m, cumulative water exudation at −0.9 m varied between 0.0019 and 0.0035 m³ m⁻³, uptake at −1.1 m varied between 0.0080 and 0.0108 m³ m⁻³, and the simulated change in the δ_soil ranged between 0.28 and 0.40 ‰, respectively.

Lateral heterogeneity in the soil water isotopic composition may as well occur at the microscopic scale. As water in micropores is less mobile than water in meso- and macropores (Alletto et al., 2006), it is likely that, in the lower half of the profile, the capillary rise of labeled water affected the composition of water in meso- and macropores more than in micropores. If roots have more access to meso- and macropore water, then the water absorbed by roots would be isotopically enriched compared with the “bulk soil water” characterized experimentally. The importance of this possible bias depends on soil texture and heterogeneity (e.g., the existence of more isolated “pockets” of soil or compact clusters) as well as on the speed of water mixing between mobile and immobile water fractions (Gazis and Feng, 2004). Including this process in the modeling would necessitate sufficient observations to estimate the aforementioned properties and ideally some quantification of the lateral heterogeneity of the soil water isotopic composition at the microscale.

The lateral heterogeneity of soil hydraulic properties and root distribution may also have participated in the generation of lateral soil water potential heterogeneities, particularly in undisturbed soils. If one had access to data on the lateral heterogeneity of soil properties and rooting density, it would be possible to simulate three-dimensional soil–root water flow with a tool such as R-SWMS (modeling “Root-Soil Water Movement and Solute transport”; Javaux et al., 2008), using a randomization technique for soil properties' distribution as in Kuhlmann et al. (2012), in order to obtain estimations of the relative importance of this type of heterogeneity on δ_tiller and ψ_leaf variability.

Unlike the tiller water isotopic composition, the leaf water potential turned out to be very sensitive to the transpiration rate in our simulations (see temporal fluctuations of gray lines in Fig. 5c) and not very sensitive to the root distribution (see small variations in leaf water potential across individuals in Fig. 5d). In this setup, this suggests that the hydraulic conductance of the soil–root system limited the shoot water supply more than the distribution of roots, as in Sulis et al. (2019). Simulated baseline (i.e., for uniform transpiration rates) leaf water potentials are shown as gray lines in Fig. 5c, and measured leaf water potentials are shown as a green line in the same panel. The fact that they match well, despite the high sensitivity of the leaf water potential to the transpiration rate, reinforces the idea that the transpiration rate was likely not spatially heterogeneous among the plant population. Therefore, the tiller water isotopic composition, whose sensitivity to the transpiration rate is already very low, was likely not affected by transpiration rate heterogeneity.

3.2.4 Do root water uptake profiles predicted by hydraulic and Bayesian models differ?

The root water uptake dynamics predicted by the mechanistic model are shown in Fig. 6a. The overall pattern of peaking water uptake in the lower part of the profile during daytime matched that of the statistical model, and the correlation coefficient of both model predictions was relatively high (R²=0.53) on average over the simulation period (see Fig. 7). The main differences were as follows: (i) in the upper soil layers where the soil water potential was lower than −1.5 MPa, the statistical model predicted water uptake, which is theoretically impossible given the leaf water potential above −0.4 MPa (van Den Honert, 1948); (ii) in the upper half of the profile, the physical model predicted exudation at a rate limited by the low hydraulic conductivity between the root surface and the bulk soil, with a peak at night, at a depth of −0.9 m (quantitative analysis in previous section); (iii) below a depth of −1.0 m, the water uptake rate predicted by the statistical model steadily increased with depth while that of the physical model was more uniform, likely due to axial hydraulic limitation (e.g., Bouda et al., 2018) counteracting the increasing soil water potential with depth. Note that the outcome of the statistical model may significantly depend on the definition of the a priori relative RWU (rRWU) profile. In the present study, we set it to follow a “flat” uniform distribution (i.e., rRWUj $=\mathrm{1}/\mathrm{10}$; see Appendix E); in other words, each layer was initially defined to contribute equally to RWU. Contrary to other studies (e.g., Mahindawansha et al., 2018), where the a priori rRWU profile was empirically constructed on the basis of soil water content and root length density profiles, we decided not to further arbitrarily constrain the Bayesian model for the sake of comparison with the physically based soil–root model.

Figure 6Time series of the profiles of root water uptake per unit soil volume (sink term, d⁻¹) computed with the physically based model. (a) Sum of sink terms across the 60 groups of the population. (b) Variability of the sink terms within the 60 groups of the population (1 standard deviation).

Download

Figure 7Time series of the profiles of root water uptake per unit soil volume (sink term, d⁻¹) computed with the statistical model SIAR (a). Panel (b) reports the variance of the estimated sink term (1 standard deviation).

Download