Asymptotic and numerical analysis of a porous medium model for transpiration-driven sap flow in trees

TL;DR

This paper presents a comprehensive 3D porous medium model for sap flow in trees, deriving analytical solutions and validating them with numerical simulations and experimental data, to better understand flow regimes and the effects of stem properties.

Contribution

It introduces an asymptotic analysis framework for a nonlinear sap flow model, highlighting the influence of stem aspect ratio and anisotropy on flow behavior.

Findings

Analytical series solutions for saturation and velocity are derived.

Flow regimes are characterized by key non-dimensional parameters.

Model validation shows good agreement with experimental data on Norway spruce.

Abstract

We develop a 3D porous medium model for sap flow within a tree stem, which consists of a nonlinear parabolic partial differential equation with a suitable transpiration source term. Using an asymptotic analysis, we derive approximate series solutions for the liquid saturation and sap velocity for a general class of coefficient functions. Several important non-dimensional parameters are identified that can be used to characterize various flow regimes. We investigate the relative importance of stem aspect ratio versus anisotropy in the sapwood hydraulic conductivity, and how these two effects impact the radial and vertical components of sap velocity. The analytical results are validated by means of a second-order finite volume discretization of the governing equations, and comparisons are drawn to experimental results on Norway spruce trees.

Tables

Table 1. Table 1: Variables and parameters for the 3D sap flow model. Listed values are for the “base case” used in the asymptotic analysis, while parameters from additional simulations are given in parentheses.

Symbol	Description	Units	Value or Formula
Independent and dependent variables:
$r$	radial coordinate	$m$
$z$	vertical coordinate	$m$
$t$	time	$s$
$s(r,z,t)$	sap volume fraction or saturation	–
$S(r,z,t)$	saturation deficit	–
$v_r(r,z,t)$	radial sap velocity	$m/s$
$v_z(r,z,t)$	vertical sap velocity	$m/s$
Solution-dependent functions:
$E\left(t\right)$	transpiration rate per unit leaf area	$m^3/m^{2}s$	Eq. 2.13b
$f\left(z\right)$	combined leaf area and shading effects	–	Eq. 2.13c
$𝗞$	hydraulic conductivity tensor	$m/s$	Eq. 2.3
$K_{r,z}$	hydraulic conductivity components	$m/s$	Eq. 2.11
$\ell \left(z\right)$	leaf area per unit height	$m^2/m$	Eq. 2.13d
$R\left(z\right)$	outer stem radius	$m$	Eq. 2.10
$\lambda \left(z\right)$	sunlight shading effect	–	Eq. 2.13e
$\psi \left(s\right)$	capillary pressure head	$m$	Eq. 2.12
Physical (dimensional) parameters:
$E_o$	transpiration flux amplitude	$m/s$	$1\times {10}^{-9}$ ($3.94\times {10}^{-8}$)
$H$	stem height	$m$	6.7
$K_o$	maximum hydraulic conductivity	$m/s$	$5.36\times {10}^{-7}$
${\ell }_o$	leaf specific area	$m$	15.3
$r_o$	maximum tree radius (at base)	$m$	0.0645
$\tau$	number of seconds per day	$s$	$8.64\times {10}^4$
${\psi }_o$	scaling constant in $\psi \left(s\right)$	$m$	$2.93\times {10}^5$
Dimensionless parameters:
$f_o$	maximum of $f\left(z\right)$		2.6
$n$	capillary pressure exponent		400
$s_o$	maximum saturation		0.574
$\alpha$	exponential stem taper rate		1.42
$\gamma$	heartwood fraction (inner : outer radius ratio)		0 (0–0.75)
$\delta$	relative change from max. saturation		0.01
$\zeta$	stem aspect ratio $=r_o/H$		0.00963
$\eta$	time parameter $=\left(2\pi n{s}_o{H}^2\right)/\left(\tau {\psi }_o{K}_o\right)$		$4.77$
$\kappa$	conductivity ratio $=K_r/K_z$		${10}^{-4}$–${10}^{-2}$
$\mu$	gravity parameter $=nH/{\psi }_o$		0.00915
$\varphi$	transpiration parameter $=2f_o{E}_o\mu /K_o\zeta$		0.00920 (0.363)

Equations

$$\gamma R(z)⩽r⩽R(z)\text{and}0⩽z⩽H.$$

$$\frac{\partial s}{\partial t}+\nabla \cdot \vec{v}=0,$$

$$\vec{v}=-\mathsf{K}\nabla (z+\psi (s)).$$

$$\mathsf{K}=\left[\begin{array}{cc}{K}_r(r,z)& 0\\ 0& K_z(r,z)\end{array}\right].$$

$$\vec{v}=(v_r,v_z)=\left(-K_r\frac{\partial \psi }{\partial r},-K_z\left(1+\frac{\partial \psi }{\partial z}\right)\right).$$

$$\frac{\partial s}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r{K}_r{\psi }^{\prime }\frac{\partial s}{\partial r}\right)+\frac{\partial }{\partial z}\left(K_z{\psi }^{\prime }\frac{\partial s}{\partial z}\right)+\frac{\partial K_z}{\partial z}.$$

$$v_z(r,H,t)=0,$$

$$s(r,0,t)=s_o.$$

$$\vec{v}\cdot \hat{n}{|}_{r=\gamma R(z)}=0,$$

$$\vec{v}\cdot \hat{n}{|}_{r=R(z)}=Q(z,t).$$

$$R(z)=r_o\exp \left(-\frac{\alpha z}{H}\right),$$

$$K_z=K_o{K}^{\ast }(r,z)\text{and}{K}_r=\kappa K_o{K}^{\ast }(r,z),$$

$$\psi (s)={\psi }_o\left[1-{\left(\frac{s_o}{s}\right)}^{1/n}\right],$$

$$Q(z,t)=f(z)E(t),$$

$$\displaystyle E(t)=E_{o}\Re\big{[}1+d_{1}\exp\left({2\pi\mathrm{i}t}/{\tau}\right)+d_{2}\exp\left({4\pi\mathrm{i}t}/{\tau}\right)\big{]}$$

$$f(z)=\frac{\ell (z)\lambda (z)}{2\pi R(z)},$$

$$\ell (z)={\ell }_o{\mathrm{sech}}^2\left(\frac{6z}{H}-2.4\right).$$

$$\lambda (z)=\frac{1}{\pi }\arctan \left(\frac{63z}{H}-50\right)+\mathrm{0.53.}$$

$$\frac{\partial }{\partial r}=\frac{1}{R(z)}\frac{\partial }{\partial \tilde{r}}\text{and}\frac{\partial }{\partial z}=-C\tilde{r}\frac{\partial }{\partial \tilde{r}}+\frac{\partial }{\partial z},$$

$$v_r=-\frac{\kappa D(s)}{R(z)}\frac{\partial s}{\partial \tilde{r}}\text{and}{v}_z=-K_z(s)+\tilde{r}CD(s)\frac{\partial s}{\partial \tilde{r}}-D(s)\frac{\partial s}{\partial z},$$

$$\frac{\partial s}{\partial t}=-\frac{1}{R(z)\tilde{r}}\frac{\partial \left(\tilde{r}{v}_r\right)}{\partial \tilde{r}}+\tilde{r}C\frac{\partial v_z}{\partial \tilde{r}}-\frac{\partial v_z}{\partial z}.$$

$$\displaystyle(\tilde{r}_{i},\,z_{k})=\Big{(}\gamma+{\textstyle\left(i-\frac{1}{2}\right)}\Delta\tilde{r},\;{\textstyle\left(k-\frac{1}{2}\right)}\Delta z\Big{)},$$

$$(v_r{)}_{i+\frac{1}{2},k}$$

$${\left(v_z\right)}_{i,k+\frac{1}{2}}$$

$${\left(v_z\right)}_{i+\frac{1}{2},k}$$

