Plant hydraulic traits regulating productivity and drought resistance in boreal crops

Plant hydraulics play a critical role in coordinating photosynthesis, water use, and hydraulic safety when facing water scarcity. The hydraulic traits of boreal crops and their relevance to productivity and drought resistance, however, remain poorly researched. In this study, we combined multi-year carbon and water flux measurements from two Finnish fields representing contrasting boreal cropping systems – a perennial forage grass mixture and an annual cereal crop (oats) – with a novel modelling framework coupling plant photosynthesis, hydraulics, and soil moisture dynamics to decipher the field-level effective hydraulic traits (i.e., apparent hydraulic traits) of the crops and their association with crop performance under varying drought conditions. The model reproduced the observed carbon and water fluxes well and revealed relatively low maximum whole-plant hydraulic conductance and conservative water-use strategies in both crops. Under drought, low hydraulic efficiency prevented soil water depletion and thereby supported high productivity through drought avoidance, whereas its advantage diminished under wetter conditions. Both crops were substantially more susceptible to atmospheric dryness than to low soil moisture, highlighting vapor pressure deficit as an important constraint on crop productivity under boreal conditions. The remaining discrepancy between the model-inferred apparent hydraulic traits and available empirical measurements for temperate crops indicates the need for more direct hydraulic measurements of boreal crops. Our findings further underscore the importance of incorporating crop hydraulic processes in process-based crop models to gain physiological insights into crop performance and to enhance their accuracy and predictive power under current and future climates.

1 Introduction

Hydraulic traits are plant properties regulating water extraction and transport from the soil to the aerial organs, which play a fundamental role in plant productivity, notably when faced with water scarcity (Cardoso et al., 2025; Torres-Ruiz et al., 2024; Sperry and Love, 2015). Among these traits, two have been most widely studied. One is the maximum hydraulic conductance, which quantifies the efficiency with which water is transported through the whole plant or its individual organs (e.g. roots, stems and leaves) to meet photosynthetic demands (Holloway-Phillips and Brodribb, 2011a; Ocheltree et al., 2016; Gleason et al., 2021). The other is the water potential at which a plant loses 50 % of its maximum conductance. It is often used to quantify plant hydraulic safety, i.e. its resistance to drought-induced embolism, loss of water-transport capacity, and the resultant mortality (Choat et al., 2012; Anderegg et al., 2016). These two traits, which characterise, respectively, the hydraulic efficiency and safety of a plant, together play a critical role in shaping the plant's water use strategy.

A trade-off between hydraulic efficiency and safety has been widely documented (Ocheltree et al., 2016; Huo et al., 2022; Zhao et al., 2025), stating that plants capable of efficiently transporting water are often more vulnerable to embolism, and vice versa, because the wider and more interconnected xylem conduits that enable such efficient transport are also more prone to cavitation. An analogous trade-off has also been proposed for stomatal regulation (Henry et al., 2019), whereby maximising stomatal conductance promotes carbon gain but also increases transpirational water loss and the risk of desiccating plant water potential to a damaging level. However, the strength of such a trade-off is still hotly debated (Gleason et al., 2016; Sanchez-Martinez et al., 2020; Liu et al., 2021a), indicating that different species adopt different positions along the efficiency–safety spectrum to improve their fitness under variable environmental conditions. For instance, common crop species (mostly annual cereal crops from temperate and tropical regions) often exhibit high hydraulic efficiency but low hydraulic safety (Holloway-Phillips and Brodribb, 2011a; Lamarque et al., 2020; Corso et al., 2020; Canales et al., 2021; Gleason et al., 2021) to obtain high crop productivity, referred to as acquisitive hydraulic strategy. This is particularly the case for annual cereals (Liu et al., 2019). On the other hand, hydraulic traits indicating high safety but low efficiency (i.e., conservative hydraulic strategy) have been widely observed in species like boreal pine and spruce (Li et al., 2023) as well as perennial crops, which are known for their conservative hydraulic strategy to survive harsh winters and tolerate more stresses (Vico and Brunsell, 2018; Liu et al., 2019). Although plant hydraulic traits are often considered relatively conserved compared with more plastic plant traits such as photosynthetic capacity (Lamy et al., 2014; Ramírez-Valiente et al., 2025), substantial phenotypical plasticity has also been found in various plants under different conditions (Pritzkow et al., 2020; Sun et al., 2021; Alves et al., 2024; Matzner et al., 2014) and growing stages (Yao et al., 2021). For example, Sun et al. (2021) showed that barley grown under conventional no-tillage management exhibited higher hydraulic vulnerability than that under organic farming with intensive tillage.

Despite their significance, it remains elusive how hydraulic traits coordinate in boreal crops and modulate their productivity and drought response. Historically, water limitation has not been a primary concern for humid boreal agroecosystems. Therefore, hydraulic traits and drought tolerance have been less of a priority when selecting or breeding crops for the region, while other traits such as early maturity, cold tolerance, and disease resistance have received more attention (Roitsch et al., 2022). However, the situation has changed in recent years (Poque et al., 2025), as more frequent and severe drought events, characterised by both limited soil water availability (soil drought) and high atmospheric evaporative demand (atmospheric drought), have been observed and are projected to occur in the future (Spinoni et al., 2018; Aalto et al., 2026). The hydraulic traits governing the crop response to drought are therefore expected to be under strong selective pressure in the region. Thus, a deeper understanding of the impacts of different types of drought on boreal crops and their association with crop hydraulic traits is urgently needed.

Process-based models with an explicit representation of plant hydraulic processes have emerged only recently and have been shown to better describe plant responses to drought, especially atmospheric drought, than models without hydraulic processes (Kennedy et al., 2019; Liu et al., 2020; Kubota et al., 2026; Hu et al., 2026). They are also shown to be valuable tools for revealing plant's hydraulic traits through model inversion or calibration approaches (Wang et al., 2020; Gleason et al., 2022; Liu et al., 2021b; Lu et al., 2022). In particular, the use of field-level flux measurements to constrain the hydraulic parameters of such models has been proven effective in retrieving plant hydraulic traits at the field or ecosystem level (Liu et al., 2020). We refer to the average of the hydraulic traits over collective individuals of the same or different species at this scale as “apparent” hydraulic traits, in contrast to the organ- or individual-level hydraulic traits measured on single plants in laboratory or field studies.

In this study, we leveraged multi-year observations of carbon and water fluxes from two Finnish crop fields under a perennial forage grass mixture and an annual cereal (oats), respectively. Combining with a novel modelling framework that explicitly couples plant hydraulics, photosynthesis, and soil moisture dynamics, we infer apparent hydraulic traits of the two dominant and contrasting cropping systems in boreal agriculture, and examine their implications for crop functioning at field level. Specifically, our objective is to (i) evaluate the ability of the novel modelling framework to reproduce observed gross primary productivity and evapotranspiration, providing a basis for hydraulic trait inference; (ii) quantify and compare the apparent hydraulic traits of perennial and annual crops and assess their influence on productivity and water-use efficiency; and (iii) determine whether the two cropping systems differ in their resistance to drought and how these differences are associated with their hydraulic traits. In the following, we first describe the modelling framework, study sites and observations, and modelling experiments and workflow for addressing these objectives. We then evaluate the performance of the model, compare inferred hydraulic traits and their effects on crop functioning and drought responses, and finally discuss the broader implications and remaining uncertainties.

2 Materials and Methods

2.1 Structure of SPY-C

We developed a comprehensive soil-plant-atmosphere modelling framework, SPY-C, illustrated in Fig. 1. SPY-C is composed of four interlinked modules mostly from existing models: plant photosynthesis and hydraulics (Phydro), soil water balance (SpaFHy), soil organic carbon decomposition (Yasso), and a simple data-driven phenology and allocation module (PA). SPY-C here stands for SpaFHy-Phydro-Yasso for Carbon calculation.

Figure 1 The model structure of SPY-C. Amber and lavender boxes denote model components and input data, respectively. The light and dark green boxes denote carbon fluxes and storages. The light and dark blue boxes denote water fluxes and storages. The outline line style denotes the time step of each module and data stream: dotted = hourly, dashed = daily, and solid = time-invariant. The arrows between boxes denote data flows. The thicker arrows denote the two-way coupling between Phydro and SpaFHy.

2.1.1 Plant hydraulic and photosynthesis module (Phydro)

Phydro is based on Joshi et al. (2022), which unifies plant photosynthesis and hydraulics through a trait-based optimality theory implemented through a single cost-benefit function. The first principle in Phydro is the steady-state leaf water balance, that is, the hydraulic flow rate at which water enters the leaf (the left side of Eq. 1) equals the transpiration rate (the right side of Eq. 1):

$$\begin{array}{}\text{(1)}& -{\displaystyle \frac{K_{\mathrm{p}}}{\mathit{\eta }}}\underset{{\mathit{\psi }}_{\mathrm{soil}}}{\overset{{\mathit{\psi }}_{\mathrm{soil}}-\mathrm{\Delta }\mathit{\psi }}{\int }}P\left(\mathit{\psi }\right)d\mathit{\psi }=\mathrm{1.6}{g}_{\mathrm{s}}D.\end{array}$$

