A model bridging waterlogging, stomatal behavior and water use in trees in drained peatland

Abstract Waterlogging causes hypoxic or anoxic conditions in soils, which lead to decreases in root and stomatal hydraulic conductance. Although these effects have been observed in a variety of plant species, they have not been quantified continuously over a range of water table depths (WTD) or soil water contents (SWC). To provide a quantitative theoretical framework for tackling this issue, we hypothesized similar mathematical descriptions of waterlogging and drought effects on whole-tree hydraulics and constructed a hierarchical model by connecting optimal stomata and soil-to-leaf hydraulic conductance models. In the model, the soil-to-root conductance is non-monotonic with WTD to reflect both the limitations by water under low SWC and by hypoxic effects associated with inhibited oxygen diffusion under high SWC. The model was parameterized using priors from literature and data collected over four growing seasons from Scots pine (Pinus sylvestris L.) trees grown in a drained peatland in Finland. Two reference models (RMs) were compared with the new model, RM1 with no belowground hydraulics and RM2 with no waterlogging effects. The new model was more accurate than the RMs in predicting transpiration rate (fitted slope of measured against modeled transpiration rate = 0.991 vs 0.979 (RM1) and 0.984 (RM2), R2 = 0.801 vs 0.665 (RM1) and 0.776 (RM2)). Particularly, RM2’s overestimation of transpiration rate under shallow water table conditions (fitted slope = 0.908, R2 = 0.697) was considerably reduced by the new model (fitted slope = 0.956, R2 = 0.711). The limits and potential improvements of the model are discussed.