$$\begin{gathered} \frac{\partial s_{i,k}}{\partial t}=-\frac{1}{{\tilde{r}}_i{R}_k\Delta \tilde{r}}\left({\tilde{r}}_{i+\frac{1}{2}}(v_r{)}_{i+\frac{1}{2},k}-{\tilde{r}}_{i-\frac{1}{2}}(v_r{)}_{i-\frac{1}{2},k}\right)+\frac{{\tilde{r}}_{i}C}{\Delta \tilde{r}}\left((v_z{)}_{i+\frac{1}{2},k}-(v_z{)}_{i-\frac{1}{2},k}\right) \\ -\frac{1}{\Delta z}\left((v_z{)}_{i,k+\frac{1}{2}}-(v_z{)}_{i,k-\frac{1}{2}}\right). \end{gathered}$$

$$r=r_o{r}^{\ast },z=H{z}^{\ast },t=\frac{\tau }{2\pi }{t}^{\ast },$$

$$R(z)=r_o{R}^{\ast }(z^{\ast }),f(z)=f_o{f}^{\ast }(z^{\ast }),E(t)=E_o{E}^{\ast }(t^{\ast }),$$

$$s(r,z,t)=s_o\left(1-\delta S(r^{\ast },z^{\ast },t^{\ast })\right),$$

$$\psi$$

$${\psi }^{\prime }$$

$$\left(\frac{2\pi n{s}_o{r}_o^2}{\tau {\psi }_o{K}_o}\right)\frac{\partial S}{\partial t^{\ast }}=\frac{1}{r^{\ast }}\frac{\partial }{\partial r^{\ast }}\left[\kappa K^{\ast }\left(1+{\psi }_1\delta S\right)r^{\ast }\frac{\partial S}{\partial r^{\ast }}\right]+{\left(\frac{r_o}{H}\right)}^2\frac{\partial }{\partial z^{\ast }}\left[K^{\ast }\left(1+{\psi }_1\delta S\right)\frac{\partial S}{\partial z^{\ast }}\right]-\left(\frac{n{r}_o^2}{\delta H{\psi }_o}\right)\frac{\partial K^{\ast }}{\partial z^{\ast }},$$

$$\zeta =\frac{r_o}{H}\approx 0.00963,\eta =\frac{2\pi n{s}_o{H}^2}{\tau {\psi }_o{K}_o}\approx 4.77,\mu =\frac{nH}{{\psi }_o}\approx \mathrm{0.00915.}$$

Full text

{AMS}

35C20, 35K61, 76M12, 76S05, 92C05, 92C80.

1 Introduction

The phenomenon of sap transport in trees has challenged plant physiologists for over a century [11]. The remarkable ability of trees to draw water upward from roots to leaves through heights of 100 metres or more is still not completely understood and still generates a great deal of interest in the biophysics and plant physiology communities [1, 6, 17]. The driving process behind sap flow is transpiration, which can be briefly explained as follows [34]. Water is lost from the leaves as a consequence of the shared path between inward diffusion of carbon dioxide (the essential precursor for photosynthesis) from the ambient air to the leaves, and the combined outward evaporation/diffusion of water due to the vapor pressure deficit. As water evaporates from the leaves, air-water surface tension is generated within the hydrophilic leaf cell interstices and is subsequently transmitted through the sap to conductive wood cells (known as sapwood or xylem) making up the stem and branches. The resultant negative sap pressure (measured relative to atmospheric) generates the driving force necessary to draw water from the roots to the leaves [27]. This transpiration-driven flow is a maximum during the day when there is a high evaporative demand, and drops to a minimum overnight when photosynthesis halts and the internal water storage within the tree is replenished [39].

The physical structure of wood plays an essential role in sap transport. Within coniferous (or softwood) tree species that are the focus of this study the primary conductive elements in sapwood are the tracheids, which are elongated and vertically-oriented dead wood cells having rigid lignified cellulose walls enclosing an empty lumen. Adjacent tracheids are hydraulically connected through paired pits, which are pores that permit flow of liquid sap but are small enough to prevent passage of any air bubbles that might be formed within the tracheids. When the sapwood is viewed as a porous medium, the combination of pit distribution and tracheid orientation engenders a high degree of anisotropy in hydraulic conductivity, which can be several orders of magnitude larger in the vertical direction than the radial [8]. Deciduous (or hardwood) tree species differ from conifers in that sapwood contains an extra class of wood cells called vessels having a much larger diameter and with greater permeability to flow, although the wood structure and transport properties remain otherwise similar.

Because typical pressures within a tree stem exceed the saturated vapor pressure, gases such as air and carbon dioxide are known to exist in a dissolved state within the sap. Under conditions of extreme dryness or freezing, pressure can fall enough to cause dissolved gases to cavitate and form air bubbles inside tracheids or vessels [39]. This process is known as embolism and can lead to blockage of affected conduits that prevents them from conducting sap. Trees have an amazing microstructure that is capable of bypassing and even eliminating such embolisms [10], but the precise mechanisms for embolism formation and recovery are still not completely understood [5]. To avoid having to model the added complexities of embolism formation, we assume in this study that the trees operate under “normal conditions” in which the sapwood remains close to full saturation so that no new embolus forms, and neither does embolism recovery play a significant role.

Numerous mathematical have been developed for studying the flow of sap within conductive sapwood. One of the most popular approaches is an electric circuit analogy, in which the porous flow elements are characterized by a flow resistance (which is inversely proportional to conductivity), combined with capacitive elements that capture storage of sap (and subsequent time lags) within the roots, leaves and various cells making up the stems and branches [14]. The earliest circuit models included only resistive elements and hence failed to capture the dynamic nature of the flow [37], whereas more recent models include water storage effects and hence capture observed lags between transpiration flux and sap flow [26]. Some authors have developed even more detailed models that capture the branching structure of the tree [36], or incorporate the added dynamic effects of radial stem growth [32] by connecting resistor-capacitor elements in a more complex branching structure. A major drawback of these circuit models is that the parameters have no direct correspondence in the context of an actual porous medium flow, and the resistance and capacitance are treated instead as fitting constants that are matched to experiments rather than being directly measured. Furthermore, the circuit analogue is constructed out of a network of resistors and capacitors connected in series, for which the resulting flow depends sensitively on the actual discretization used. Finally, the circuit parameters are usually treated as constants although the transport properties of actual sapwood change in time owing to the local saturation state. Having said that, circuit models are still applied widely in the tree physiology literature because of their simplicity and straightforward algebraic structure.

Another type of model that overcomes many of these deficiencies is the class of porous medium models. The study of unsaturated porous flow in porous media is very well developed [20], especially in the context of groundwater transport in soils [9]. Sap flow is especially suitable for treatment using continuum porous medium models because of the simple repeating microstructure of wood. In this class of models, sap flow is driven by pressure gradients according to Darcy’s law, and the governing equations consist of nonlinear partial differential equations (PDEs) that capture spatial and temporal variations in variables such as water content (saturation) and pressure. Chuang et al. [7] developed a simple 1D porous medium model for transpiration driven flow in a conifer stem, and used numerical simulations to fit their results to experimentally measured sap fluxes. Bohrer et al. [4] extended these results to include the effect of complex branching structure within the crown [7]. Aumann and Ford [2, 3] applied a different approach focused on the wood microstructure in Douglas fir by developing a detailed model for transport within a tracheid network including water-air interface dynamics and many microstructural parameters. Such a model is useful for uncovering detailed aspects of flow within individual wood cells, but it is not very practical for use at the scale of an entire tree. A number of other related PDE models have been also developed [13, 21, 25, 30], but many are one-dimensional and so ignore effects such as radial variations within the stem, not to mention that to date little mathematical analysis has been done to determine the character of the solutions.