Here, K_p (unit: m) is the maximum whole-plant conductance per unit leaf area, which lumps the entire soil-to-leaf hydraulic pathway (including root, stem, and leaf) into a single effective conductance. η (unit: Pa s) is the dynamic viscosity of water, which is a function of air temperature (T) and pressure (P_atm). Note that T, P_atm, VPD, PAR, LAI, and c_a, which are defined in the following sections, are all input variables to the model. ψ_soil (Pa) is the soil water potential. Δψ (Pa) is the difference in water potential between the soil and the leaf (i.e., ψ_soil − ψ_leaf). g_s (mol CO₂ m⁻² s⁻¹) is the stomatal conductance to CO₂. D is the vapor pressure deficit (VPD) divided by P_atm. P(ψ) is the vulnerability curve (Eq. 2) that describes the decline in plant water conductance K(ψ) with its water potential (ψ) (Eq. 3):

$$\begin{array}{}\text{(2)}& {\displaystyle }P\left(\mathit{\psi }\right)=(\mathrm{1}/\mathrm{2}{)}^{(\mathit{\psi }/{\mathit{\psi }}_{\mathrm{50}}{)}^b},\text{(3)}& {\displaystyle }K\left(\mathit{\psi }\right)=K_{\mathrm{p}}P\left(\mathit{\psi }\right),\end{array}$$

where ψ₅₀ is the water potential at which 50 % of hydraulic conductivity is lost, and b determines the sensitivity of conductivity loss to water potential. Note that although the hydraulic framework underlying Phydro was developed primarily in the context of woody plants, its principles are also applicable to grass-like crops, where water is transported through functional xylem vessels within vascular bundles that provide a continuous hydraulic pathway from roots to leaves, despite the hollow stem structure and absence of secondary woody tissues (Heckwolf et al., 2015).

The second principle is the photosynthetic coordination hypothesis, which states that under typical day time conditions, the photosynthetic capacity is optimized so that the carboxylation-limited photosynthesis rate A_c (i.e., the left side of Eq. 4) and light-limited photosynthesis rate A_j (i.e., the right side of Eq. 4) are equal:

$$\begin{array}{}\text{(4)}& V_{\mathrm{cmax}}{\displaystyle \frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}^{*}}{c_{\mathrm{i}}+K_{\mathrm{M}}}}-R_{\mathrm{d}}={\displaystyle \frac{J}{\mathrm{4}}}\left({\displaystyle \frac{c_{\mathrm{i}}-{\mathrm{\Gamma }}^{*}}{c_{\mathrm{i}}+\mathrm{2}{\mathrm{\Gamma }}^{*}}}\right)-R_{\mathrm{d}}.\end{array}$$

Here, V_cmax (µmol CO₂ m⁻² s⁻¹) is the maximum carboxylation capacity of leaves. c_i is the leaf-internal CO₂ partial pressure divided by P_atm. K_M and Γ^* are the Michaelis–Menten coefficient for C₃ photosynthesis, and the photosynthetic light-compensation point, respectively. Both are a function of T and P_atm. R_d is the dark respiration, which is assumed to be proportional to V_cmax (Eq. 5). J (µmol CO₂ m⁻² s⁻¹) is the effective electron transport capacity of leaves, which is described by Eq. (6). We note that there are still uncertainties on the time scale on which the optimisation of photosynthetic capacity operates, spanning from instantaneous to weekly time scales (Joshi et al., 2022). Here, we assume that the optimization occurs at an hourly time scale to be tightly coupled with soil moisture dynamics during the day, which is also supported by a recent study by Drake et al. (2026).

$$\begin{array}{}\text{(5)}& R_{\mathrm{d}}=b_r{V}_{\mathrm{cmax}},\end{array}$$

where b_r is a unitless coefficient set to 0.015 in this study (Atkin et al., 2015).

$$\begin{array}{}\text{(6)}& J={\displaystyle \frac{\mathrm{4}{\mathit{\varphi }}_{\mathrm{0}}{I}_{\mathrm{abs}}}{\sqrt{{\left(\frac{\mathrm{4}{\mathit{\varphi }}_{\mathrm{0}}{I}_{\mathrm{abs}}}{J_{max}}\right)}^{\mathrm{2}}+\mathrm{1}}}},\end{array}$$

where ϕ₀ is the intrinsic quantum-yield efficiency, which is a function of T. J_max (µmol CO₂ m⁻² s⁻¹) is the maximum electron transport capacity of leaves under light saturation. I_abs ($\mathrm{\mu }\mathrm{mol}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{2}}$_leaf) is the incident photosynthetic photon-flux density per unit leaf area, which is calculated following Eq. (7):

$$\begin{array}{}\text{(7)}& I_{\mathrm{abs}}=\text{PAR}{f}_{\mathrm{apar}}/\text{LAI},\end{array}$$

where PAR is photosynthetically active radiation ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), LAI (m² m⁻²) is the total one-sided green leaf area per unit ground surface area. f_apar is the fraction of absorbed photosynthetically active radiation by the plant, which can be approximated following Beer's law:

$$\begin{array}{}\text{(8)}& f_{\mathrm{apar}}=\mathrm{1}-e^{-k\text{LAI}},\end{array}$$

where k is the “effective” canopy light attenuation parameter (including the effect of leaf reflectance and transmittance), which is assumed to be 0.5 in this study. For simplicity, the model does not resolve sunlit/shaded leaf fractions or multiple canopy layers.

The third principle is the steady-state leaf CO₂ balance, which states that the rate of photosynthesis (the left side of Eq. 9, based on Eq. 4), equals the rate at which CO₂ diffuses into the leaves (the right side of Eq. 9):

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{J}{\mathrm{4}}}\cdot {\displaystyle \frac{c_{\mathrm{i}}(\mathrm{1}-b_r)-({\mathrm{\Gamma }}^{*}+b_r{K}_{\mathrm{M}})}{c_{\mathrm{i}}+\mathrm{2}{\mathrm{\Gamma }}^{*}}}=g_{\mathrm{s}}{c}_{\mathrm{a}}(\mathrm{1}-\mathit{\chi }).\end{array}$$

Here, c_a is the atmospheric CO₂ partial pressure divided by P_atm and χ is the ratio of c_i and c_a.

Furthermore, a cost-benefit function F ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) is defined for plants as the benefit of net photosynthetic assimilation (A_j) minus the costs of maintaining the photosynthetic capacity (J_max) and the hydraulic pathway (Δψ), as given by Eq. (10):

$$\begin{array}{}\text{(10)}& F=A_j-\mathit{\alpha }{J}_{max}-\mathit{\gamma }\mathrm{\Delta }{\mathit{\psi }}^{\mathrm{2}}.\end{array}$$

Here, α is a unitless coefficient denoting the unit cost of maintaining photosynthetic capacity, such as the regeneration of RuBP and electron transport proteins. γ ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{2}}$) denotes the unit cost of maintaining the plant hydraulic pathway, such as the construction cost to build xylem vessels and fibres. Both are conceptual parameters that are inferred rather than directly measured. It is assumed that plants have adapted to their living environment in such a way that they can modulate their stomatal conductance g_s (and thus Δψ according to Eq. 1) and their electron-transport capacity J_max independently to maximise their net profit F.

Using Eqs. (1)–(9), the variable A_j in Eq. (10) can be expressed as a function of J_max and Δψ. As a result, finding the values of J_max and Δψ that maximize F becomes a two-dimensional bound-constrained optimisation problem and can be solved using a two-dimensional root finding algorithm (Morales and Nocedal, 2011).

With J_max and Δψ being determined, plant gross primary productivity (GPP, µmol CO₂ m⁻² s⁻¹) and transpiration (Tr, µmol H₂O m⁻² s⁻¹) can be calculated according to Eqs. (11) and (12):

$$\begin{array}{}\text{(11)}& {\displaystyle }\text{GPP}=\left(A_j+R_{\mathrm{d}}\right)\cdot \text{LAI},\text{(12)}& {\displaystyle }\text{Tr}=\mathrm{1.6}{g}_{\mathrm{s}}D\cdot \text{LAI},\end{array}$$

where factor 1.6 accounts for the relative diffusivity of water vapor to that of CO₂. V_cmax, R_d, and g_s can also be derived according to Eqs. (4), (5), and (1), respectively, once J_max and Δψ are determined.

2.1.2 Soil water balance module SpaFHy

The soil water balance in SPY-C describes the processes such as infiltration, drainage, and evapotranspiration (ET), the latter consisting of plant transpiration (Eq. 12) and soil evaporation. We adopted representations of rainfall interception, snow accumulation and melt, and a single-layer soil water balance from the SpaFHy (Spatial Forest Hydrology, Launiainen et al., 2019) model. The hydrological processes can be simulated at hourly resolution to enable close coupling with Phydro. The simulated volumetric soil water content is converted to soil water potential using the van Genuchten soil water retention function (Schaap and Van Genuchten, 2006) and passed to Phydro as an input variable (Fig. 1 and Eq. 13):