Introduction
Waterlogging, or flooding in soil characterized by high soil water content (SWC) and shallow water table depth (WTD), drives oxygen from soil pores and leads to hypoxia or anoxia (i.e. insufficient or absent oxygen) in the belowground parts of plants. This phenomenon inhibits aquaporin activity and mitochondrial respiration and causes cytosol acidosis in the roots of nearly all trees, thereby significantly decreasing root hydraulic conductance (Lambers et al. 2008, Parent et al. 2008, Kreuzwieser and Rennenberg 2014. Following the inhibited water uptake in roots, a variety of tree species showed stomatal closure in waterlogged conditions (Pezeshki and Chambers 1986, Pezeshki et al. 1996, Ferner et al. 2012. This consequent decrease of stomatal conductance occurs likely because of whole-plant hydraulic signaling for maintaining water homoeostasis (Jackson 2002, Christmann et al. 2013 and is more severe in floodsensitive individuals or species (Pezeshki et al. 1996, Pearson et al. 2013, Kreuzwieser and Rennenberg 2014. Ultimately, severe symptoms may occur, including leaf necrosis, bark shedding and whole-tree dieback, and the tree's production and growth (above-and belowground) may be strongly impaired (Parolin 2001, Alaoui-Sossé et al. 2005, Ferry et al. 2010, Kreuzwieser and Rennenberg 2014, Repo et al. 2020. Despite their wide occurrence and drastic consequences, the waterlogging effects on whole-tree hydraulics have not been clearly quantified due to difficulties in conducting continuous observations of their whole hydraulic system. Nevertheless, the similar responses of stomatal behavior to excess soil water and to air/soil dryness suggest that the existing models of root and stomatal conductance may be modified to represent waterlogging effects. More specifically, the recent observation of the similar impacts of drought and flooding on root conductance has suggested the possibility of simulating the effects of both disturbances in similar mathematical structures (Domec et al. 2021). Such a model, once verified by data, should describe the soil-tree-atmosphere continuum more accurately on a wider range of water conditions, and shed light upon whole-tree hydraulics and leaf gas exchange responding simultaneously to above-and belowground factors in a quantitative manner.
Among the models of stomatal behavior in relation to air dryness, irradiance and temperature, the Lagrangian optimization model is one of the most widely tested and successful (Cowan 1977, Cowan and Farquhar 1977, Hari et al. 1986, Hari and Mäkelä 2003. The core of this type of model is that the steady state of stomatal behavior realizes the optimization (maximization) of the integrated difference between carbon gain and water cost over a given time, i.e. max A(t) − λE(t) dt (Hari et al. 1986), which yields the marginal water-use efficiency (MWUE, λ) at the steady state. In later studies, the model has been developed to cover more physiological processes and/or wider scales (Dewar 2002, Katul et al. 2010, Medlyn et al. 2011, Manzoni et al. 2013, Prentice et al. 2014, Bell et al. 2015 and used for quantifying specific ecophysiological issues, such as effects of severe drought, trees' competition for water and global functional diversity of plants (Lin et al. 2015, Wolf et al. 2016, Lu et al. 2020. Nonetheless, the relationship between soil water and MWUE in the related variants of the model is always monotonic, and excess water (i.e. waterlogging or flooding) turning into a constraint to soil-to-root and stomatal conductance has not been incorporated. Not only is such an incorporation important for quantifying the ecophysiological processes, but it should also improve our ability to account for the effects of increasing extreme events of precipitation on stomatal behavior in largescale predictions (Palmer and Räisänen 2002, Christensen and Christensen 2003, Frei et al. 2006, IPCC 2021. Considering such potential benefits and previous studies, we hypothesized that the impact of waterlogging in decreasing trees' water uptake rate is similar in magnitude to that of air/soil dryness. In other words, similar algebraic expressions should be employed to describe the effects of both waterlogging and air/soil dryness on soil-to-root conductance, water use and stomata behavior. To test the hypothesis, we linked the Lagrangian optimal stomata model with a soil-to-leaf conductance model (Duursma et al. 2008, Nikinmaa et al. 2013, Hölttä et al. 2017, Dewar et al. 2018, which was modified to describe the waterlogging effects on water use in trees. This bridging was constructed by writing MWUE, instead of a single parameter λ, as a segmented (and thus 'non-monotonic') function of SWC and soil-to-leaf conductance affected by WTD. The underlying assumption of the decreased conductance was that the inhibited oxygen diffusion under hypoxic conditions reduces soil-to-root conductance, which was implemented in an algebraically simplified manner. The model was tested using data of sap flow density in Scots pine (Pinus sylvestris L.) trees, WTD and other environmental factors collected from a subarctic drained peatland site in the Finnish Lapland over four growing seasons.

Study site and data acquisition
The study site was in Sattasuo (66 • 30 N, 26 • 42 E, 165 m above sea level), a subarctic flat peatland drained by open ditch in the Finnish Lapland. The mean air temperature in July over 1981-2010 was 15.0 • C, and mean precipitation sum 505 mm per year (Stenberg et al. 2018). Typically, the ground was snow-covered from October through May. The first drainage was conducted in 1934, and the latest one before the studied period (2008-11) was in 2006. In 2006, the standing volume of the dominant Scots pine (Pinus sylvestris L.) was 91 m 3 ha −1 . The site was classified as medium in productivity (Vasander andLaine 2008, Stenberg et al. 2018).
Six naturally grown Scots pine trees (aged 66-87 years, counted at the breast height) were chosen as sample trees, of which two (Nos 2 and 3) were noticeably taller than the other four (Table 1). To determine their sapwood thickness without damaging them, another eight Scots pine trees of similar size at the site were bored, and sapwood thickness ϑ SW , mm measured on the cores was correlated with diameter at breast height (DBH, d, mm) as (1) ϑ SW of the six sample trees was calculated by Eq. (1) ( Table 1). Their sap flow was measured using thermal dissipation probes (TDPs; Granier 1987) through the unfrozen season (May to October) of 2008-11. A pair of TDP, comprising a heated (HP) and a reference (RP) probe, both c. 3.5 cm long, were installed on the north side of each sample tree's stem at the height of c. 1.3 m. The HP was heated with a constant power of 0.2 W (Lu et al. 2004), and the RP was directly below HP at a distance of c. 10 cm. All the probes were covered by aluminium foils with ventilation holes. The voltage difference between HP and RP ( U) was measured every 30 s, recorded as 10-min average, and converted to sap flow density using Baseliner 4.0 (Oishi et al. 2016). The sap flow signals of trees nos 2 and 5 were extremely noisy since June 2011, and thus their data were discarded thereafter. On a forest clearing 190 m north from the study site, air temperature, humidity (CS215 probe, Campbell, UT, USA) and photosynthetic photon flux density (PPFD; SP1110 pyranometer, Campbell) were recorded every 10 min. The vapor pressure deficit (VPD, in Pa) was calculated as (Campbell and Norman 1998) where T is air temperature ( • C), and h r relative humidity. For the consistency with the process model (Eq. (4)), VPD was converted to the unit of mol H 2 O m −3 using the ideal gas law: To obtain the daily-level data for the subsequent analyses, the median of the largest 10% data by ranking of each day was gathered for each of the variables. Evenly distributed in the site, 10 1-m-long tubes were mounted into the ground for manually measuring water table depth (WTD). The manual measurement was conducted at an interval of 1-2 weeks through the same period as of the other data. A water height meter (Intech Instruments ltd New Zealand) was installed at the approximate center of the site, recording WTD and water temperature (WTT) every hour. The meter's temporal pattern of WTD was used for interpolating between the manual measurement data points to obtain the hourly WTD at each measurement tube. As for the other variables, the median of the largest 10% data by ranking of each day was taken as the daily WTD for the subsequent analyses.

The hierarchical model
A hierarchical model was constructed, which comprises two levels, namely, process and data models. The process model is an optimal stomata model associated with a module on the dependence of soil-to-root conductance on WTD via SWC. The data model contains the probability density function (PDF) of the error between modeled (E (M) ) and observed (E (O) ) transpiration rate.
Optimal stomata model In the optimal stomata model, the stomatal behavior maximizes the difference between the integrated carbon gain and water loss over a given time (Cowan andFarquhar 1977, Hari et al. 1986). At the steady state, transpiration rate (E (M) ) is modeled proportional to the stomatal conductance for water (g σ , m s −1 ), which in turn is dependent on VPD (D), respiration rate (R, mol CO 2 m −2 s −1 ) and PPFD (I, mol m −2 s −1 ) where λ is MWUE (mol CO 2 mol −1 H 2 O), C a the atmospheric CO 2 concentration (mol m −3 ), 1.6 the ratio of water vapor to CO 2 diffusion rates, and R and f (I) (irradiance reaction curve), respectively, are modeled as follows Mäkelä 2003, Mäkelä et al. 2004): where T l is leaf temperature (Eq. (9)), Q 10 relative increase of R per 10 • C, R 0 the value of R at 0 • C; and where γ is the saturation (asymptote) of the curve (m s −1 ), and ι (initial slope, m 3 mol −1 ) is where c is the slope coefficient (m 3 (mol • C) −1 ), and S is the photosynthetic acclimation of foliage to temperature above the threshold (S 0 ), of which the marginal change with respect to time (day) follows: where τ is the time constant (day) of the acclimation, and T l , leaf temperature, was calculated from air temperature (T, • C) and as Tree Physiology Volume 42, 2022 (Kolari et al. 2007), where I is as in Eqs (4) and (6). The 'S'model (Eqs (6-8)) has been widely applied to physiological and growth acclimation to temperature and irradiance at a range of temporal and spatial scales (Suni et al. 2003, Mäkelä et al. 2004, Mäkelä et al. 2008, Schiestl-Aalto et al. 2015. Water use and soil-to-leaf conductance Different from evaluating MWUE as a parameter (λ in Eq. (4)) directly (Hari and Mäkelä 2003, Mäkelä et al. 2004, Liu et al. 2020), we modeled λ as a function of soil-to-leaf conductance (k sl , mol H 2 O m −2 leaf s −1 Pa −1 ), which has two components in series, i.e. soil-to-root (k sr ) and root-to-leaf (k rl ) conductance, following When stomatal aperture is optimized and only constrained by nonstomatal factors (viz. carboxylation capacity and/or mesophyll conductance), k sl and λ are inversely correlated (Dewar et al. 2018), which can be expressed empirically as , and z 0 and z 1 are parameters evaluated in this study. The tree-specific root-to-leaf conductance (k rl in Eq. (10)) was assumed constant through the study period and estimated as where max J ij is the maximum sap flow density of tree i observed over the study period (i.e. data points j; Table 1), ρ leaf-to-sapwood area ratio, and lmin the absolute value of the lowest leaf water potential (Duursma et al. 2008). The trees were assumed isohydric through the study period, and thus lmin was assumed constant at −2 MPa, a typical value for Scots pine (Martínez-Vilalta et al. 2009). The other component of k sl , k sr (Eq. (10)) was designed in a segmented manner to reflect its non-monotonic dependence on soil water content (SWC, θ, m 3 m −3 ).

Soil-to-root conductance and water table depth
We segmented the soil-to-root conductance (k sr in Eq. (10)) into increasing, stable and decreasing phases with respect to SWC (calculated from WTD; Eq. (17)) to reflect the switch from water (deep WTD) to oxygen (shallow WTD) limitation via a relatively stable stage. In the increasing phase of SWC, k sr is a power function of SWC relative to its saturation (θ sat ), i.e.
where ξ m (mol m −2 s −1 Pa −1 ) is a parameter depending on root length index (m root m −2 ground surface), rhizosphere and fine root radii, and saturated soil conductivity, and ξ p is related to the soil water retention curve (Duursma et al. 2008, Nikinmaa et al. 2013). This form of correlation is also supported experimentally (Poyatos et al. 2018) at low to medium SWC. We were not aware of any experimental studies on the responses of mature Scots pine to excess water featured with a gradient of SWC or WTD, of which the results would help hypothesize a continuous k − sr for the decreasing phase. However, the flooding-sensitive loblolly pine (P. taeda L.) has shown similar declines in root conductance under drought and flooding treatments compared with control (Domec et al. 2021). Thus, we assumed k − sr θ to be also a power function with the base modified to represent the available soil porosity not occupied by water (i.e.1 − θ/θ sat ). Furthermore, to utilize prior knowledge on the optimal WTD (δ * ; Hökkä et al. 2021), we introduced the optimal SWC (θ * ) corresponding to δ * (Eq. (17)) as a parameter. Hence, where η m and η p are parameters. This formula resembles the soil-type-specific diffusivity coefficient of oxygen (Moldrup et al. 1996(Moldrup et al. , 1997. The optimal conductance k * sr was assumed to be maintained through a range of 20 cm of WTD between the increasing and decreasing phases, i.e. where θ • is SWC as a function of WTD (Eq. (17)). Finally, over the range of the WTD data, Throughout this study, all values of θ were calculated from WTD (δ, cm) by where θ res is the residual SWC (Leppä et al. 2020). From Eqs (4) through (17), transpiration rate and stomatal conductance are linked with WTD via MWUE in trees ( Figure 1).

Data model and parameterization
The data model or the PDF of the error between tree-specific observed and modeled transpiration rate is Tree Physiology Online at http://www.treephys.oxfordjournals.org  Table 2 for meanings and prior ranges of the evaluated parameters, and where i is tree number, j the index of data points in the time series, J ij and ρ same as in Eq. (12), and E (M) ij follows Eq. (4). Judged by the initial runs using the parameter priors, the error ε ij calculated by Eq. (18) was assumed to follow the normal distribution, of which the standard deviation (SD) parameter was correlated with E (M) for the observed heteroscedasticity. This correlation was exponential (i.e. semi-log linear) in two trees (nos 4 and 6, noted collectively W 1 ) and simple linear without significant intercept in the other four (no. 1-3 and 5, noted collectively W 2 ). Thus, the likelihood function of ε is where ε is the error matrix ε def = ε ij , α and β, respectively, are the intercept (if used) and slope coefficients in the semi-log or simple linear regression for the SD parameters, i and j as in (Eq. (18)), N i the total number of data points of tree i, and the subscripts denote the ordinal numbers of the trees.
As k + sr has been found correlated with tree height (Martínez-Vilalta et al. 2007a), its multiplier coefficient (ξ m in Eq. (13)) was tree specific, each tree still sharing the same prior range (Table 2). Considering the possible difference among trees in the acclimation to temperature regarding reaction to irradiance (Eq. (7)), c was also tree specific with the same prior range. All the other parameters were shared by all the trees, including η m and η p (Eq. (14)) due to limited observed evidence. The model had 22 estimated parameters (including 19 for the process model; Table 2) and four direct input variables D, I, T l and δ. All the parameters were assumed independently and uniformly distributed. The full list of symbols is attached in Tree Physiology Volume 42, 2022 Intercept (z 0 ) and slope (z 1 ) of the log-log linear regression correlating marginal water-use efficiency and soil-to-leaf conductance  Table 1). 3 ξ mi and ξ p are products of multiple parameters in Refs 2 and 3, and the reference range of each component parameter was used for calculating the priors of ξ mi and ξ p with adjustment based on Ref. 1 specifically on peatland. The prior ranges of η m and η p were similar to those of ξ mi and ξ p , respectively, but slightly broadened due to lack of direct observation in literature. 4 c is denoted c 1 in Ref. 9, which estimated its values using shoot-level data. Considering the possible heterogeneity at whole-crown level, we broadened its prior range for the current study.
Supplementary data (Table S1 available as Supplementary data at Tree Physiology Online). The maxima a posteriori (MAP) estimates of the parameters were sought by adaptive Markov chain Monte Carlo (MCMC) algorithm DREAM (ZS) (Vrugt et al. 2009) in R 4.0.3 (R Development Core Team 2020) with package 'BayesianTools' (Hartig et al. 2019). Convergence was defined atR < 1.1, and only the second halves of the chains (i.e. without 'burn-in') were used for the convergence diagnostic and the subsequent analyses (Gelman et al. 2013). The 95% credible interval of predictive uncertainty was generated by Eq. (19), defined as the interval between the 2.5 and 97.5% quantiles of ε calculated with 6000 parameter vectors randomly sampled from the posterior distribution.

Model performance assessment
To assess the improvement brought about by the new module linking MWUE to WTD with waterlogging effects (Eqs (10-17)); hereafter termed full model, FM), two reference models were calibrated for comparison: (i) the base optimal stomata model (Eqs (4-9)) with λ estimated directly as a temporally constant tree-specific parameter (RM1), and (ii) a bridged model including the soil-to-leaf conductance but with a monotonic correlation between k sr and SWC (i.e. without considering waterlogging effects; RM2). The prior range of λ in RM1 was according to previous studies by similar models on Scots pine (Hari and Mäkelä 2003, Mäkelä et al. 2004, Liu et al. 2020). The heteroscedasticity descriptions in the data models (Eq. (19)) of RMs were adjusted accordingly. Otherwise, FM and RMs shared the same structure and prior ranges of parameters, MCMC method and input data.
The models' overall performances were judged by linear regression with zero intercept of E (O) to E (M) computed with the parameters' MAP estimates (hereafter simply E (M) ). To compare the performances of FM and RM2 particularly on shallow WTD, this linear regression was performed with E (M) using the full data as well as those in the range of k − sr (corresponding to WTD < δ * − 10 cm; Eq. (15)).
The tree-specific performance of FM was assessed by the normalized root-mean-square error (NRMSE) of the corresponding tree, defined as denotes the mean of E (O) of tree i over the study period, and the rest of the notation is the same as in Eqs (18) and (19). To clarify the possible causes of the residual, correlation coefficient (Pearson's r) was calculated between FM's residual (ε = E (O) − E (M) , mol H 2 O m −2 leaf s −1 ) and WTT. The WTT was measured within the same tubes for the manual measurement of WTD. All the statistical analyses were performed in R 4.0.3.

Results
With the MAP estimates of the parameters (Table 3) and the full data, FM performed well in simulating transpiration rate, presenting a fitted slope of 0.991 (R 2 = 0.801, Figure 2). In contrast, the fitted slope for RM1 (0.979) deviated more from 1, Tree Physiology Online at http://www.treephys.oxfordjournals.org Table 3. MAP estimates and 95% credible intervals of the estimated parameters of the process model. Detailed information of the parameters (including their prior ranges) is shown in and its R 2 was lower (0.665, Figure 3A). The performance of RM2 (slope = 0.984, R 2 = 0.776) was close to that of FM when using the full data, but on shallow WTD (< δ * − 10 cm = 37.98 cm; Table 3) RM2's overestimation of transpiration rate (fitted slope = 0.908, R 2 = 0.697; Figure 3B) was larger than FM's (fitted slope = 0.956, R 2 = 0.711; Figures 2). All the fitted slopes were significant (P < 0.001). The tree-specific NRMSE of FM was low to moderate, ranging from 20.98 to 35.34% (Figure 4). The optimal WTD (δ * , i.e. the WTD through the stable phase of k sr θ ) was estimated to be 46.9 cm (Table 3), and thus, a decline in k sr occurred from WTD = 33.4-39.3 cm (varying across trees) through shallower WTD levels ( Figure 5). Correspondingly, the simulated maximum k sr at the stable phase ranged between 2.72 × 10 −9 and 3.55 × 10 −9 mol H 2 O m −2 leaf s −1 Pa −1 (Figure 5). The simulated λ at δ * ranged between 3.94 × 10 −3 and 5.73 × 10 −3 mol CO 2 mol −1 H 2 O ( Figure 6). The values of λ in all the trees increased sharply in waterlogged conditions (i.e. during the decrease of k sr θ , Eq. (14), Figure 5), while a moderate rate of increase was found when WTD was deep. As a result, similar MWUE values were yielded at WTD = 0.1 m and WTD ≈ 1.5 to 2 m (varying across trees; Figure 6). The two taller trees (Nos 2 and 3, Table 1) had the highest λ values (Figure 6), but their lowest k sr values did not differ notably from those of the other trees ( Figure 5).
Of the full data, FM's residual was significantly but not strongly correlated with WTT (r = 0.223), and the correlation was higher when WTT ≤ 5 • C (r = 0.224, Figure 7). Under shallow WTD (< 36.9 cm) and WTT ≤ 5 • C, the residual and Figure 2. Overall performance of the full bridged model with waterlogging effects (FM). Lines without intercept were fitted to the data of observed transpiration rate (E (O) , converted from measured sap flow density per sapwood area; y) and simulated results (E (M) ; x). Dark gray and dark blue lines, respectively, denote the fittings using the full data (light gray circles) and those with shallow WTD (< δ * − 10 cm = 37.98 cm; Table 3; blue circles). Their respective equations are y = 0.991 x and y = 0.956 x, R 2 = 0.801 and 0.711, and both slopes' P < 0.001. The dashed red line marks the slope of 1.
WTT had r = 0.234, but on the full WTT range the correlation was low (r = 0.215, Figure 7). In all cases, correlations were highly significant (P < 0.001).
Tree Physiology Volume 42, 2022 base optimal stomata model (Eqs (4-9)) with λ directly estimated without soil-to-leaf conductance, and (B) a bridged model including soil-to-leaf conductance with a monotonic soil-to-root conductance (i.e. without accounting waterlogging effect). Dark gray and dark blue lines, respectively, denote the fittings using the full data (light gray circles) and those with shallow WTD (< δ * − 10 cm = 37.98 cm; Table 3; blue circles). The data set was not so separated for RM1 as it lacks the belowground structure. The dark gray fitted line in (A) is y = 0.979 x (R 2 = 0.665), and the dark gray and dark blue fitted lines in (B) are, respectively, y = 0.984 x and y = 0.908 x, R 2 = 0.776 and 0.697. All slopes' P < 0.001.

Discussion
A soil-to-leaf hydraulic conductance model was modified to reflect waterlogging effects, linked with stomatal optimization model, and parameterized with measured data from a drained peatland. Compared with the reference models (RM1 and RM2), the full model (FM) using the full data presented an E ( Figure 2). FM was also accurate at individual tree level, presenting root-mean-square errors c. 1/5 to 1/3 of the observed means of respective trees (Figure 4). Thus, the hypothetical model structures of non-monotonic k sr (Eqs (13-16), Figure 5) and λ k sl (Eq. (11), Figure 6) gained support from the data and the comparison with RMs without such designs.

Model structure
The adverse effects of excess soil water on hydraulics and gas exchange in plants vary across taxa and soil types (Pezeshki and Chambers 1986, Oren et al. 2001, Ferner et al. 2012, Pearson et al. 2013, and we lack experimental observations of such effects under a gradient of SWC or WTD for continuous quantification. Thus, no precedent process-based models tackled this problem before the current attempt. Nevertheless, there has been consensus that hypoxia or anoxia is the key phenomenon that strongly alters the physiochemical properties and processes in soil related to plant physiology (Kreuzwieser and Rennenberg 2014). Furthermore, the consequent changes following hypoxia in the number and/or the status of aquaporins have been found to be generally the pivotal factor impacting soil-to-root conductance (Clarkson et al. 2000, Chaumont et al. 2005, Lambers et al. 2008, Johnson et al. 2014, Vandeleur et al. 2014, Domec et al. 2021. Particularly in flooding-sensitive pine species, it is likely that the decline of aquaporin activity under flooding is similar to that under drought (Domec et al. 2021). Therefore, we assumed the decreasing phase of soil-to-root conductance k − sr to occur due to excess soil water, which is related to the oxygen availability controlled by the available soil air space and diffusion (Eq. (14)). Despite the algebraically simplistic design, this hypothetical structure successfully simulated the transpiration rate of Scots pine trees in a drained peatland over four growing seasons (Figures 2  and 4). Also, the performance was better than that of RMs representing (i) stomatal behavior influenced by D, I and T l alone or (ii) with monotonic correlation between SWC and k sr as in previous studies (Duursma et al. 2008, Nikinmaa et al. 2013, Hölttä et al. 2017, Figure 3). With regard to these results, we Tree Physiology Online at http://www.treephys.oxfordjournals.org conclude that the new model corroborates the prior findings on hypoxia and aquaporins in the perspective of soil-treeatmosphere continuum and provides new quantitative insights into the water uptake and use of trees under waterlogged conditions in drained peatland.

Soil-to-root conductance, water use and water table depth
The maximum simulated k sr ( Figure 5) did not show significant negative correlation with tree height (correlation coefficient r < 0 but P > 0.5), in contrast with the belowground conductance significantly declining with tree height in previous Tree Physiology Volume 42, 2022 Figure 5. Simulated soil-to-root conductance (k sr , Eqs (13-16)) in relation to water table depth (WTD). Black lines are k sr using the MAP estimates of the related parameters. From high to low at the maximum k sr ('plateau') are trees Nos 4 (dot-dashed), 6 (long-dot-dashed), 5 (long-dashed), 2 (dashed), 3 (dotted) and 1 (solid line). Red solid segments demarcate the range of measured WTD. Figure 6. Simulated marginal water-use efficiency (λ, Eqs (10-16)) in relation to water table depth (WTD). Red solid segments demarcate the range of measured WTD. The parameterization of the simulation and the legend of lines (with tree numbers marked) are the same as in Figure 5.
observations (Martínez-Vilalta et al. 2007a). This distinction might be related to the different modeling method and more homogeneous sample trees in the current study. Instead of the simple one-membrane two-compartment (above-and belowground) model in the reference study (Martínez-Vilalta et al. 2007a), the current model employed a root structure module with more details, which was compacted into the treespecific parameter ξ mi (Eq. (13)). Ecophysiologically, ξ m = R L 2π ln r cyl /r root K sat , where R L is root length to all-sided leaf area ratio (m root m −2 leaf), r cyl /r root the ratio of hydraulically active rhizosphere cylinder to average fine root radii (both in m), and K sat the saturated conductivity of the soil (mol H 2 O m −1 s −1 Pa −1 ) (Nikinmaa et al. 2013). As the sample trees grew . The gray box refers to WTT ≤ 5 • C. cyan crosses and black circles, respectively, indicate the data points within and outside the range of decreasing soil-to-root conductance (k − sr θ , Eqs (14) and (16)). The color legend of data points (cyan and black) is the same as that of the Pearson correlation coefficient (r), which are labeled on the bottom right of the box (using the data with WTT ≤ 5 • C) and the panel (using the full data). All P < 0.001. closely in one stand with narrow height and age spans, it is reasonable to assume that none of K sat , R L or r cyl /r root varied significantly across the trees. In uneven-aged stands, r cyl /r root may particularly be positively correlated with tree age/height as a result of greater belowground carbon allocation of old/taller trees, thereby presenting k sr declining with tree age/height as well. However, new observations are required to verify this hypothesis.
As k sr declines in waterlogged conditions, λ steeply increases, and its values of the two taller trees (Nos 2 and 3, Table 1) are distinct from those of the other four ( Figure 6). This difference, more significant than in k sr , was partially caused by a higher hydraulic resistance between roots and leaves (i.e. lower k rl , Eqs (10) and (12)), represented by their lower J O max (Table 1), E (O) and E (M) (Figure 4). These results are in accordance with the hypothesis on hydraulic resistance being the main limiting factor of trees' assimilation and height growth, namely, the hydraulic limitation hypothesis (HLH). Widely supported by modeling and observations (Hubbard et al. 1999, West et al. 1999, McDowell et al. 2002b, Koch et al. 2004, Martínez-Vilalta et al. 2007b, Drake et al. 2010, Olson et al. 2014, HLH states that increasing hydraulic resistance against water transport from soil to leaves is the key constraint on the transpiration and photosynthesis of old/tall trees Yoder 1997, Ryan et al. 2006). The tests on this hypothesis have always been on very tall trees (Koch et al. 2004) and/or under drought (Sevanto et al. 2014), while the current study suggests that lower transpiration rate and higher MWUE related to lower soilto-leaf conductance are noticeable also in short (10-15 m) trees with excess soil water. However, nutrient cycles might have affected the trees' growth too regarding the slowed decomposition under subarctic conditions. Further research is required to clarify this issue.

Residual and water table temperature
Suggested by the correlation between FM's residual (ε) and WTT (r = 0.228 using the full data and 0.237 in shallow WTD; Figure 7), the model slightly overestimated transpiration rate when WTT was close to 0 • C, but the residual diminished as WTT increased to 5 • C. This correlation may be explained by (i) lower diffusion coefficient of oxygen in colder soils (Chapman-Enskog theory) and (ii) higher plasma membrane resistance in roots to water at low temperature (Lambers et al. 2008), and the main reason behind the latter is likely the reduced activity of aquaporins (Javot and Maurel 2002, Murai-Hatano et al. 2008, Ionenko et al. 2010; but see section Caveats for the potential impacts of seasonality absent from the current structure). Thus, including WTT as an input variable in the model may potentially improve the model's performance on soil-to-root hydraulic conductance. However, on the whole-tree scale, the effects of belowground temperature barely over the freezing point on hydraulic conductance or transpiration are yet unclear (but see Running and Reid (1980), Mellander et al. (2004) and Lintunen et al. (2020) for the influences of wider ranges of soil temperature). Also, the current k − sr θ does not yet explain detailed ecophysiological processes associated with waterlogging or WTT (see the next section for these limitations of the model). Therefore, more observations and process-based studies are necessary before the implementation of WTT.

Caveats
There were several structural or parametric assumptions in the current modeling to overcome the difficulties in measurements and/or computation. Firstly, due to lack of previous observations over a gradient of shallow WTD or high SWC, k − sr θ (Eq. (14)) was modeled using simplified coefficients η m and η p , of which the ecophysiological interpretation is not straightforward, unlike that of ξ m and ξ p for k + sr θ (Eq. (13); Duursma et al. 2008, Nikinmaa et al. 2013. Despite the formula's superficial resemblance to the model of oxygen diffusion in soil (Moldrup et al. 1997), η m and η p lump together a variety of factors, e.g. aquaporin dynamics, respiration-related enzymes dynamics, and soil gas composition. Therefore, this part of the model provides only a heuristic quantification, and more mechanistic understanding is needed for modeling these processes in the future.
Secondly, the parameters not directly related to the focal properties and processes (MWUE, hydraulic conductance and waterlogging) were kept constant through the study. Of the estimated parameters, only ξ m and c were tree specific. The constant or non-tree-specific parameters may have introduced errors into the results due to site-or tree-level variances in ecophysiological properties, e.g. lowest water potential at tree top ( lmin in Eq. (12)), saturation of irradiance reaction (γ in Eq. (4)) and leaf-to-sapwood area ratio (ρ in Eqs (12) and (18)). These errors, in addition to the novel model structures, may have contributed to the ranges of k sr and/or λ results different from those of previous modeling studies (Mäkelä et al. 2004, Duursma et al. 2008, Liu et al. 2020. Particularly, leaf-to-sapwood area ratio in Scots pine has been found decreasing with increasing tree height (Vanninen et al. 1996, McDowell et al. 2002a. Therefore, assuming the ratio constant (ρ = 2500 m 2 leaf m −2 sapwood; Table S1 available as Supplementary data at Tree Physiology Online) should have induced an error in converting sap flow density to leaf-areaspecific transpiration rate (Eq. (18)), albeit the error should be minor due to the limited difference among the heights of the sample trees (Table 1).
Moreover, no processes except k sr were assumed directly correlated with SWC or WTD, while, for instance, the maximum photosynthetic rate (reflected by parameter γ in the current model) may be negatively correlated with SWC under waterlogged conditions (Sojka 1992, Carter et al. 2006. The trees might have acclimated to waterlogging or flooding conditions so that their impacts on hydraulic conductance and MWUE might have been mediated to some extent. Such an acclimation may be similar to the flooding-tolerance observed in Scots pine seedlings grown on peat (Pearson et al. 2013). Therefore, cautious scrutiny on ecophysiological processes and phenomena is necessary when the current modeling method is applied to the trees in different landscapes or soil conditions. A number of measurement points scatter outside the 95% credible interval of predictive error (Figure 4), suggesting considerable measurement error additional to the aforementioned structural or parametric limits. Particularly, the measurement error in WTD and J should be emphasized. In drained peatland, the soil water distribution may be highly correlated with the distance between the points of interest and the nearest draining ditch . Thus, the automatically recorded WTD at the approximate center of the site may have induced error when used as input data supplementary to the tree-specific manual measurement. The error in the thermal dissipation method of measuring sap flow density has been widely discussed, such as on the sensitivity to direct sunshine and the thermal sensitivity of the device (Köstner et al. 1998, Hölttä et al. 2015. These issues in WTD and sap flow density data may have caused the heteroscedasticity in the data model (Eq. (19)), which was ubiquitous across trees but in different Tree Physiology Volume 42, 2022 correlations with E (M) . Considering the limitations in modeling methodology and data, we expect to improve the model with measurements of more environmental factors (e.g. oxygen content and/or diffusion rate in soil) and by more accurate measurement methods.

Conclusion
Our semi-mechanistic model of soil to leaf hydraulic continuum performed well on data collected from six Scots pine trees in a drained peatland with seasonal waterlogging. It performed better than did the original optimal stomata model alone and the other reference model (RM2) structured similarly but without the waterlogging module. Particularly, the full new model's simulated transpiration rate of the sample trees had a much less deviation from the observations than the results of RM2. In the future, the model could be improved on providing more mechanistic expressions of related processes, e.g. oxygen diffusion in soil, roots' water uptake and respiration. Additionally, continuous measurements of more variables, e.g. respiration rate, soil oxygen concentration, and water potentials at canopy and soil, should improve the model's performance and applicability.

Supplementary data
Supplementary data for this article are available at Tree Physiology Online.