The main goal of this paper is to extend the 1D porous medium model from [7] to a more realistic 3D cylindrical (rotationally symmetric) model of a tree trunk, with an imposed transpiration flux distributed along the outer surface. Along with realistic coefficient functions fit to data on Norway spruce, this model will permit study of the radial flow patters that develop with the stem. We will then develop an approximate analytical solution with the aid of asymptotic analysis, identifying the different parameter regimes where anisotropy begins to dominate which provides an alternate explanation for the observed radial variation in the vertical velocity [15, 28] that isn’t due to loss of hydraulic conductivity. We study the relative importance of gravity and transpiration as well as investigating the nature and relative importance of the radial/vertical sap fluxes, which is significant in light of the recent experimental advances that permit measurements of the relatively small radial sap velocity components [12, 35]. We further demonstrate how temporal and spatially localized disturbances in saturation propagate along the tree. Throughout, we use a finite volume discretization of the governing nonlinear PDE in order to solve the problem numerically and validate the asymptotic results.

2 Model for 3D Axisymmetric Flow in a Tapered Cylindrical

Annulus

Our work is based on the one-dimensional sap flow model of Chuang et al. [7], who considered only vertical variations within a tapered cylindrical tree stem. We extend their model to a 3D axisymmetric stem geometry by incorporating the effect of radial variations as well as a core heartwood region that is impermeable to flow. Except for slight changes in notation and the need for several additional boundary conditions, the model is much the same as the 1D analogue. A description of all parameters and solution-dependent coefficient functions is provided in Section 2.2, which are obtained from the literature or fit to experimental measurements from Norway spruce (Picea abies) provided in [7].

2.1 Governing Equations and Boundary Conditions

The Norway spruce is a conifer species whose stem is well-approximated by a right circular cylinder that tapers from base to crown. For a tree of height $H$ (in units of $\mathrm{m}$) the outer stem radius may be written as $r=R(z)$ $[\mathrm{m}]$, where the axial coordinate $z$ satisfies $0\leqslant z\leqslant H$ and $R(z)$ is a decreasing function of $z$ (which will be specified later in Section 2.2). We extend the 1D model from [7] by taking a more realistic geometry pictured in Fig. 1a where the stem consists of an outer layer of conducting sapwood that surrounds a core region of non-conducting heartwood. The heartwood is assumed to take up some fraction $0<\gamma<1$ of the stem cross-section corresponding to $0\leqslant r\leqslant\gamma R(z)$. We are therefore concerned with capturing sap transport within the outer conductive portion which has the shape of a tapered cylindrical annulus defined by

[TABLE]

In the absence of any directional forcing around the stem, we can suppose rotational symmetry and neglect any dependence on the polar angle.

The sapwood is treated as a variably-saturated medium whose porous structure contains a mixture of two phases: liquid (sap) and gas (mostly air). The saturation or local pore volume fraction containing liquid is denoted by $s(r,z,t)$, and depends on location and time $t$ $[\mathrm{s}]$. Enforcing conservation of liquid yields the continuity equation

[TABLE]

where the sap velocity $\vec{v}(r,z,t)$ $[\mathrm{m/s}]$ obeys Darcy’s law

[TABLE]

Here $\psi(s)$ $[\mathrm{m}]$ is the hydrostatic pressure head, which is a saturation-dependent function that is specified later in Section 2.2. The hydraulic conductivity tensor $\boldsymbol{\mathsf{K}}$ $[\mathrm{m/s}]$ depends in general on both the location within the tree stem (due to spatial variations in pore structure) and the local saturation state. The porous structure of wood is highly anisotropic since it is composed of dead, hollow, elongated cells that are directed vertically within the xylem; in conifers, these cells are tracheids, whereas in deciduous trees the primary route for sap transport is through much larger cells called vessels. Microscopic “pits” or pores connect adjacent wood cells in the radial direction, as well as a much smaller number of radially-directed “ray cells”; however, these structures are much less permeable to flow than the vertically-oriented tracheids, so it is reasonable to assume that the anisotropic hydraulic conductivity has the form of a diagonal tensor

[TABLE]

Note that $\boldsymbol{\mathsf{K}}$ depends on position but not on saturation, which we justify based on the earlier assumptions that sapwood remains nearly saturated and no significant build-up of gas occurs owing to absence of embolism formation. This should be contrasted with Chuang et al. [7] who took the conductivity to be a function of saturation; however, since we are primarily concerned with trees near the fully saturated state where the hydraulic conductivity is relatively insensitive to changes in $s$, taking $\boldsymbol{\mathsf{K}}$ to be independent of $s$ is a reasonable approximation (to be discussed more detail in the next section).

In most trees, the dominant vertical orientation of conductive cells induces anisotropy in $\boldsymbol{\mathsf{K}}$ that ranges from moderate to extreme. Therefore, the analytical derivations in this paper focus on two limiting cases: an isotropic conductivity with $K_{r}\sim K_{z}$, for which a general series solution is relatively straightforward to derive; and high anisotropy with $K_{r}\ll K_{z}$, which can be treated analytically under certain restrictions. In either case, numerical simulations are used to study the solution over the entire range of anisotropy.

Using the above definitions, the velocity can be expressed in cylindrical coordinates as

[TABLE]

After substituting these two expressions into the continuity equation Eq. 2.2, we obtain the following nonlinear parabolic PDE for saturation

[TABLE]

Moving on to the boundary conditions, no sap flows through the top of the stem so that

[TABLE]

whereas the base of the stem (at the roots) is assumed to be fully saturated with

[TABLE]

On the inner sapwood boundary, a zero flux is imposed

[TABLE]

where $\hat{n}(z)=(1,-R^{\prime})/\sqrt{1+(R^{\prime})^{2}}$ is the unit outward-pointing surface normal. This last condition reflects the fact that no sap exchange occurs with the non-conducting heartwood (and reduces to a simple radial symmetry condition in the case $\gamma=0$).

The final boundary condition on the outer stem surface derives from the transpiration flux, which is the process whereby water is drawn upward from the roots, through the stem and branches to the leaves (or needles). Transpiration is driven by leaf evaporation and the subsequent vapor pressure deficit between the leaf interior and the surrounding atmosphere. A particular feature of species like Norway spruce is that branches are distributed densely along the entire stem so that the transpiration flux can be specified as a corresponding distribution in the axial direction (this should be contrasted with many deciduous species in which branches are concentrated within the crown at the top of a long, bare trunk). Therefore, we impose a sap outflow due to branches distributed continuously along the stem by means of another flux boundary condition

[TABLE]

The transpiration flux $Q$ $[\mathrm{m^{3}/m^{2}\,s}]$ depends on time owing to diurnal variations in transpiration, as well depending on the branch distribution along the stem, and we provide a specific functional form for $Q(z,t)$ in the next section. A primary difference from the 1D model in [7] is that the transpiration flux was incorporated as a source term in the saturation equation Eq. 2.5 rather than as a boundary condition. A graphical summary of the geometry and boundary conditions is provided in Fig. 1b.

2.2 Parameters and Coefficient Functions

We now provide parameter estimates and functional forms for the model coefficients, which are chosen to match as closely as possible the data for Norway spruce provided in [7]. All symbols are listed in Table 1 along with units and values. We emphasize that neither the model nor the asymptotic analysis in the following sections is restricted to a particular species of tree, and that the specific functional forms introduced here are not critical for the analysis. Rather, we exploit these parameters as a convenient illustration of our 3D model, which has the added benefit of allowing us to draw a concrete comparison with the results in [7].

We first specify the form of the outer tree stem radius that tapers with height according to

[TABLE]

where $\alpha$ controls the rate of taper from roots to crown. This choice of exponential function was motivated in [7] for reasons of mathematical convenience. However, there is an extensive literature on more complicated diameter-versus-height relationships [22, 24], from which we observe that many coniferous and deciduous tree species have a small enough taper rate that such an exponential function provides a reasonable approximation of stem shape. However, we still develop the majority of our analytical results for the general function $R(z)$ and impose Eq. 2.10 when we need to exploit additional simplifications in Section 4.6, and in numerical simulations when a specific form of $R(z)$ is required.