$$\begin{array}{}\text{(13)}& {\mathit{\psi }}_{\mathrm{soil}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_v}}{\left(S_{\mathrm{e}}^{-\frac{\mathrm{1}}{m}}-\mathrm{1}\right)}^{\frac{\mathrm{1}}{n}},\end{array}$$

where ψ_soil (kPa) is the soil water matric potential and α_v, n, and m are empirical shape parameters. The parameter m is commonly expressed as m = 1 − $\frac{\mathrm{1}}{n}$. S_e is the effective saturation as given by Eq. (14):

$$\begin{array}{}\text{(14)}& S_{\mathrm{e}}={\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}},\end{array}$$

where θ is the volumetric soil water content (m³ m⁻³), θ_s is the saturated soil water content, and θ_r is the residual water content.

To be consistent, the soil hydraulic conductivity (K_soil) affecting drainage is also derived using the van Genuchten model (Schaap and Van Genuchten, 2006), which can be written as:

$$\begin{array}{}\text{(15)}& K_{\mathrm{soil}}=K_{\mathrm{sat}}{S}_{\mathrm{e}}^{\mathrm{0.5}}{\left(\mathrm{1}-(\mathrm{1}-S_{\mathrm{e}}^{\mathrm{1}/m}{)}^m\right)}^{\mathrm{2}}.\end{array}$$

Here, K_sat (m s⁻¹) is the saturated hydraulic conductivity. When soil water potential is above field capacity, free drainage according to K_soil can occur at the lower boundary of the soil in SpaFHy. There is no explicit estimation of root water potential. Root water conductance has also been lumped into whole-plant hydraulic conductance (K_p) in Phydro. Root water uptake from soil is merely the transpiration demand given by Phydro met from the total soil water storage simulated in SpaFHy (Fig. 1). The calculation of canopy interception evaporation and soil evaporation follow the SpaFHy formulations, using Penman–Monteith equation and soil-water availability restricted Priestley–Taylor approach, respectively (Launiainen et al., 2019). Total ET is the sum of canopy and soil evaporation given by SpaFHy and transpiration flux given by Phydro (Fig. 1).

2.1.3 Phenology and allocation module

The phenology and carbon allocation module (PA) directly allocates daily assimilated carbon ($\stackrel{\mathrm{‾}}{\text{GPP}}$) into above- (i.e., leaf) and below-ground (i.e., root) carbon pools according to specific allometric ratios, following the concept from Sierra et al. (2022) (Eqs. 16 and 17).

$$\begin{array}{}\text{(16)}& {\displaystyle }{\text{CF}}_{\mathrm{gpp}:\mathrm{leaf}}=\stackrel{\mathrm{‾}}{\text{GPP}}\cdot a_{\mathrm{leaf}},\text{(17)}& {\displaystyle }{\text{CF}}_{\mathrm{gpp}:\mathrm{root}}=\stackrel{\mathrm{‾}}{\text{GPP}}\cdot a_{\mathrm{root}},\end{array}$$

where $\stackrel{\mathrm{‾}}{\text{GPP}}$ denotes the daily accumulation of hourly GPP flux derived by Eq. (11), CF_{gpp : leaf} and CF_{gpp : root} (µmol m⁻² d⁻¹) are the carbon fluxes from $\stackrel{\mathrm{‾}}{\text{GPP}}$ to leaf and root, and a_leaf and a_root are the respective proportions of $\stackrel{\mathrm{‾}}{\text{GPP}}$ allocated to leaf and root carbon storage at each daily time step. For cereal crops, grains can constitute a significant proportion of the aboveground biomass. The grain filling flux (CF_{gpp : grain}) is thus considered and described following Eq. (18):

$$\begin{array}{}\text{(18)}& {\text{CF}}_{\mathrm{gpp}:\mathrm{grain}}=\stackrel{\mathrm{‾}}{\text{GPP}}\cdot a_{\mathrm{grain}},\end{array}$$

where a_grain denotes the proportion of $\stackrel{\mathrm{‾}}{\text{GPP}}$ allocated to grain carbon pool at each daily time step. The details on how a_leaf, a_root, and a_grain are derived and constrained by observation can be found in Methods S1 in the Supplement.

The newly allocated carbon is then subject to autotrophic respiration (AR), consisting of growth (CF_gr,leaf, CF_gr,root, CF_gr,grain) and maintenance respiration (CF_mr,leaf, CF_mr,root, CF_mr,grain) of the plant organs. In addition, it is assumed that both leaf and root lose their carbon through turnover (CF_{leaf : litter}, CF_{root : litter}), which produces litter carbon input to soil (Fig. 1). More detailed formulation of these fluxes can also be found in Methods S1.

The changes in carbon storage of leaf (ΔCS_leaf), root (ΔCS_root), and grain (ΔCS_grain) at each time step are given by Eqs. (19), (20), (21), respectively:

$$\begin{array}{}\text{(19)}& {\displaystyle }\mathrm{\Delta }{\text{CS}}_{\mathrm{leaf}}={\text{CF}}_{\mathrm{gpp}:\mathrm{leaf}}-{\text{CF}}_{\mathrm{leaf}:\mathrm{litter}}-{\text{CF}}_{\mathrm{mr},\mathrm{leaf}}-{\text{CF}}_{\mathrm{gr},\mathrm{leaf}},\text{(20)}& {\displaystyle }\mathrm{\Delta }{\text{CS}}_{\mathrm{root}}={\text{CF}}_{\mathrm{gpp}:\mathrm{root}}-{\text{CF}}_{\mathrm{root}:\mathrm{litter}}-{\text{CF}}_{\mathrm{mr},\mathrm{root}}-{\text{CF}}_{\mathrm{gr},\mathrm{root}},\text{(21)}& {\displaystyle }\mathrm{\Delta }{\text{CS}}_{\mathrm{grain}}={\text{CF}}_{\mathrm{gpp}:\mathrm{grain}}-{\text{CF}}_{\mathrm{mr},\mathrm{grain}}-{\text{CF}}_{\mathrm{gr},\mathrm{grain}}.\end{array}$$

Note that root carbon is tracked as a prognostic pool (Eq. 20) subject to allocation, respiration, and turnover, and it supplies litter to soil decomposition model Yasso. However, this root carbon does NOT dynamically change rooting depth or hydraulic traits: the root-zone depth (0.6 m) and K_p are fixed within a simulation, hence is decoupled from the plant hydraulic routine in the model.

In the current model version, phenology is prescribed rather than dynamically simulated: the seasonal canopy development is imposed through the daily LAI interpolated from Sentinel-2 product (see more details in Sect. 2.2.2), and reported management events such as harvest, grazing, fertilisation (see more details in Methods S2).

2.1.4 Soil decomposition module Yasso

We describe soil carbon balance dynamics and heterotrophic respiration (HR) using a daily timestep version of Yasso20 soil carbon model, which formulation has been described in previous studies (Viskari et al., 2022; Heimsch et al., 2024). In our modelling framework, the litter input (CF_{leaf : litter}, CF_{root : litter}) to Yasso is provided by the PA module (Eqs. S8 and S9 in the Supplement) and the litter quality (i.e., carbon fractions) is prescribed depending on crop types following the study by Palosuo et al. (2015).

Combining the HR derived from Yasso with the autotrophic respiration flux (AR) calculated by the PA module (Eq. S7) sum up to the total ecosystem respiration (RECO):

$$\begin{array}{}\text{(22)}& \text{RECO}=\text{AR}+\text{HR}.\end{array}$$

Further, combining $\stackrel{\mathrm{‾}}{\text{GPP}}$, the net ecosystem exchange (NEE) at a daily time step (Fig. 1) can be expressed as:

$$\begin{array}{}\text{(23)}& \text{NEE}=\text{RECO}-\stackrel{\mathrm{‾}}{\text{GPP}}.\end{array}$$

2.1.5 Coupling between model components

As shown in Fig. 1, SPY-C runs on two nested time steps. Phydro and SpaFHy are integrated hourly and two-way coupled: each hour, SpaFHy's soil moisture is converted to soil water potential before passing to Phydro; Phydro solves stomatal conductance and thus transpiration (Eq. 12); the transpiration is then sent to SpaFHy and is withdrawn from the soil water store together with soil evaporation, updating soil moisture for the next hour. The PA and Yasso modules run daily: hourly GPP from Phydro is accumulated to a daily total, allocated by PA to leaf, root, and grain pools (Eqs. 16–18), and the resulting litter drives Yasso, yielding RECO and NEE on daily scale (Eqs. 22–23). The hourly to daily link is GPP accumulation, while the daily to hourly link is the prescribed canopy state (i.e., LAI). Daily leaf area growth enters the hourly Phydro calculations through absorbed PAR (Eqs. 7–8) and is used to scale GPP and transpiration to the canopy level (Eqs. 11–12), modulating both carbon and water fluxes every hour. In the current version of SPY-C, daily LAI is an observation-driven boundary condition (Fig. 1), rather than a prognostic variable predicted in the PA module. Consequently, two-way coupling between PA and Phydro/SpaFHy is not yet fully implemented.