Another important geometric parameter is the heartwood ratio $\gamma$ that determines the thickness of the annular sapwood region. For young trees that have not yet developed a well-defined heartwood region, taking $\gamma=0$ is a reasonable choice. For spruce trees of the age and height considered in [7] the typical heartwood fraction is $\gamma\approx 0.5$ [31], while for other species $\gamma$ can be as high as 0.75 [38]. For most of this paper including the asymptotic developments we assume that $\gamma=0$, although we explain later how our results can be extended to the case when the heartwood fraction is much larger.

Based on our earlier assumption that the hydraulic conductivity is a diagonal tensor with entries that depend only on position, we take the axial and radial conductivities to have the form

[TABLE]

where $K_{o}$ $[\mathrm{m/s}]$ is a constant equal to the maximum value of axial hydraulic conductivity, and $K^{*}$ is a dimensionless function that is strictly positive. The dimensionless factor $\kappa$ is the ratio of radial to axial conductivity that captures the degree of anisotropy in the sapwood. Typical values of $\kappa$ lie between $10^{-4}$ to $10^{-2}$ [8, 29], although we will also consider the case when $\kappa\sim{1}$ (where we use $\sim$ to denote asymptotic equivalence in which the two quantities have the same order of magnitude).

Many models of variably-saturated porous media specify the hydraulic conductivity as a function of saturation; indeed, Chuang et al. [7] imposed a Weibull-type function of the form $K^{*}(s)=\exp(-a|\psi(s)|^{b})$. However, their simulations remained within 20% of the fully saturated state $s\approx s_{o}$, for which $K^{*}(s)$ has a nearly linear dependence on $s$ with a very small (negative) slope. The reason for this behaviour is that the primary cause of hydraulic conductivity loss in trees is embolism formation, which only occurs when saturation drops sufficiently below $s_{o}$. Since we expect that variations in saturation remain relatively small, it is reasonable to approximate variations in conductivity using a simpler spatially dependent function $K^{*}(r,z)$ for which we derive most of our asymptotic results (and which avoids a nonlinear $K^{*}$). After investigating the asymptotic solution for spatially-dependent conductivity, we then consider in Section 4.6 the special case of constant conductivity and a periodic transpiration rate, which allows us to derive a simpler closed–form solution for the leading order term.

The hydraulic pressure head is taken to depend on saturation according to

[TABLE]

where $\psi_{o}$ and $n$ are fitting parameters. This is similar to the van Genuchten model commonly used for capillary pressure in soil and rock [33], and also applied to drying of lumber [18]. Note that the head is a negative quantity because $s<s_{o}$, reflecting the understanding in the sap hydraulics literature that sap within a tree stem is under tension. Our asymptotic derivation is not specific to this or any other particular form of the capillary pressure function, but does rely on two essential features: namely, that $\psi(s)$ is a smooth and monotone increasing function in a neighbourhood of $s=s_{o}$.

The final ingredient in the model specification is the transpiration source term, which we assume takes the separable form

[TABLE]

The height-dependent transpiration factor is decomposed as

[TABLE]

where $2\pi R(z)$ is a geometric scaling factor (equal to stem perimeter) and $\ell(z)$ is leaf area density (units of $\mathrm{m^{2}/m}$) that captures the impact of sun exposure on transpiration and is given in [7, Fig. 4] as

[TABLE]

The remaining factor $\lambda(z)$ is a dimensionless quantity called transpiration flux density that captures shading effects due to branches/leaves located above a given height. Although $\lambda(z)$ was not provided in [7], it can be approximated using experimental data in [7, Fig. 5] as

[TABLE]

To summarize, the 3D sap flow model corresponds to solving Eqs. 2.5 to 2.9 along with a suitable initial condition on saturation and the function definitions in Eq. 2.10–Eq. 2.13.

3 Numerical Method

We next describe a numerical method for solving the governing equations based on a cell-centered finite volume approximation. This will be used to validate our asymptotic results in the case of an isotropic conductivity, and to produce comparisons for the anisotropic case. To simplify the discrete equations, it is helpful to first transform the radial coordinate for the tapered annular cylindrical domain.

3.1 Coordinate Transformation

The radius of the tapered cylinder obeys $\gamma R(z)\leqslant r\leqslant R(z)$, which suggests defining a transformed radial coordinate $\tilde{r}=r/R(z)$ that is bounded between $\gamma\leqslant\tilde{r}\leqslant 1$ for all $0\leqslant z\leqslant H$. The spatial derivatives within the governing equations can then be transformed via

[TABLE]

where $C:=R^{\prime}/R=-\alpha/H$ is a constant owing to the special exponential form Eq. 2.10 assumed for $R(z)$. Applying these transformations to the velocity components Eq. 2.4 yields

[TABLE]

where $D(s):=K_{z}(s)\psi^{\prime}(s)$. Here the velocities are written in terms of a general saturation dependent hydraulic conductivity function $K_{z}(s)$, which is useful for later comparison to the results in [7] (see Section 3.3); however, the same formulas extend easily to a spatially-dependent conductivity $K_{z}(r,z)$. Making use of these velocity components, the continuity equation Eq. 2.5 becomes

[TABLE]

The primary advantage to transforming the radial coordinate in this manner is that the equations above are now imposed on a rectangular computational domain in $(\tilde{r},z)$ space, for which standard finite difference approximations can be applied.

3.2 Cell-centered Finite Volume Discretization

We now discretize the transformed governing equations using a cell-centered finite volume scheme, so that we preserve as accurately as possible the conservation of mass embodied in the continuity equation. The computational domain is divided into an $N_{r}\times N_{z}$ rectangular grid of cells having centers

[TABLE]

with $\Delta\tilde{r}=(1-\gamma)/N_{r}$ and $\Delta z=H/N_{z}$ for $i=1,2,\dots,N_{r}$ and $k=1,2,\dots,N_{z}$. As shown in Fig. 2, the saturation $s_{i,k}$ is approximated at cell centers, whereas the velocity components are located at the center of each cell edge, denoted for example by $(v_{r,z})_{i\pm\frac{1}{2},k}$ and $(v_{r,z})_{i,k\pm\frac{1}{2}}$. The spatial derivatives in Eq. 3.1 are then replaced using compact centered difference formulas to obtain edge-centered velocities

[TABLE]

Note that to maintain both stencil compactness and second-order accuracy we have had to introduce approximate values of saturation at cell edges and corners along with corresponding approximations for $D(s)$, all of which are computed using appropriate arithmetic averages of cell-centered saturation values. The velocity components are then substituted into the discrete continuity equation to get

[TABLE]

To discretize the boundary conditions, we set the normal velocity at the center of each boundary cell edge to the velocity specified in the corresponding boundary conditions Eq. 2.6, Eq. 2.8 and Eq. 2.9. As for the remaining lower boundary, we introduce a band of fictitious cells with centers located one-half grid spacing below the boundary and define values of saturation at these fictitious points. Using the boundary cell values and the lower boundary edge value in Eq. 2.7, a linear extrapolation is used to calculate the saturation in the fictitious cells. Complete details of the spatial discretization are provided in [16]. To integrate the equations in time, we use a method-of-lines approach in which the spatially-discrete equations are treated as a system of time-dependent ODEs, which is then integrated in time using Matlab’s stiff ODE solver ode15s with tolerance values $\mathtt{ABSTOL}=\mathtt{RELTOL}=1\times 10^{-10}$.

In order to verify that our Matlab implementation yields the desired order of spatial accuracy, we performed a numerical convergence study by choosing a sequence of spatial grids with $N_{r}=N_{z}=32,64,128,256,512$ and computing the solution on a fixed time interval that roughly reaches a steady state. The error for each simulation is estimated by treating the fine-grid solution as the “exact solution” and then calculating the 1-norm difference. The resulting errors exhibit a convergence rate of approximately 1.99 which is strongly indicative of second-order accuracy.

3.3 Comparison With Experimental Data on Norway Spruce

We now compare results from the 3D numerical scheme with the experimental data on vertical sap flux provided by Chuang et al. [7]. For this simulation, we use parameters listed in Table 1, except that in order to be consistent with the 1D model we assume the entire stem is made of conductive sapwood ($\gamma=0$) that is isotropic ($\kappa=1$). Figure 3a depicts the measured sap flux (in the vertical direction) alongside our finite volume simulations at two different times (noon and 4:00 pm). The model results clearly capture the overall solution behaviour, exhibiting the same “double-peak” behaviour in sap flux that can be attributed to the bimodal behaviour of the transpiration function $f(z)$ from Eqs. 2.13c to 2.13e (shown in Fig. 3b). In addition to this qualitative agreement, we note that the magnitude of the sap flux is also reasonably well-approximated by computations. Simulations of the 1D model equations of Chuang et al. were also performed using an analogous finite volume scheme and the results were indistinguishable to the eye from the results in Fig. 3.

The primary reasons for developing this 3D sap flow model are to capture radial velocity (as well as radial variations in the solution) and to investigate the impact of including a non-conducting heartwood region with $\gamma>0$. With this in mind, we performed a series of three simulations with different heartwood thickness ($\gamma=0$, 0.5 and 0.75) and plotted the resulting radial velocities in Fig. 4. Two positive peaks appear in $v_{r}$ which clearly derive from the local maxima in the transpiration flux, and these are offset by a comparatively large negative radial velocity at the tree base due to root influx. This is a geometric effect that mimics the inward radial tilt of sapwood vessels, which due to stem taper is largest at the base. The effect of this radially-inward flow is accentuated as the thickness of the annulus decreases (i.e., as $\gamma$ increases) in order to maintain a total mass balance that matches the specified outward transpiration flux.

4 Asymptotic Analysis

The asymptotic derivation in this section bears some resemblance to Kevorkian and Cole’s analysis of heat conduction in a long circular rod [19, chapter 4]. Our results hold for the general case of sap flow in an annular stem cross-section Eq. 2.1 as well as a saturation-dependent conductivity, and we make no assumption on the functional form of either $R(z)$ or the transpiration flux defined via $f(z)$ and $E(t)$. However, for the time-varying solution we will exploit simplifications that arise from assuming an exponential taper Eq. 2.10 and constant hydraulic conductivity. We begin by considering the special case when the radial and axial conductivities $K_{r}$ and $K_{z}$ are the same order of magnitude so that $\kappa\sim{1}$, and defer study of the anisotropic case $K_{r}\ll K_{z}$ to Section 4.5.

4.1 Non-Dimensionalization in the Nearly-Saturated Regime

We next recast the 3D model in dimensionless form in order to reduce the number of free parameters and identify important dimensionless groupings (an analogous version of the 1D model can be derived as well, be we don’t present it here). Introduce the following non-dimensional variables

[TABLE]

where $\delta$ measures the relative deviation from full saturation and typically satisfies $0<\delta\ll 1$. In practice, $\delta$ increases with transpiration rate $E_{o}$, although its precise value is not so important since $\delta$ ultimately cancels out from the final result in our asymptotic derivation. A rough estimate for $\delta$ can be found by recognizing that the threshold pressure for embolism in spruce has been measured as \psi\approx-255\;\mbox{\mathrm{m}} [23], corresponding to a saturation of $s\approx 0.405$. Because we are interested in “normal” flow conditions that are well-removed from any embolized state, we should thus constrain $\delta\lesssim 0.1$. In particular, using a “base value” of $E_{o}=1\times 10^{-9}$, simulations yield a $\delta=\mathcal{O}\!\left(10^{-2}\right)$ so that saturation is guaranteed to remain close enough to $s_{o}$ (further discussion of $\delta$ is provided at the end of this section).

This rescaling has been performed such that $S$ and all “starred” variables can be considered $\mathcal{O}\!\left(1\right)$ quantities. Before rescaling the governing equations, we note that the hydrostatic pressure derivative $\psi^{\prime}(s)$ appears in the saturation equation Eq. 2.5 and the velocity boundary conditions through the Darcy velocities Eq. 2.4. Because $\psi$ depends nonlinearly on saturation in Eq. 2.12, we will ultimately need to expand $\psi$ as a power series as well. Therefore we perform this expansion step for $\psi$ now, which will also permit us to determine an appropriate set of dimensionless parameters that characterize the problem. To this end, we change variables using Eq. 4.1 and expand $\psi$, $\psi^{\prime}$ as power series in $\delta$ to obtain

[TABLE]

where $\psi_{1}=1+\frac{1}{n}$ and $\psi_{2}=\left(1+\frac{1}{n}\right)\left(2+\frac{1}{n}\right)=2+\mathcal{O}\!\left(\frac{1}{n}\right)$ (recalling that $n=400$ is large). Next, substitute the rescaled variables Eq. 4.1 into the saturation equation Eq. 2.5, which after rearranging yields

[TABLE]

noting that only the first two terms in the series Eq. 4.2 are required. Upon careful consideration of the various factors multiplying each term, we are led to introduce the following three dimensionless ratios

[TABLE]

The parameter $\zeta\ll 1$ has an obvious physical interpretation as the stem aspect ratio and will play a central role as the primary expansion parameter in our asymptotic analysis. The rescaled saturation equation may be rewritten in terms of these parameters (after dividing by $\zeta^{2}$) as

[TABLE]

where all “stars” have been omitted on dimensionless quantities to streamline notation. For transpiration-driven flow we are interested only in the case where gravitational effects are significant at leading order, which means that the final term in Eq. 4.4 must have $\mu/\delta\sim{1}$. This implies that all terms are balanced except possibly the second, whose relative importance depends on the anisotropy parameter through the ratio $\kappa/\zeta^{2}$ – in practice, $\kappa$ varies between $1$ and $\zeta^{-1}$ although we will begin considering the lower limit $\kappa\sim{1}$ in the asymptotics.

We next consider the boundary conditions and assume for the present that $\gamma=0$ in Eq. 2.1, which corresponds to a stem consisting entirely of conducting sapwood (i.e., no heartwood). This does not restrict the generality of our results since we can still capture the effect of $\gamma>0$ by introducing an appropriate spatial variation in $K$. Taking $\gamma=0$ introduces a major simplification in that the non-dimensional form of Eq. 2.8 becomes

[TABLE]

The corresponding top and bottom boundary conditions Eqs. 2.6 and 2.7 become respectively

[TABLE]

It is in this second boundary condition that the earlier assumption of $\mu\sim\delta$ ensures that the gravitational effects also contribute to the top boundary condition at leading order.

The derivation of the transpiration flux condition Eq. 2.9 at the outer trunk surface is complicated by the presence of the curved boundary where the normal direction is not aligned with a coordinate axis. After changing variables and expanding terms involving $\psi^{\prime}$, this boundary condition reduces to

[TABLE]

where we have introduced the new dimensionless parameter

[TABLE]

which means that $\phi\sim\zeta$. Thus, based on the “base value” for $E_{o}=1\times 10^{-9}$, the primary dimensionless parameters obey the equivalence $\delta\sim\zeta\sim\mu\sim\phi$. We note nonetheless that the asymptotic solution still yields an accurate approximation for much larger values of $E_{o}$ that violate this equivalence, which we will see later in Fig. 9.

Because of the small aspect ratio, the radial and vertical components of the outward-pointing normal may be expanded as power series in small $\zeta$:

[TABLE]

These expressions are exploited in the next section to simplify the boundary condition Eq. 4.7.

4.2 Asymptotic Expansion: General Case

Based on the assumption that the saturation in Eq. 4.1c stays close to its maximum value $s_{o}$, we seek a regular power series expansion of $S$ in terms of the small parameter $\zeta$ as

[TABLE]