2.2 Model experiments

2.2.1 Study sites and observation data

Two agricultural fields on mineral (i.e., non-peat) soils in Southern Finland were used in this study (Fig. 2a). The sites were selected to represent two agronomically important and contrasting boreal cropping systems: a perennial forage grass mixture and an annual cereal crop. They differ in crop life form, soil texture, and management, and both are equipped with multi-year eddy covariance flux measurements that support model evaluation and inversion. The Qvidja site (60.29550° N, 22.39281° E; elevation 5 m) is on the southwestern coast of Finland, which has continuous eddy covariance (EC) flux measurements since May 2018 (Heimsch et al., 2021). During the study period, the footprint area of flux measurement was covered by a mixture of perennial grasses and legumes, dominated by timothy (Phleum pratense) and meadow fescue (Festuca pratensis). The soil is a relatively fine textured clay loam. The Hauho site (61.13977° N, 24.58401° E; elevation 89.5 m) has continuous EC measurements since June 2022 (Vira et al., 2025b). The footprint area of the flux measurement was primarily covered by oats (Avena sativa); in 2023, Italian ryegrass (Lolium multiflorum) was additionally sown as a cover crop during the growing season. Unlike Qvidja, Hauho has a coarser-textured silty loam soil. Together, these contrasting sites provide a useful test of the model across different crop types and growing conditions and allow us to examine whether field-scale hydraulic traits and water-use strategies differ between perennial and annual boreal cropping systems.

Figure 2 (a) Locations of the study sites; (b) five representative soil water retention curves for Finnish forest soils (dashed lines) (Launiainen et al., 2022), together with the observed ranges of volumetric soil moisture at the study sites (colour bars); (c, d) Sentinel-2 leaf area index (LAI) observations (dots) and daily LAI interpolated using Gaussian process regression (line) (Vira et al., 2025b) for Qvidja and Hauho, respectively, and (e, f) daily vapor pressure deficit and soil moisture observed at Qvidja and Hauho, respectively. Management events included in the simulations are indicated in panels (c) and (d): H, Harvest; O, Organic fertilizer; G, Grazing; M, Mowing. Drought events identified using the Standardised Precipitation-Evapotranspiration Index (SPEI < −1.5) from the global SPEI database (https://spei.csic.es/database.html, last access: 15 May 2025) are denoted by brown bars in panels (e) and (f).

The EC instrumentation at both sites was the same, measuring ${\mathrm{CO}}_{\mathrm{2}}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ mixing ratio (LI-7200, LI-COR Biosciences, USA), and three-dimensional wind at 10 Hz frequency (u-Sonic3 Scientific, METEK GmbH, Germany). Besides flux measurements, auxiliary meteorological measurements close to the EC tower were also performed, such as air temperature and relative humidity (Humicap HMP155, Vaisala Oyj, Finland), photosynthetically active radiation, soil temperature (5 or 10 cm below the surface; Pt100 IKES sensors, Nokeval Oy, Finland), and soil moisture (5 or 10 cm below the surface; ML3 ThetaProbe sensor, Delta-T Devices Ltd., UK). The half-hourly turbulent ${\mathrm{CO}}_{\mathrm{2}}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ fluxes were calculated using the Eddypro software (v.7.0.9, LI-COR Biosciences, USA), following standard quality control and filtering protocols (Vira et al., 2025b). The gap-filling of the flux data was performed using the method developed by Vekuri et al. (2023). As environmental drivers for gap-filling, we used air and soil temperature, photosynthetically active radiation, soil moisture, and VPD at both sites, with LAI and precipitation additionally included for the Qvidja site. NEE was partitioned into GPP and RECO. RECO was assumed to follow an Arrhenius-type temperature response function (Lloyd and Taylor, 1994) and parameters were fitted following Reichstein et al. (2005). The parametrisation of the GPP was based on a non-linear radiation response function as in Gerin et al. (2023). For more detailed information about the sites and data, we refer to the aforementioned references, as well as the Field Observatory (https://fieldobservatory.org, last access: 9 November 2025), which is an open platform for real-time monitoring of crop fields in Finland (Nevalainen et al., 2022a).

2.2.2 Model initialization

All input data required to initialize SPY-C are summarised in Methods S2. Specifically, the meteorological forcing data for the two sites (e.g., T, P_atm, PAR, and VPD) were obtained directly or derived from ERA5-Land reanalysis data (Muñoz-Sabater et al., 2021) using their geographical locations. ERA5-Land was used because it provides a gap-free and internally consistent set of all meteorological drivers required by SPY-C over the entire simulation period. Its use also facilitates the transferability of the modelling framework to sites and periods where complete in situ meteorological observations are unavailable, supporting its intended future applications. ERA5-Land showed good agreement with the on-site measurements for key meteorological variables at our study sites, such as air temperature, precipitation, and VPD (data not shown).

The annual LAI cycle was prescribed by interpolating the LAI extracted from Sentinel-2 remote sensing products to daily values (Fig. 2c) using a Gaussian process regression as described in Vira et al. (2025b). The retrieval of LAI was based on Sentinel-2 level-2A (bottom-of-atmosphere reflectance) products and follows the LAI algorithm (v1.1) developed by Weiss et al. (2020) for the European Space Agency's SNAP software (Science Toolbox Exploitation Platform) and implemented in the satellitetools Python package (https://doi.org/10.5281/zenodo.5993291, Nevalainen and Vira, 2022). The retrieval provides a total single-sided green LAI at 20 × 20 m, which was then aggregated over the EC flux footprint at each site to obtain the field-level total LAI that drives the model. This aggregation is consistent with the model's big-leaf canopy representation and with the spatial scale at which modelled fluxes were compared with EC observations. Gaussian process regression was used to interpolate satellite-derived LAI to a daily time step because it accommodates irregularly spaced observations while providing both smooth interpolation and an explicit estimate of prediction uncertainty. In the presence of data gaps, the interpolation gradually approaches the mean of the Gaussian process, which was set to zero. This approach was found to perform well under Finnish conditions (Vira et al., 2025b), since extended data gaps were almost exclusively during winter, when the fields were snow-covered and LAI was effectively zero. These gaps are reflected in the larger prediction uncertainty of the interpolated time series.

All model parameters used in SPY-C and their values for the two study fields are listed in Table S2 in the Supplement. Most parameters were kept at their default values and were identical between sites, except for the SpaFHy soil hydraulic parameters, which were adjusted to represent site-specific soil texture, and the Phydro plant hydraulic parameters, which were calibrated for each crop type (Sect. 2.2.3). The soil depth for root water uptake is set to 0.6 m, which was found to represent available root-zone water storage (Methods S3) and reproduce observed surface soil moisture dynamics at both sites (Fig. S1 in the Supplement). Parameters for the soil water retention function (Eq. 13) were selected by comparing observed volumetric soil moisture ranges with soil water retention curves (SWRCs) established for major Finnish soil types by Launiainen et al. (2022) (Fig. 2b). SWRC5 and SWRC4 best encompassed the observed soil moisture ranges at Qvidja and Hauho, respectively, and were therefore used in the baseline simulations. The lower wilting point of SWRC4 is consistent with the coarser soil texture at Hauho. SWRC3 was additionally tested at Qvidja because, at a given volumetric soil moisture, it produces lower soil water potentials and thus stronger water stress than SWRC5 (Fig. 2b), allowing crop responses to stronger soil water stress to be evaluated.

The model simulations started on 1 January 2018 for Qvidja and 1 January 2022 for Hauho, and both ended on 31 December 2023. As the simulations were initialized in winter, soil moisture was set to full saturation, consistent with observations at both sites (Fig. 2e, f). Soil organic carbon pools in SPY-C were initialized in two steps. First, the relative fractions of the four faster-cycling Yasso carbon pools were assumed to be at steady state under a fixed litter carbon composition according to Palosuo et al. (2015) and mean climatological conditions. The slowest-cycling pool (humus) was assigned a prescribed fraction of legacy carbon that does not equilibrate with contemporary litter inputs. Second, the relative carbon pool sizes were rescaled to match the observed total soil organic carbon stock in the 0–100 cm soil layer: 16.0 kg m⁻² at Qvidja (measured range: 14.3–18.8 kg m⁻²) (Heimsch et al., 2021), and 12.0 kg m⁻² at Hauho (measured range: 8.5–15 kg m⁻²) (Fig. S2). We tuned the legacy carbon fraction to match the observed RECO flux at both Qvidja and Hauho. The resulted values were 0.3 for Qvidja and 0.1 for Hauho. The lower legacy carbon fraction at Hauho is in line with its coarser soil texture and thus lower capacity to stabilize carbon through mineral protection (Haddix et al., 2020).

2.2.3 Parameter perturbation ensemble and sensitivity analysis

Although the studied fields (i.e., EC flux footprints) contained mixtures of crop species, these species were spatially well mixed and were therefore represented in the model as a single effectively homogeneous canopy. This treatment is consistent with the model's field-scale representation and allows the inference of apparent hydraulic traits, whereby one effective trait set is estimated for each footprint and study period. To infer these traits and assess their impacts on productivity and drought resistance at each field, we selected four tunable Phydro parameters for parameter perturbation ensemble (PPE) experiments. Three parameters, K_p, ψ₅₀, and γ, characterize plant hydraulic properties or strategies, whereas α represents the photosynthetic cost associated with maintaining photosynthetic capacity.

Figure 3 Workflow for using the SPY-C parameter perturbation ensemble (PPE) experiments at two boreal crop sites to address the objectives of this study. QV = Qvidja; HA = Hauho; noVPD = No atmospheric drought; fullSM = Fully saturated soil.

For each PPE, four evenly spaced values were drawn from the prior range of each parameter and combined factorially, yielding 256 ensemble members. The prior ranges of each parameter were defined based on the literature. For the conceptual Phydro parameters α and γ, the ranges were directly from Joshi et al. (2022), who examined a broad set of plant species: [0.03,0.12] for α and [0.1,5] ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{Pa}}^{-\mathrm{2}}$) for γ. For the hydraulic parameters K_p and ψ₅₀, Joshi et al. (2022) reported ranges of $[\mathrm{0.2},\mathrm{15.9}]\times {\mathrm{10}}^{-\mathrm{16}}$ m and $[-\mathrm{2.2},-\mathrm{0.4}]$ MPa, respectively. Because relatively few herbaceous species were included in that study, we refined these priors using observational studies of herbaceous crops (Table S3). Accordingly, K_p was set to $[\mathrm{0.5},\mathrm{10}]\times {\mathrm{10}}^{-\mathrm{16}}$ m and ψ₅₀ to $[-\mathrm{3},-\mathrm{1}]$ MPa. These ranges encompass values reported for crops, while remaining sufficiently broad to avoid overly constraining the parameter space.