Our aim is to derive equations for the first two terms $S_{0}$ and $S_{1}$ so as to capture the effect of the nonlinearity in pressure head Eq. 4.2a. We will also require two additional terms up to $\mathcal{O}\!\left(\zeta^{3}\right)$ (involving $S_{2}$ and $S_{3}$), which are needed for matching purposes to obtain a closed set of equations for $S_{0}$ and $S_{1}$. Furthermore, we will see shortly that the two leading order solutions are independent of $r$ so that $S_{2}$ is also required to determine the leading order term in the radial sap velocity. Taking the above expansion for saturation, the leading order equation from Eq. 4.4 at $\mathcal{O}\!\left(\zeta^{-2}\right)$ is simply

[TABLE]

Because conductivity $K$ must be nonzero, it follows that $S_{0}$ is independent of the radial coordinate and so $S_{0}=S_{0}(z,t)$. The same equation governs $S_{1}$ at $\mathcal{O}\!\left(\zeta^{-1}\right)$, so we conclude likewise that $S_{1}=S_{1}(z,t)$. Considering the boundary condition Eq. 4.7 at $r=R(z)$, if transpiration is to have any effect on saturation at the first two orders $S_{0}$ and $S_{1}$ (where the radial derivative term vanishes) then the remaining terms should balance. This is equivalent to requiring $\phi\sim\delta$, which is satisfied for the parameters in Table 1 with $E_{o}=1\times 10^{-9}$ as discussed at the end of Section 4.1.

The next order equation in Eq. 4.4 involving terms at $\mathcal{O}\!\left(\zeta^{0}\right)$ is

[TABLE]

which can be integrated in $r$ and then evaluated at $r=R(z)$ to obtain

[TABLE]

Taking the same order terms arising in the transpiration boundary condition Eq. 4.7 yields

[TABLE]

where the last term is included at this order because $\frac{\phi}{\delta}\sim 1$. Eliminating $S_{2}$ from these last two equations leads to

[TABLE]

which simplifies to

[TABLE]

where we have defined

[TABLE]

Equation Eq. 4.12 can then be solved subject to the leading order boundary conditions from Eq. 4.6:

[TABLE]

Next we consider the $\mathcal{O}\!\left(\zeta\right)$ terms in Eq. 4.4 which lead to the equation

[TABLE]

where we used the fact that $\delta\sim\zeta$. Integrating from [math] to $r$ and evaluating at $r=R(z)$ yields

[TABLE]

and a second expression for $S_{3}$ comes from the boundary condition at the corresponding order in Eq. 4.7

[TABLE]

These last two equations can be combined to eliminate $S_{3}$, and then further simplified by applying Eq. 4.11 to obtain

[TABLE]

Finally, using the fact that $S_{0}$ and $S_{1}$ are independent of $r$, the integral terms can be evaluated and the governing PDE written in the more compact form

[TABLE]

This PDE can be solved for $S_{1}$ provided that we impose appropriate boundary conditions, which are obtained from the $\mathcal{O}\!\left(\zeta\right)$ terms in Eq. 4.6 as

[TABLE]

Note that the two terms in the top boundary condition are in balance only if $\mu\sim\zeta$, which is satisfied to a very good approximation.

4.3 Steady State Solution (Constant Transpiration Rate)

Explicit solutions can be derived for the two leading order asymptotic terms in the special case where the transpiration rate is constant ($E(t)\equiv 1$) and the saturation has reached a steady state. In this case, the $S_{0}$ equation Eq. 4.12 can be integrated twice using the boundary conditions Eq. 4.14 to obtain

[TABLE]

The first term in $\bar{S}_{0}$ represents the effect of gravity due to changes in pressure with height. The second term captures the net effect of transpiration, where the expression in square brackets denotes total transpiration flux due to branches located above height $z^{\prime}$, evaluated per unit stem conductivity. The next order saturation correction is obtained by integrating Eqs. 4.15 to 4.16 in a similar manner, yielding

[TABLE]

It is insightful at this point to draw an analogy between the formula for the leading order saturation $\bar{S}_{0}$ and electric circuit representations commonly used to model tree sap hydraulics. To this end, we neglect the effects of gravity in Eq. 4.17 and rewrite the remaining integral term as

[TABLE]

can be interpreted as an average resistance to flow within the cross-sectional stem slice lying between $z^{\prime}$ and $z^{\prime}+dz^{\prime}$. This expression has the form of a harmonic average of conductivities $K$, which is analogous to the formula relating electrical resistances and conductances arranged in a series circuit. The corresponding “current” is

[TABLE]

which represents the total transpiration rate drawn through branches between $z=z^{\prime}$ and $z=1$. The product $\mathcal{R}\cdot\mathcal{J}$ is integrated in Eq. 4.19 along the height of the tree in order to obtain the sap potential (or “voltage”), which when properly scaled yields the local saturation state.

4.4 Steady State Solution With Constant Conductivity

Our aim in this section is to estimate the ratio of radial to axial sap velocity and show that the reduction in radial-versus-axial flow is primarily due to the small stem aspect ratio ($\zeta$) and not the conductivity ratio ($\kappa$) as one might expect. We have already exploited the fact that variations in hydraulic conductivity with saturation are small by assuming that $K(r,z)$ varies only with location. We now assume further that spatial variations in $K$ are likewise small so that it is reasonable to take the conductivity function $K\equiv 1$.

Based on this assumption, Eq. 4.10 may be integrated to obtain $S_{2}(r,z)$, which is the first term contributing to the radial velocity and also the lowest-order term containing any radial dependence. The resulting equation at steady state obeys

[TABLE]

which may be substituted into Eq. 2.4 along with earlier approximations and (4.17) to obtain a leading order expression for the (dimensional) radial velocity

[TABLE]

Note that the conductivity ratio is absent from this expression for $\bar{v}_{r}$, which explains our earlier remark that the anisotropy only influences the radial velocity via higher order terms in the asymptotics.

Moving on to the axial velocity, we make use of the fact that when $K$ is constant equation Eq. 4.13 reduces to $G(z)=\frac{1}{2}R^{2}(z)$, which leads to a simpler form of $\bar{S}_{0}$ in Eq. 4.17 that gives the leader order saturation derivative as

[TABLE]

This expression can then be substituted into the axial component of Eq. 2.4 and simplified to obtain the estimate

[TABLE]

The relative magnitude of the two velocity components may then be approximated by

[TABLE]

after dropping any $\mathcal{O}\!\left(1\right)$ terms. Consequently, the velocity ratio at steady state depends to leading order solely on stem aspect ratio, and furthermore the radial velocity is a factor of roughly 100 times smaller than the vertical component. The effect of material anisotropy on sap flow (through small $\kappa$) has no impact at leading order, instead entering only via higher order terms in the asymptotics. This is the main reason that in earlier sections we restricted the asymptotic analysis to the case $\kappa=\mathcal{O}\!\left(1\right)$, since taking $\kappa\ll 1$ only pushes the effects of the anisotropy to higher order without impacting the leading order solution.

To illustrate these asymptotic results, we present in Fig. 5a,b plots of the vertical velocity and saturation, determined using the two leading order terms in the steady state asymptotic solution. In both cases, the solution variables are averaged across the stem cross-section, and the corresponding finite-volume numerical solution is included in order to demonstrate the close correspondence. The vertical velocity exhibits the characteristic double peak that was observed for the variable-transpiration problem in Fig. 3. To illustrate the relative magnitude of the velocity components, Fig. 5c depicts the log of the ratio $v_{r}/v_{z}$. Except for a thin boundary layer adjacent to the top boundary the radial component is at least a factor of 100 smaller than the vertical component, which is consistent with our asymptotic estimate of $\bar{v}_{r}/\bar{v}_{z}$ in Eq. 4.22. The final plot in Fig. 5d provides a clearer picture of the actual flow direction within the stem by depicting both streamlines and direction field arrows.