The main PPE experiments performed in this study and their purposes are summarised in Table 1 and Fig. 3, respectively. For Qvidja, a baseline PPE experiment using SWRC5 (QV_ctrl) and a sensitivity PPE experiment using SWRC3 (QV_swrc3) were conducted. Since both sites experienced drought events of varying severity during the study periods (Fig. 2e, f), two additional PPE experiments were conducted to better assess the impacts of drought on crop performance: one with soil moisture maintained at full saturation throughout the simulation period (QV_fullSM) and one with VPD constrained to very low values (≤ 50 Pa; QV_noVPD). The effects of soil (atmospheric) drought can then be deduced by the differences between QV_fullSM (QV_noVPD) and QV_ctrl. The PPE experiments for Hauho (HA_ctrl, HA_fullSM, HA_noVPD) followed the same design, except that SWRC4 was used in all PPEs (Table 1).

Table 1 Summary of the settings for the major parameter perturbation ensemble (PPE) experiments performed and analyzed in this study.

The PPE experiments QV_ctrl and HA_ctrl were used for model inversion to infer the apparent hydraulic traits of the crops at each study site by identifying the ensemble members that best reproduced the observed daily GPP and ET fluxes (Fig. 3). Because plant photosynthesis and hydraulics are tightly coupled, we employed both GPP and ET fluxes to provide stronger constraints on the four perturbed Phydro parameters than either flux alone. Model performance for each PPE member was evaluated using the Pearson correlation coefficient (cor) and normalized mean square error (NMSE) (Eqs. 24 and 25). The five members that ranked highest in cor and lowest in NMSE jointly for GPP and ET were selected, and their corresponding Phydro parameter sets were interpreted as the best-fit apparent hydraulic traits of each crop. Details of the ranking procedure are provided in Methods S4. Both cor and NMSE were calculated using daily data to emphasize the model's ability to reproduce daily variability and subseasonal drought responses. Using monthly data yielded similar parameter rankings and is therefore not shown here.

$$\begin{array}{}\text{(24)}& {\displaystyle }\text{cor}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left(X_{\mathrm{obs},i}-{\mathit{\mu }}_{\mathrm{obs}}\right)\left(X_{\mathrm{sim},i}-{\mathit{\mu }}_{\mathrm{sim}}\right)}{\sqrt{{\sum }_{i=\mathrm{1}}^n{\left(X_{\mathrm{obs},i}-{\mathit{\mu }}_{\mathrm{obs}}\right)}^{\mathrm{2}}{\sum }_{i=\mathrm{1}}^n{\left(X_{\mathrm{sim},i}-{\mathit{\mu }}_{\mathrm{sim}}\right)}^{\mathrm{2}}}}}\text{(25)}& {\displaystyle }\text{NMSE}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n{\left(X_{\mathrm{obs},i}-X_{\mathrm{sim},i}\right)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n{X}_{\mathrm{obs},i}^{\mathrm{2}}}},\end{array}$$

where i is the time step and n the total number of time steps; X_obs,i and X_sim,i are the observed and simulated values at each time step, respectively; and μ_obs and μ_sim are their respective mean values over all time steps. Here, the periods used to calculate these metrics correspond to the model inversion periods for inferring apparent hydraulic traits and were therefore selected separately for each crop. Specifically, the inversion periods were January 2019–December 2022 for the perennial grass mixture at Qvidja and June 2022–December 2023 for oats at Hauho. For the cover crop at Hauho, the inversion period was September–October 2023, when it was actively growing and had become dominant in the field. At Qvidja, year 2018 and 2023 were excluded from the inversion period because sowing events in those years may have altered species composition. Instead, these years were used as independent tests of whether the inferred apparent hydraulic traits were transferable beyond the inversion period. Additional analyses using alternative inversion periods are described in Methods S4.

To detect the effects of crop hydraulic traits on crop performance and drought response, the Morris sensitivity analysis (Morris, 1991) was also conducted using the PPE experiments for both sites, including both the control and counterfactual PPEs (Fig. 3). The sensitivities of monthly GPP, ET, and water use efficiency (WUE) to the 4 Phydro parameters perturbed were quantified. Here, WUE was defined at both the plant level (pWUE, unit: µmol CO₂ (mmol H₂O)⁻¹) as given by:

$$\begin{array}{}\text{(26)}& \text{pWUE}=\text{GPP}/\text{Tr},\end{array}$$

and the leaf level, often called intrinsic WUE (iWUE, unit: µmol CO₂ (mol H₂O)⁻¹) as given by:

$$\begin{array}{}\text{(27)}& \text{iWUE}=\left(A_j+R_{\mathrm{d}}\right)/\mathrm{1.6}{g}_{\mathrm{s}}.\end{array}$$

The Morris method is an efficient global sensitivity analysis technique for parameter screening (Campolongo et al., 2007), especially when the sampling size is limited (Gan et al., 2014). It is conceptually similar to the commonly used one-at-a-time sensitivity analysis, yet also considers the interactions among parameters. The overall effect and interaction effect of each parameter are approximated by the mean (μ) and standard deviation (σ) of the elementary effects of each parameter sampled. In addition, the mean of the absolute effects (μ^*) is also used to provide a more robust measure of overall parameter importance when elementary effects have opposite signs (Campolongo et al., 2007). The analysis was performed using the Python package SALib version 1.4.8 (Herman and Usher, 2017).

3 Results

3.1 Apparent hydraulic traits at Qvidja and Hauho

We identified five Phydro parameter sets that best reproduced daily GPP and ET at both sites. The results are shown in Table 2. For the perennial grass mixture at Qvidja and oats at Hauho, these best-fit parameter sets, interpreted as apparent field-scale hydraulic traits, were consistently characterized by low K_p and high γ, indicating relatively low hydraulic efficiency and a conservative water-use strategy. A low α (i.e., the cost of maintain photosynthetic capacity) was also required to reproduce the observed fluxes for both crops, implying relatively high photosynthetic capacity (Fig. S3). Compared to the perennial grass mixture, oats showed slightly lower best-fit α and γ values (Table 2), pointing to higher photosynthetic capacity and a slightly more acquisitive hydraulic strategy. The cover crop at Hauho exhibited higher inferred K_p, indicating greater hydraulic efficiency than oats or the perennial grass mixture. In contrast, no clear preference for ψ₅₀, a commonly used indicator of hydraulic safety, emerged across the studied crops, whereas the other three Phydro parameters were more strongly constrained.

Table 2 Summary of the prior ranges of the Phydro parameters used in the parameter perturbation ensemble (PPE) experiments, and the five best-fit PPE members for the forage grass mixture at Qvidja (based on QV_ctrl) and for oats (cash crop) and Italian ryegrass (cover crop) at Hauho (based on HA_ctrl). These parameter sets reflect effective field-scale hydraulic traits (i.e., apparent hydraulic traits) of the dominant crop community within each EC footprint and model-inversion period. K_p: maximum whole-plant conductance per unit leaf area (10⁻¹⁶ m). ψ₅₀: water potential at which 50 % of plant hydraulic conductivity (MPa) is lost. α: unit cost of maintaining photosynthetic capacity. γ: the unit cost of maintaining plant hydraulic pathways. cor_gpp and cor_et: correlation coefficients between model and observation for GPP and ET, respectively (Eq. 24). nmse_gpp and nmse_et: normalized mean square error between modelled and observed GPP and ET, respectively (Eq. 25).

Figure 4a–h shows how the five best-fit PPE members selected for the perennial grass mixture and oats performed in simulating daily GPP and ET at the study sites. Overall, the simulated GPP and ET matched the observed fluxes very well, with cor ≈ 0.94 and NMSE ≤ 0.14 for the Qvidja site, and cor ≥ 0.91, NMSE ≤ 0.11 for the Hauho site. Only a slight overestimation of daily ET was observed at Qvidja (Fig. 4c, d), and a mismatch in ET occurred at Hauho in early spring of 2023 when the model failed to capture the observed high ET (Fig. 4g), largely due to a bias in the timing of snowmelt onset in the meteorological forcing data (Fig. S7). The observed diurnal cycles of GPP and ET in summer were also well captured by the best-fit PPE members (Fig. 5), except for a slight underestimation (overestimation) of GPP (ET) in the early morning and late afternoon, probably resulting from a different light absorption rate from what was assumed in the model (see Eq. 8) under low sun angles.

Figure 4 (a–d) Daily mean gross primary production (GPP) and evapotranspiration (ET) modelled by the five best-fit PPE members from QV_ctrl and HA_ctrl (see Table 2) and their comparison with observation at Qvidja (a–d) and Hauho (e–h), respectively. The modelling results based on the five best-fit PPE members of QV_swrc3 can be found in Fig. S4.

Figure 5 Diurnal cycle of gross primary productivity (GPP) and evapotranspiration (ET) fluxes in summer (June, July and August) modelled by the average of the five best-fit PPE members from QV_ctrl and HA_ctrl and their comparison with observations at Qvidja (a, b) and Hauho (c, d), respectively. The violin denotes the distribution of hourly values across all summer days, and the horizontal lines within the violins indicate the 25th, 50th (median), and 75th percentiles.

We further examined the results of the five best-fit PPE members in the context of the full PPE experiments spanning a wide hydraulic and photosynthetic trait space (Fig. 6). For GPP, the best-fit PPE members generally fell in the upper range of the full PPE (Fig. 6a, b), except in May at Qvidja, when they shifted towards the lower end of the spread. But for ET, the best-fit PPE members tended to lie in the lower range (Fig. 6c, d), except during dry summer months, such as June 2018, 2019, and 2023 and July 2021 at Qvidja and July–August 2022 at Hauho. Consequently, their monthly mean pWUE and iWUE were close to the upper limit of the PPE, except in those same drought periods (Fig. 6e–h). Notably, the soil moisture simulated by the best-fit PPE members consistently lies in the upper range of the PPE (Fig. 6i–j), far above the wilting point.

Figure 6 Violin plots showing the spread of monthly mean gross primary production (GPP), evapotranspiration (ET), water use efficiency (pWUE and iWUE), and bulk root-zone soil moisture in the PPE for Qvidja (QV_ctrl) (a, c, e, g, i) and Hauho (HA_ctrl) (b, d, f, h, j). The width of each violin denotes the probability density of the ensemble members yielding a given values. The red circles represent the mean of the five best-fit PPE members. The red bars in panels (a)–(d) denote the observations. The brown boxes in panels (i) and (j) denote drought events defined by a Standardised Precipitation-Evapotranspiration Index (SPEI) below −1.5. The dotted lines in panels (i) and (j) denote the wilting point (i.e., the soil moisture at a matric potential of −1.5 MPa). Results from the additional Qvidja PPE experiment, QV_swrc3, can be found in Fig. S5.

3.2 Sensitivity of crop performance to hydraulic traits

The Morris sensitivity analysis of GPP, ET, pWUE, and iWUE to the four Phydro parameters in QV_ctrl and HA_ctrl is shown in Fig. 7a–h. The sensitivities are context-dependent and vary seasonally, but show broadly similar patterns at the two study sites. GPP is dominated by a negative sensitivity to α (the cost of maintaining photosynthetic capacity) throughout the growing season, except in June–July, when K_p and γ exert strong negative and positive effects, respectively, suggesting that lower hydraulic efficiency may promote crop productivity during these months. ET, pWUE, and iWUE exhibit more pronounced sensitivities to the hydraulic parameters. In particular, K_p (γ) negatively (positively) influences pWUE and iWUE in May, September, and October, while the effects often reverse in summer, suggesting a positive effect of hydraulic efficiency on pWUE and iWUE during peak growth, opposite to its effect on GPP.

Figure 7 Morris sensitivity analysis of monthly mean gross primary production (GPP) (a, b), evapotranspiration (ET) (c, d), plant-level water use efficiency (pWUE) (e, f), intrinsic water use efficiency (iWUE) (g, h) and bulk root-zone soil moisture (i, j) to Phydro parameters (K_p, ψ₅₀, α, and γ) in the PPE experiments for Qvidja (QV_ctrl) and Hauho (HA_ctrl) during the growing season (May–October). The mean of the absolute elementary effects (μ^*) is denoted by the size of the circles. The mean of the elementary effects (μ) is denoted by the colour shading of the circles. The standard deviation of the elementary effects (σ) is denoted by the colour shading of the rectangles. Results from the additional Qvidja PPE experiment, QV_swrc3, are shown in Fig. S6.

We note that while the effects of α are relatively independent of the other parameters (low σ), the effects of K_p and γ are more dependent on parameter interactions (high σ). We observe little influence of ψ₅₀ on GPP, ET, and WUE at either study sites, but a noticeable positive effect on soil moisture in summer (Fig. 7i, j). A consistently negative (positive) effect of K_p (γ) on soil moisture (Fig. 7i, j) suggests that lower crop hydraulic efficiency may reduce crop water extraction and thereby conserve soil water.

3.3 Response to soil and atmospheric drought

As illustrated in Fig. 8, the removal of soil drought and VPD increases GPP during the summer months in most PPE members (Fig. 8a–d) with the magnitude of the increase positively related to K_p. In contrast, transpiration responds differently to the absence of soil and atmospheric droughts (Fig. 8e–h): removing soil drought produces a consistently positive change, whereas reducing atmospheric drought induces a large spread with both positive and negative responses. Changes in transpiration are again most sensitive to K_p with positive effects mainly in summer.

Figure 8 Modelled response of monthly mean gross primary production (GPP) and transpiration to the removal of soil drought (FullSM-CTRL) (a, b, e, f) and atmospheric drought (NoVPD-CTRL) (c, d, g, h) at Qvidja (a, c, e, g) and Hauho (b, d, f, h). FullSM: QV_fullSM and HA_fullSM; NoVPD: QV_noVPD and HA_noVPD; CTRL: QV_ctrl and HA_ctrl. The upper panel of each subplot shows the spread of the PPE members (violins), and the mean of the five best-fit PPE members (red dots). The lower panel of each subplot shows the Morris sensitivity of the simulated responses to the Phydro parameters (K_p, ψ₅₀, α, and γ).

The five best-fit PPE members display little change in GPP and transpiration when soil drought is removed (Fig. 8a, b, e, f), indicating that crop photosynthesis and transpiration may have not suffered from soil drought, even though drought events were widely observed at the study sites (Fig. 2e, f). In contrast, a discernible rise in GPP and a reduction in transpiration during summer is shown in these best-fit PPE members when atmospheric drought is removed (Fig. 8c, d, g, h), suggesting that VPD may have exerted a more prevalent effect on crop performance than soil drought at the study sites.

The effect of VPD is further manifested in Fig. 9, where the five best-fit PPE members from QV_ctrl, HA_ctrl and QV_swrc3 display a clear impact of VPD on photosynthetic assimilation rate, pWUE, and iWUE (Fig. 9b, n, p) at the daily time scale. This can be associated with its strong regulation of stomatal conductance (Fig. 9d) and particularly the soil-to-leaf water potential difference (Fig. 9l), which further impact leaf internal-to-external CO₂ ratio (χ) (Fig. 9h) and thus carboxylation and electron-transport capacity (Fig. 9f, j). The modelled effect of VPD appears to be strongly nonlinear: while the impact on hydraulic behaviours, such as stomatal conductance, pWUE, and iWUE, mainly occurs during the initial increase in VPD (VPD ≤ 0.8 kPa) (Fig. 9d, h, i, n, p), the downregulation of photosynthetic activity emerges primarily at higher VPD (VPD ≥ 0.8 kPa) (Fig. 9b, f, j).

Figure 9 Daily response of photosynthetic and hydraulic variables to soil water potential (first and third columns) and vapor pressure deficit (second and fourth columns) simulated by the five best-fit PPE members selected from QV_ctrl, QV_swrc3, and HA_ctrl. (a, b) Assimilation rate (A, µmol CO₂ m⁻² s⁻¹); (c, d) stomatal conductance (g_s, mol CO₂ m⁻² s⁻¹); (e, f) carboxylation capacity (V_cmax, µmol CO₂ m⁻² s⁻¹); (g, h) leaf internal-to-external CO₂ ratio (χ); (i, j) electron-transport capacity (J_max, µmol CO₂ m⁻² s⁻¹); (k, l) soil-to-leaf water potential difference (Δψ_soil−leaf, MPa), (m,n) plant-level water use efficiency (pWUE, µmol CO₂ (mmol H₂O)⁻¹), (o, p) intrinsic water use efficiency (iWUE, µmol CO₂ (mol H₂O)⁻¹).