4.5 Anisotropic Conductivity With $\boldsymbol{\kappa\ll 1}$

Our earlier assumption that the hydraulic conductivity ratio is $\kappa=K_{r}/K_{z}=\mathcal{O}\!\left(1\right)$ led to an asymptotic result in which the two leading order saturations $S_{0}$ and $S_{1}$ were independent of the degree of anisotropy in $K$. However, sapwood is in reality a highly anisotropic material with $10^{-4}\lesssim\kappa\lesssim 10^{-2}$ [8, 29], which leads us to ask how taking values of $\kappa\ll 1$ might alter the asymptotic solution behaviour. In particular, the physical values of anisotropy and expansion parameter $\zeta=\mathcal{O}\!\left(10^{-2}\right)$ suggest considering two asymptotic limits: $\kappa=\mathcal{O}\!\left(\zeta\right)$ and $\kappa=\mathcal{O}\!\left(\zeta^{2}\right)$. For reasons of simplicity, we will continue assuming that the conductivities are constant and the solution is at steady state.

First, consider the case of moderate anisotropy with $\kappa=\zeta$, which corresponds to taking $K_{z}=K_{o}$ and $K_{r}=\zeta K_{o}$ (which we note is only an order of magnitude estimate and doesn’t presume any explicit dependence of $\kappa$ on $\zeta$). The saturation equation Eq. 4.4 then becomes

[TABLE]

and expanding $S$ in the form of a power series Eq. 4.9 yields the leading order solution $S_{0}=\bar{S}_{0}(z)$, which is identical to that obtained in the previous section. However, the next order equation for $S_{1}$ contains an additional term that introduces a radial dependence of the form

[TABLE]

where the function $h(z)$ is determined similarly as in the previous section (so we omit the details here).

In the second case of a more extreme anisotropy with $\kappa=\zeta^{2}$, the saturation equation Eq. 4.4 becomes

[TABLE]

so that the radial and vertical dependence are now fully coupled at all orders. Assuming for simplicity that the stem has no taper (i.e., $\alpha=0$ and $R(z)\equiv 1$) we may transform the leading order solution using $\hat{S}_{0}(r,z)=S_{0}(r,z)-\frac{\mu z}{\delta}$, which yields Laplace’s equation $\Delta\hat{S}_{0}=0$ in cylindrical coordinates, along with boundary conditions

[TABLE]

The advantage of transforming $\hat{S}_{0}$ in this manner is that the $z$ boundary conditions become homogeneous, and hence separation of variables may be applied to obtain the series solution

[TABLE]

where $\lambda_{n}=\pi(n+1/2)$, $I_{0}$ is the zero’th order modified Bessel function of the first kind and

[TABLE]

The effect of anisotropy on the solution is investigated in Fig. 6 where we compare the simulated vertical velocity profiles for $\kappa=\zeta^{p}$, using the three exponents $p=0,1,2$ and taking 6 terms in the Fourier–Bessel series for $p=2$. These results are computed assuming a tree with no taper ($\alpha=0$) and constant transpiration rate. For each $\kappa$, we plot $v_{z}$ as a function of radius at five heights corresponding to the points labelled A–E in Fig. 3b. For the isotropic or moderately anisotropic cases ($p=0,1$) the velocity remains essentially constant with radius, whereas the extreme case of $p=2$ exhibits significant radial variations. This is consistent with our asymptotic results which show that radial dependence only enters the leading order solution when $\kappa=\mathcal{O}\!\left(\zeta^{2}\right)$, and may help to explain the radial dependence in velocity that was observed experimentally in [15, 28].

This sequence of simulations was then repeated for a tapered stem with $\alpha=1.42$ and the corresponding velocity plots are shown in Fig. 7. We observe similar behaviour to the previous cases except that the $\kappa=\zeta^{2}$ results have a more pronounced radial variation. Even for the moderately anisotropic case ($\kappa=\zeta$), there is a slight radial dependence visible in the bottom-most $v_{z}$ plot (location A). It is also interesting to note that introducing stem taper causes a significant drop in vertical velocity near the tree base owing to the increase in sapwood cross-section there; this should be comtrasted with the untapered case where the vertical velocity increases monotonically with height.

We performed one further validation of in the extreme case $\kappa=\zeta^{2}$ by comparing a numerical simulation with the 6-term series solution Section 4.5. The relative difference in the saturation deficit $s_{o}-s=\delta S$ between the asymptotic and numerical solutions is of order $\mathcal{O}\!\left(10^{-2}\right)$, which is the same order as the missing correction term in the asymptotic expansion of $S$, as expected.

4.6 Time-Dependent Transpiration Source

Next, consider a more general time-dependent source where the transpiration rate $E(t)$ is expanded as a Fourier series

[TABLE]

where $d_{0}=1$ and the remaining $d_{m}$ are the complex Fourier coefficients. In order to obtain a closed-form solution, we continue to exploit simplifications that arise from assuming a constant, isotropic conductivity ($K_{r}=K_{z}=1$ in dimensionless variables). We also make explicit use of the exponential form $R(z)=\exp(-\alpha z)$ for the stem taper function, so that equation Eq. 4.13 yields simplified formulas for both $G(z)=\frac{1}{2}R(z)^{2}$ and $G^{\prime}(z)=-\alpha R(z)^{2}$. The leading order saturation $S_{0}$ is then taken to have an analogous Fourier mode decomposition

[TABLE]

where $S_{0}^{m}(z)$ are unknown functions and the transpiration coefficients $d_{m}$ are introduced as scaling factors to simplify later expressions. Substituting these two series into the leading order equation Eq. 4.12 and collecting terms in $\exp({\mathrm{i}mt})$ yields a sequence of linear ODEs for $S_{0}^{m}(z)$

[TABLE]

where

[TABLE]

The corresponding boundary conditions for $S_{0}^{m}$ are obtained from Eq. 4.14 as

[TABLE]

The leading order term $S_{0}^{0}(z)$ has already been determined as the steady solution $\bar{S}_{0}(z)$ in Eq. 4.17. For the remaining ODEs with $m\geqslant 1$, we split the general solution of Eq. 4.26 into the sum of homogeneous and particular solutions as

[TABLE]

Proceeding to the next order in the asymptotic solution for $S_{1}$, a similar series expansion

[TABLE]

is substituted into Eq. 4.15, and Eq. 4.12 is used to simplify the right hand side involving $S_{0}$. The resulting equation involves three extra nonlinear terms that require individual Fourier series expansions:

[TABLE]

After some further simplification, the ODEs for $S_{1}^{m}(z)$ can be written as

[TABLE]

which are identical to the $S_{0}^{m}(z)$ equations except that the right hand side is given by Eq. 4.12 as

[TABLE]

The corresponding boundary conditions from Eq. 4.16 are

[TABLE]

We have already obtained the first term ($m=0$) in the $S_{1}$-series as Eq. 4.18 from the steady state solution, while for $m\geqslant 1$ we proceed as before by splitting

[TABLE]

where the homogeneous solution is

[TABLE]

and the particular solution $S_{1,p}^{m}$ is identical to Eq. 4.31d with $H_{0}^{m}$ replaced by $H_{1}^{m}$. Finally, applying the boundary conditions yields the coefficients

[TABLE]

Although the formulas for $S_{0}^{m}$ and $S_{1}^{m}$ are somewhat complex, some insights can be drawn about the behaviour of solutions by concentrating on the underlying structure. It is clear from Eq. 4.31 that the leading order solution involves terms of the form $\Re[\exp(\varrho_{m}^{\pm}z)\exp(\mathrm{i}mt)]$, which when summed give rise to upward- and downward-travelling saturation waves moving at speed

[TABLE]

and with amplitude that decays with $z$ in the direction of travel. Furthermore, in the simple case of zero forcing (i.e., no transpiration and $E_{o}=0=\phi$) we can show that the leading order term in Eq. 4.12 undergoes a simple exponential decay process from the initial state over a relaxation time scale