In comparison, soil drought and its impacts are barely evident in the five best-fit PPE members, as soil water potential in these simulations remains well above the wilting point (Fig. 9a, c), except for those from QV_swrc3 with a more drought-prone SWRC. When soil moisture approaches the wilting point, photosynthetic capacity, assimilation rate, and stomatal conductance (Fig. 9a, c, e, i) decrease drastically, whereas pWUE and iWUE increase with lower soil water potential (Fig. 9m, o). Such response contrasts with the effect of increasing VPD, which tends to promote iWUE but reduce pWUE (cf. Fig. 9m, o and n, p).

4 Discussion

4.1 Apparent hydraulic traits of boreal crops

Crop hydraulic traits have often been overlooked in boreal crop studies (Roitsch et al., 2022). By employing SPY-C and EC flux measurements, we were able to infer apparent hydraulic traits at field scale and assess their impacts on crop performance under typical boreal conditions. Constrained by both observed carbon and water fluxes, the maximum whole-plant hydraulic conductance (K_p) inferred for both the perennial forage grass mixture and the annual cereal (oats) was toward the lower end of values reported for crop species (cf., Tables 2 and S3), suggesting relatively low hydraulic efficiency and a conservative water-use strategy. This contrasts with the relatively high hydraulic efficiency reported for many crops (Table S3), but is consistent with observations of relatively lower hydraulic efficiency in plants from cold boreal environments (He et al., 2020). This similarity supports the role of environmental water stress in shaping hydraulic traits across species (Bartlett et al., 2016), implying that environmental stresses associated with the boreal climate may favour relatively conservative water-use strategies in boreal crop species.

The inferred apparent hydraulic traits of both crops show little variation when using different model-inversion periods (Methods S4), consistent with evidence that hydraulic traits can be relatively conserved within a species (Lamy et al., 2014). We also inferred a slightly more acquisitive water-use strategy (e.g., lower γ or higher K_p) for the annual cereal crop (oats) and cover crop (Italian ryegrass) than for the perennial forage grass mixture (timothy, meadow fescue, and white clover) (Table 2), which is in broad agreement with the empirical evidence summarized in Table S3.

To further evaluate the plausibility of the inferred apparent hydraulic traits, we compared the modelled responses of stomatal conductance and whole-plant hydraulic conductance to changes in leaf water potential with available empirical observations. As shown in Fig. 10, the model inferred plant hydraulic conductance (K) and stomatal conductance for the perennial grass mixture are generally lower than those reported empirically (Fig. 10a, d). Agreement improved at lower leaf water potentials. One possible reason for such discrepancy can be the difference in scale: model-inferred K reflects whole-plant hydraulic conductance at the field scale, whereas empirical studies typically measure hydraulic conductance at the organ or individual plant scale (Wang and Frankenberg, 2022; Jiménez-Rodríguez et al., 2024). Differences in species, cultivars, and climatic origin may also contribute, as the available measurements were all from temperate rather than boreal cultivars. Additional hydraulic measurements of Nordic crops are therefore needed to better constrain and evaluate the model-inferred hydraulic traits and behaviours.

Figure 10 Modelled response of stomatal conductance (1.6 × g_s, mol H₂O m⁻² s⁻¹) (a–c) and hydraulic conductance (K, 10⁻¹⁶ m) (d–f) to declining leaf water potential (ψ_plant, MPa) and comparison with observations for similar crops reported in the literature. Panels (a) and (d) are for forage grasses, (b) and (e) are for oats and (c) and (f) are for the cover crop (Italian ryegrass). The modelled results (small circles) correspond to the QV_swrc3 PPE members using the five best-fit parameter sets listed in Table 2 for each crop, with best-fit members 1–5 following the same order as in the table. We use QV_swrc3 here because it produces a wider range of leaf water potentials, allowing hydraulic responses to be more clearly evaluated. The empirical data are denoted by triangles, and were obtained from Holloway-Phillips and Brodribb (2011b) (Lolium perenne), Jacob et al. (2022) (Festuca arundinacea, Phalaris aquatica), Canales et al. (2021) (Avena sativa) and Holloway-Phillips and Brodribb (2011a) (Lolium multiflorum cv. Feast). Hourly model output from the growing season (May–October, 2018–2023) that matched the environmental conditions of the corresponding leaf gas-exchange measurements were selected. L. perenne and L. multiflorum: VPD 1.35–1.65 kPa, PAR 900–1100 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and temperature 18–22 °C; F. arundinacea and P. aquatica: VPD 1–2 kPa, PAR 1620–1980 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and temperature 26–30 °C; A. sativa: VPD 0.72–0.88 kPa, PAR 225–275 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and temperature 18–22 °C.

For oats and Italian ryegrass, the model generally captured the observed stomatal responses to declining leaf water potential (Fig. 10b, c), except that the observation exhibited more rapid stomatal closure during dehydration than predicted by the model for Italian ryegrass, indicating more isohydric behaviour. A higher ψ₅₀ (i.e., −1 MPa) produced a closer match to the observed response to leaf dehydration. As values above −1 MPa have recently been reported for some crops (Cardoso et al., 2025), increasing ψ₅₀ above −1 MPa may further improve agreement between the model and observations. Increasing γ (the cost of maintaining hydraulic pathways) also shifted the simulated response toward more isohydric behaviour (data not shown). These sensitivities highlight the need for direct hydraulic measurements of boreal crops to better constrain both ψ₅₀ and γ and to assess the remaining model–observation discrepancies.

4.2 Role of hydraulic coordination in crop productivity and drought resistance

Previous studies have widely reported positive relationships between hydraulic efficiency (as indicated by maximum plant hydraulic conductance) and productivity (Holloway-Phillips and Brodribb, 2011a; Ocheltree et al., 2016; Gleason et al., 2021). Such a relationship was also detected in our study, but only during cooler and wetter months, such as May (Fig. 7a, b). During summer, when drought events were often present (Fig. 2e, f), however, we found a pronounced negative relationship between hydraulic efficiency (K_p) and productivity (GPP) across the full PPE (Fig. 7a, b) in both perennial and annual boreal crops. Therefore, the low hydraulic efficiency inferred for the studied crops appears to be particularly important for maintaining high productivity under boreal summer conditions. This can be attributed to the fact that plants with lower hydraulic efficiency withdraw soil water more slowly, thereby preventing soil moisture from being depleted to a critical threshold at which productivity becomes strongly constrained (Figs. 6i, j and 7i, j). A similar mechanism has been reported in forest ecosystem model simulations (Jiménez-Rodríguez et al., 2024). Consistently, previous empirical studies also suggest that crops with a conservative water-use strategy (i.e., lower hydraulic efficiency) can sustain higher yields than those with an acquisitive water-use strategy (i.e., higher hydraulic efficiency) (Holloway-Phillips and Brodribb, 2011a; Gleason et al., 2022; Schell et al., 2026). Importantly, the negative effect of hydraulic efficiency on productivity largely disappeared when either soil or atmospheric drought was removed (Figs. S10 and S11), indicating that the advantage of low hydraulic efficiency for maintaining high productivity is manifested primarily when soil and atmospheric drought co-limit crop productivity.

Overall, for the five best-fit PPE members, our modelling results indicate that atmospheric drought exerted a much stronger effect on crop productivity than soil drought at the two study sites (Fig. 8). Little evidence of soil-water limitation was found in the best-fit members from the baseline PPEs (i.e., QV_ctrl and HA_ctrl), whereas substantial soil drought effects emerged only in the best-fit members of the PPE with a more drought-prone soil water retention curve (i.e., QV_swrc3) (Fig. 9). The limited role of soil drought was further supported by the absence of a critical soil-moisture threshold in either the flux observations or the best-fit PPE members from QV_ctrl and HA_ctrl (Fig. S8b–d), following the diagnostic proposed by Fu et al. (2024). The hydraulic safety margin of the studied crops, defined as the difference between the minimum leaf water potential experienced during the growing season (ψ_min) and ψ₅₀ (Choat et al., 2012; Anderegg et al., 2016), remained well above zero despite the large uncertainties in ψ₅₀ (Fig. S9). Together, this evidence suggests that although soil water stress was indeed a potential constraint on crop productivity during drought events under current boreal climate conditions, the inferred conservative hydraulic strategies of the studied crops appear to have helped them largely avoid this constraint, thereby exhibiting high resistance to the exposure of drought events primarily through drought avoidance rather than tolerance (Haugum et al., 2021; Alves et al., 2024; Pereira et al., 2024; Cardoso et al., 2025). We note that the current version of SPY-C does not account for potential impacts of soil drought on above- and below-ground carbon allocation, which have been widely recognised in previous studies (Qiao et al., 2021; Xia et al., 2017). Whether such responses occurred at our study sites remains to be investigated. Incorporating a dynamic carbon allocation scheme responsive to soil water availability would help address this limitation.

In contrast to the limited effect of soil drought, atmospheric drought exerted a clear and consistent constraint on crop productivity at both sites (Figs. 8c, d and 9). The extensive periods of high VPD may have contributed substantially to the yield reductions observed during drought years at Qvidja (Heimsch et al., 2021), consistent with evidence for a dominant role of atmospheric dryness in regulating plant productivity in boreal regions (Mirabel et al., 2023). The model inferred relatively low stomatal conductance (Fig. 9c, d) and a faster decline in leaf water potential than in soil water potential (Fig. S8a), an indicator of anisohydric stomatal behaviour (Hochberg et al., 2018; Ratzmann et al., 2019). Stomatal closure generally occurred only after significant loss of hydraulic conductivity (Fig. S9), implying a relatively “risky” stomatal regulation and a low stomatal safety margin (Gleason et al., 2022; Martin-StPaul et al., 2017; Pereira et al., 2024). These model-inferred behaviours broadly concur with empirical evidence from perennial forage grasses (e.g. Holloway-Phillips and Brodribb, 2011b; Jacob et al., 2022) and some cereal crops (e.g. Corso et al., 2020), and may explain the strong sensitivity to atmospheric drought found in our simulations. The relatively low hydraulic efficiency, however, appeared to partly buffer the crops against atmospheric drought, as evidenced by a positive effect of hydraulic efficiency (K_p) on the response of GPP to the removal of VPD (Fig. 8c, d). The buffering was nevertheless insufficient at daily VPD ≥ 0.8 kPa, when hydraulic adjustment could no longer sustain leaf water balance and both stomatal conductance and photosynthesis declined (Fig. 9b, d, f, j). This is in agreement with the close coordination between plant hydraulic function, stomatal regulation, and photosynthesis under drought conditions discussed by Flexas et al. (2018).

There has been controversy over how plant hydraulic traits and their coordination affect its water use efficiency. While high (low) hydraulic efficiency (safety) has often been suggested to promote plant-level WUE (Holloway-Phillips and Brodribb, 2011b; Gleason et al., 2022), low hydraulic efficiency has also been found to promote intrinsic water use efficiency for alpine shrubs (Yao et al., 2024). Our results showed that the effect of hydraulic efficiency (K_p) on plant water-use efficiency (both pWUE and iWUE) was generally opposite to its effect on GPP (cf. Fig. 7a, b and e–h). Accordingly, hydraulic strategies that favoured GPP did not necessarily maximize pWUE or iWUE (Fig. 6), especially under drought conditions. In our simulations, higher hydraulic efficiency (K_p) tended to increase pWUE and iWUE during summer water deficit but decrease it during wetter spring and autumn conditions (Fig. 7e–h). This is in agreement with the contrasting effects of hydraulic efficiency reported in a previous study on annual ryegrass (see Table 4 in Holloway-Phillips and Brodribb, 2011a). Under summer drought, plants with higher hydraulic efficiency (K_p) may accelerate soil-water depletion and lower soil and plant water potentials (Fig. 11a), giving rise to a substantial decline in realized plant hydraulic conductance (K) (Eq. 3). This reduction in realized conductance decreases both GPP and transpiration, but transpiration declines more strongly than GPP, thereby increasing pWUE (Fig. 11b–d). In contrast, under wetter spring and autumn conditions, plants with higher hydraulic efficiency can sustain higher realized hydraulic conductance, and thus exhibit lower pWUE than their less hydraulically efficient counterparts.

Figure 11 Relationship between whole-plant hydraulic conductance (K, 10⁻¹⁶ m) and soil water potential (ψ_soil, MPa) (a), and between the ratio of whole-plant hydraulic conductance and the corresponding ratios of GPP (b), transpiration (c), and pWUE (d). The ratios were calculated between simulations with higher K_p and the simulation with the lowest K_p (0.5 × 10⁻¹⁶ m). All results are from the QV_noVPD PPE. The intersection of the horizontal and vertical grey dashed lines in panels (b), (c), (d) denotes the reference point where both ratios equal 1.

4.3 Model uncertainties and future development

Several structural simplifications in the current version of SPY-C may contribute to uncertainties in the simulated carbon and water fluxes and the inferred apparent hydraulic traits. First, SPY-C adopts a big-leaf canopy representation, in which prescribed total one-sided green LAI is combined with within-canopy light attenuation using Beer's law. The model therefore does not explicitly distinguish sunlit and shaded leaves or resolve multiple canopy layers. This simplification may contribute to biases in the diurnal partitioning of GPP and ET, particularly during the early morning and late afternoon (Fig. 5). A two-leaf or multilayer canopy representation could improve the treatment of within-canopy radiation and associated carbon and water fluxes.

A second important simplification concerns the representation of roots dynamics and its influence on whole-plant hydraulics. In the current framework, K_p represents the maximum whole-plant hydraulic conductance per unit leaf area and therefore integrates the entire soil-to-leaf hydraulic pathway, including roots and the root–soil interface. Root hydraulic traits are thus represented only implicitly rather than as separate root, stem, and leaf conductances. Moreover, water uptake occurs from a fixed effective root zone (0.6 m in this study), while root biomass is tracked as a carbon pool but does not feed back on rooting depth, water uptake, or hydraulic conductance. Dynamic rooting depth, root distribution, hydraulic redistribution, and drought-induced changes in root hydraulic conductivity are therefore not represented. These processes can play an important role in drought adaptation of cereals and grasses (e.g. Canales et al., 2021). Coupling dynamic root growth to water uptake and hydraulic conductance, together with prognostic leaf growth, is therefore an important direction for future development. Such extensions would also allow drought avoidance mediated through root-zone dynamics to be represented more explicitly.

Uncertainty also arises from model structural biases, parameterization, and model inversion approach. We note indications that canopy and ground evaporation may be overestimated leading to the overestimation of ET (Fig. 4d, h). However, sensitivity experiments with reduced evaporation did not alter the inferred K_p (data not shown), suggesting that this bias has little influence on the inferred apparent hydraulic traits. Although the PPE approach enabled us to constrain apparent hydraulic traits using both carbon and water fluxes, the present ensembles varied only a limited set of plant hydraulic and photosynthetic parameters. Future applications should propagate uncertainty from a broader range of parameters, particularly soil hydraulic properties, rooting characteristics, and canopy processes. More efficient multi-objective calibration and inversion approaches (Denager et al., 2023; Bell et al., 2025) may further improve both parameter identifiability and robustness. An additional source of uncertainty also lies in the representation of mixed cropping systems. At Qvidja, the forage grass and legume species are represented as a single effectively homogeneous canopy, although individual species may differ in hydraulic traits, stomatal regulation, and drought response. Resolving such within-canopy hydraulic heterogeneity would therefore be a valuable direction for future development.

Despite these uncertainties, the strong influence of plant hydraulics on crop productivity, water use, and drought responses demonstrated here highlights the value of explicitly representing hydraulic processes in modelling boreal crops. SPY-C also reproduced net ecosystem exchange reasonably well at the study sites (Fig. S12), suggesting potential for applications in monitoring, reporting, and verification (MRV) of agricultural carbon balances (Smith et al., 2020; Guan et al., 2023; Wijmer et al., 2024). Developing SPY-C toward a robust MRV tool will nevertheless require a more complete representation of management and nutrient constraints. In particular, improved descriptions of sowing, harvesting, fertilization, and crop rotation, together with coupling of nitrogen dynamics to the carbon and water cycles (e.g. Stocker and Prentice, 2024), would allow crop growth and carbon–water exchange to respond more realistically to management. Combined with improved canopy, root, and parameter-inference schemes, these developments would strengthen the transferability of SPY-C across crop types, sites, and management regimes.

5 Conclusions

We examined the apparent hydraulic traits of typical annual and perennial boreal crops and their impact on crop productivity, water use efficiency, and drought response by combining multi-year carbon and water flux measurements from two Finnish crop fields with the newly developed model framework SPY-C, which explicitly couples plant hydraulic, photosynthesis, and soil moisture dynamics. SPY-C reproduced the observed carbon and water fluxes well to support field-scale inference of apparent hydraulic traits. Both oats and the perennial forage grass mixture exhibited relatively low apparent hydraulic efficiency and conservative hydraulic strategies, which are crucial for reducing rapid soil-water depletion and partly buffered crop functioning against atmospheric drought, thereby promoting drought avoidance and productivity under boreal summer conditions. Here, we present SPY-C as an efficient tool for inferring apparent hydraulic traits and assessing their effects on crop performance across contrasting pedoclimate conditions. This can be fundamental for selecting and developing better crops and management to support climate-smart and sustainable agriculture. We also underscore the good potential of SPY-C for broader applications in agricultural carbon monitoring, reporting, and verification to enhance accuracy and predictive power in both current and future boreal climates.