[TABLE]

In order to test these observations on a concrete example, we consider a special time-varying transpiration source term $E(t)$ as pictured in Fig. 8a that begins at zero, jumps suddenly to a relatively large non-zero constant value (here $E_{o}$), and then after some delay returns to a “normal” diurnal periodic cycle (a pure $m=1$ mode). This source may be viewed physically as arising from a severe weather event or other rapid change in ambient conditions. To allow clear travelling waves to develop along the stem, we also concentrate the transpiration source term in a small region near the tree-top using the function $f(z)$ as shown in Fig. 8b.

Numerical simulations yield the time-dependent saturation profiles in Fig. 8c at two different stem heights. In response to the initial transpiration jump, there is a clear relaxation phase where saturation decays gradually back to a constant steady state. The zoomed plot in Fig. 8d shows that the relaxation time is roughly 3.5 h, which compares well with the estimate of $T_{r}\approx 4.06$ h obtained using Eq. 4.34, converted to dimensional time. Moving next to the diurnal variations in transpiration, the zoomed plot of the two saturation profiles in Fig. 8e demonstrates the existence of travelling waves of saturation moving down the stem, shown as a time shift between the two corresponding peaks. The shift can be estimated as roughly $1$ h, whereas the asymptotic wave speed formula Eq. 4.33 can be used to estimate a shift of $1.44$ h. A second simulation is shown for a larger value of $\eta=14.3$ in Fig. 8f, which exhibits a computed shift of 2.7 h as compared with the asymptotic estimate of 2.86 h, which is a significantly better agreement.

One further verification of the asymptotic solution is now performed in which we compare the radially-averaged saturation profile from simulations to that obtained from the two-term asymptotic expansion over a one-day period. Taking $E_{o}=1\times 10^{-9}$ so that $\mu\sim\phi\sim\zeta$, the relative difference between numerical and asymptotic values of the saturation deficit is $\mathcal{O}\!\left(10^{-4}\right)$. It is only when the transpiration rate is increased to $E_{o}=3.94\times 10^{-8}$ that we begin to move out of the asymptotic regime and differences in saturation become visible to the naked eye as shown in Fig. 9.

4.7 Discussion: Physical Relevance of Dimensionless Ratios

To conclude this section, we highlight four dimensionless parameters (or ratios) that play a prominent role in the asymptotic solution just derived, and which also have straightforward physical interpretations.

Stem aspect ratio, $\zeta$

This parameter plays a central role in the asymptotic analysis as the power series expansion parameter. One of our main conclusions is that the ratio of the radial and vertical velocity components is ${\bar{v}_{r}}/{\bar{v}_{z}}=\mathcal{O}\!\left(\zeta\right)$, which is small regardless of the conductivity ratio $\kappa$. Indeed, this is what allows us to assume $\kappa=1$ in our asymptotic derivation and still obtain a solution whose leading order term is relevant to the anisotropic case.

Transpiration–flux ratio, $\chi=\phi/\delta$

This ratio appears in several key places throughout the asymptotic derivation wherever transpiration terms appear in the solution, including Eqs. 4.7, 4.17 and 4.26. Substituting the parameter definitions into $\chi$ yields

[TABLE]

and based on the right-most expression $\chi$ may be interpreted as the ratio of transpiration flux through the stem surface to vertical sap flux through a circular stem cross-section. For the parameters of interest in this study $\chi\sim{1}$, which reflects the balance that must exist between these two fluxes under “normal” daytime conditions. Other limits could nonetheless be considered, such as $\chi\ll 1$ for which the transpiration rate is insufficient to generate an appreciable change in saturation and hence the impact of transpiration will only be felt in higher order terms. On the other hand, imposing a higher transpiration rate with $\chi\gg 1$ could be viewed as shifting the tree into an embolism regime for which saturation is no longer a smooth function, violating a fundamental assumption in our model.

Gravity–saturation ratio, $\xi=\mu/\delta$

This ratio also appears in the governing equations Eqs. 4.4 and 4.12 and the leading order saturation boundary conditions Eqs. 4.6 and 4.14. In terms of dimensional parameters, $\xi=\frac{nH}{\psi_{o}\delta}$, which can be viewed as a balance between the driving force due to gravity, and the corresponding (saturation-dependent) capillary forces acting on the pore scale in both stem and roots. At night when transpiration is a minimum these two forces must be in balance to prevent water loss into the soil, which is reflected in the fact that $\xi\sim{1}$. Recall that our analysis requires $\mu\sim\zeta$ (when $\psi_{1}\sim 1$, see Eq. 4.16) which places a restriction on the model parameters. For example, if tree height and radius are scaled up by the same factor so that $\zeta$ remains fixed then the ratio $\psi_{o}/n$ must also increase, meaning that larger trees may develop larger tensions for a given saturation deficit.

Time parameter, $\eta$

Consider the formulas for the travelling wave speed Eq. 4.33 and relaxation time Eq. 4.34 derived in the previous section. For simplicity, consider a tree with no taper ($\alpha=0$) in which case the wave speed formula reduces to $(2m/\eta)^{1/2}$ and the relaxation time to $T_{r}=4\eta/\pi^{2}$. Clearly, the parameter $\eta$ is intimately tied to the time variation of the solution both through the speed of propagation of saturation disturbances along the stem (with characteristic time proportional to $\eta^{1/2}$) and the time for decay of disturbances (proportional to $\eta$).

5 Conclusions

We have extended and generalized the 1D porous medium model for transpiration-driven tree sap flow developed by Chuang et al. in [7] to a tapered 3D axisymmetric stem geometry. Methods of asymptotic analysis are used to derive the first two terms in a regular asymptotic series expansion in powers of the aspect ratio $\zeta$, which are then used to obtain formulas for the spatial and temporal variations of the saturation and velocity components. Various possible flow regimes are studied through the use of several dimensionless ratios. The results are illustrated using a set of physical parameters and nonlinear coefficient functions corresponding to Norway spruce trees, although the analysis applies to a much more general class of parameters. One interesting conclusion of our analysis is that moderate levels of anisotropy $\kappa$ in the hydraulic conductivity do not induce any radial flow (or radial solution dependence) at leading order so that the simpler isotropic porous medium solution is a reasonable approximation for moderate $\kappa$. Indeed, introducing an anisotropy $\kappa\sim\zeta^{p}$ that depends on powers $p=0,1,2$ has the effect of successively “shifting” radial variations up to the higher order terms in the asymptotic solution.

The asymptotic results are verified using a second order finite volume approximation of the original governing equations, showing the results to be accurate for a relatively large range of saturations (as long as we avoid the low pressure regime where embolisms are likely to form, at which point the model assumptions are no longer valid). An interesting correspondence is drawn between asymptotic results for the steady state case and the more pervasive circuit model representation for tree sap flow. Furthermore, the vertical mass flux was shown to agree with experimental results from [7].

In future, we plan to extend the model to a more general non-symmetric 3D geometry where the solution either experiences angular variations or a more complicated branching distribution along the stem. This model also forms an ideal platform from which to study the interplay between transpiration and embolism formation under more extreme conditions.

Our asymptotic analysis will also facilitate the study of inverse problems related to estimating the transpiration functions $f(z)$ and $E(t)$, as well as model parameters such as $\mu$, $\eta$ and $\phi$. Given a set of noisy measurements of sap velocity, we can use our formulas to recover estimates of $f(z)$ and $E(t)$. Then if $\eta$ is small enough that saturation relaxes to the zero transpiration steady-state at night, we may estimate $\mu$ using Eq. 4.17. The parameter $\eta$ can then be determined using measurements of the relaxation time at night in Eq. 4.34. Finally, assuming that data for $K_{o}$ is available from vulnerability curve measurements, $\phi$ can be estimated making use of values for $f_{o}$ and $E_{o}$. The details of this estimation procedure can be found in [16].