## Abstract

Sap transport in trees has long fascinated scientists, and a vast literature exists on experimental and modelling studies of trees during the growing season when large negative stem pressures are generated by transpiration from leaves. Much less attention has been paid to winter months when trees are largely dormant but nonetheless continue to exhibit interesting flow behaviour. A prime example is sap exudation, which refers to the peculiar ability of sugar maple (*Acer saccharum*) and related species to generate positive stem pressure while in a leafless state. Experiments demonstrate that ambient temperatures must oscillate about the freezing point before significantly heightened stem pressures are observed, but the precise causes of exudation remain unresolved. The prevailing hypothesis attributes exudation to a physical process combining freeze–thaw and osmosis, which has some support from experimental studies but remains a subject of active debate. We address this knowledge gap by developing the first mathematical model for exudation, while also introducing several essential modifications to this hypothesis. We derive a multiscale model consisting of a nonlinear system of differential equations governing phase change and transport within wood cells, coupled to a suitably homogenized equation for temperature on the macroscale. Numerical simulations yield stem pressures that are consistent with experiments and provide convincing evidence that a purely physical mechanism is capable of capturing exudation.

## 1. Introduction

The study of tree sap flow has a long history that has given rise over time to the concept of the hydraulic architecture of trees [1]. Despite the extensive literature on this subject, several aspects of sap transport remain controversial, including the cohesion–tension theory of sap ascent [2–4]; embolism formation and recovery [5,6], which is ubiquitous in species subject to drought- or freezing-induced stresses; and sap exudation in maple and related species such as walnut, butternut and birch [7]. Furthermore, there is a great deal of current interest in the possible effects of recent changes in weather patterns on both individual trees and forest ecosystems [8,9], and their connections with sap hydraulics [10]. The problems just described involve complex interactions between sap flow and other phenomena such as nutrient transport, photosynthesis, soil physics, atmospheric dynamics, cell growth, etc. Despite the extensive work to date on mathematical and computational modelling of trees and their interactions with the environment, many open questions remain that can only be addressed by considering sap flow coupled with other processes and building integrated models that connect flow and structure at different spatial scales and levels of organization [11].

Sugar maple is a keystone species in the forests of central and eastern North America [12] and so is worthy of special attention. Members of the maple family are distinguished from other hardwoods by a number of unusual structural and functional features that allow them to exude sap during winter [13,14], to generate unusually high rates of nitrification [15] or to recover from freeze-induced embolism [16,17]. The potential impacts of climate change on maple have also attracted recent attention [12,15], motivated by the economic importance of the maple syrup industry, not to mention maple's high timber value. In particular, maple sap yields are sensitive to even small variations in temperature or snow cover during the harvest season, so that recent unusual weather patterns underscore the importance of developing a better understanding of the effects of local environmental conditions on sap flow [18,19].

Hundreds of scientific papers have addressed the phenomenon of sap exudation during winter when maple trees are leafless and yet still exhibit pressure variations that range over 150–180 kPa [7,14,20,21]. However, the precise mechanism driving the generation of heightened exudation pressure is still not fully understood [1]. The first systematic study appeared in an 1860 article by Sachs [22], who attributed exudation pressure to thermal expansion of gas within sapwood or *xylem*. The next major advance in understanding followed from the exhaustive study of Wiegand [23], who found to the contrary that thermal expansion of gas, water or wood has minimal impact on exudation. Instead, Wiegand proposed a vitalistic or ‘living cell’ hypothesis wherein sugar is released into the sap by some cellular activity, which gives rise to elevated pressure from osmotic gradients across selectively permeable membranes separating wood cells. Subsequently, this osmotic mechanism figured prominently in the literature, although experimental studies have continued to yield conflicting results that in turn stimulated development of new theories advocating alternate (bio-)physical mechanisms. For example, some authors continued to support the thermal expansion hypothesis [24], while others advocated the various roles of gas dissolution [25], cryostatic suction due to freezing [26] or temperature-induced changes in bark thickness [27]. More recent studies have led to a new understanding of exudation as a physical process deriving from a combination of freezing and thawing of sap [13] with osmosis [28]. Although some experimental evidence supports this hypothesis [7] the precise mechanisms behind the exudation phenomenon is still not fully understood.

We aim to resolve this long-standing open question by developing the first mathematical model for the freeze–thaw process in maple. We uncover the essential role played by two physical mechanisms whose significance has not yet been recognized—namely, root water uptake and freezing point depression (FPD) due to sap sugar content. Using numerical simulations of repeated freeze–thaw cycles, we obtain computed exudation pressures that are consistent with experimental results.

Although the focus of this paper is on developing a complete and physically consistent model for sap exudation, our results also have more far-reaching consequences. This work affords new insights into the complex multi-physics processes occurring in trees and also provides a framework for studying other practical questions of importance to tree physiologists and maple syrup producers. Our model also provides a platform for studying related phenomena such as embolism that occur in a much broader range of tree species, as well as evaluating the response of trees to changes in environmental variables such as temperature and soil moisture arising from various climate change scenarios.

## 2. Physical mechanism for sap exudation

### 2.1. Milburn and O'Malley's hypothesis

Experimental work up to the 1980s demonstrated that no single physical mechanism is capable of capturing measured winter stem pressures [29], and sap exudation remained an unsolved puzzle until the ground-breaking study of Milburn & O'Malley [13]. They proposed a physical mechanism based on freezing and thawing of sap, motivated by the unique structural characteristic of xylem in maple (and related trees) that a significant proportion of the *libriform fibres* (or simply fibres) are primarily gas-filled rather than being liquid-filled as in most other hardwood species [23]. This peculiar feature of the fibres should be contrasted with the two other cell types that play an active role in sap transport—*vessels* and *tracheids*—which are mostly sap-filled and are connected hydraulically to each other via paired pits (figure 1*a*). Indeed, recent experiments [7,30] suggest that fibres are essentially non-conductive in comparison with other xylem elements because they lack end-to-end cell connections and their lateral walls contain mostly unpaired (blind) pits that are smaller and fewer in number than the more conductive vessels and tracheids.

Milburn and O'Malley focused on the dynamics of a single fibre–vessel pair as shown in figure 1*b*, and ignored other xylem elements such as parenchyma and ray cells. Although fibres had previously been thought to play a purely structural role, Milburn and O'Malley proposed that as temperature falls below freezing, liquid from the vessel is drawn by cryostatic suction through the porous fibre–vessel wall to freeze on the inside of the fibre (figure 1, stages 1–2–3). As a result, any gas contained within the fibre is compressed and acts as a pressure reservoir. When temperature subsequently rises and ice thaws, the process reverses and the compressed gas bubble forces melt-water into the vessel, thereby re-pressurizing the vessel sap (figure 1, stages 3–4–1).

### 2.2. Tyree's modified hypothesis, with gas dissolution and osmosis

This freeze–thaw hypothesis was critically evaluated by Tyree [28], who proposed a modified hypothesis featuring two important additions. First, he recognized that gas under pressure will dissolve within an adjacent liquid [31] and that pressures encountered in maple xylem are high enough that any bubbles should dissolve completely given sufficient time [32]. Therefore, some additional mechanism is required to sustain gas bubbles in the fibres. Tyree's second observation was that measured xylem pressures depend strongly on sugar concentration in the vessel sap, 80% of which derives from sucrose [21], which led him to conclude that sucrose is required for exudation. He recognized that although the axial conductivity of fibres is negligible in comparison with vessels and tracheids, the lignified cellulose making up secondary cell walls should admit a small radial conductivity. He then hypothesized that the fibre–vessel wall forms an osmotic barrier that allows water to penetrate but not the larger sucrose molecules. Consequently, an additional osmotic pressure difference exists between the sweet vessel sap and pure fibre water, which he argued is responsible for preventing fibre gas bubbles from completely dissolving.

Tyree's modified freeze–thaw hypothesis includes water phase transitions, gas dissolution and osmosis, and is currently the prevailing hypothesis for sap exudation [28]. It depends strongly on the existence of a hydraulically isolated system of fibres and a selectively permeable fibre–vessel wall, both of which have since been confirmed experimentally [7]. Although this evidence is compelling, there has been no attempt yet to model this process mathematically (except for a related process without phase change in the context of embolism recovery in maple [17]), and so it remains unclear whether this physical description is capable of capturing exudation.

### 2.3. Three essential physical mechanisms

Before proceeding further, we extend the freeze–thaw hypothesis just described by incorporating three additional mechanisms described in §2.3.1.–2.3.3.

#### 2.3.1. Gas bubbles in the vessel

Sap (like water) is an incompressible fluid so that in the rigid, closed vessel network of a leafless tree there is no mechanism for fibre–vessel mass transfer if vessels are completely saturated with sap. However, the existence of gas bubbles within vessels is well documented in maple [23,33] and other hardwood species [28]. Even if xylem pressures were high enough to dissolve such bubbles, gas would eventually be forced out of solution upon freezing and so at least a transient presence of gas bubbles is unavoidable. Therefore, introducing a gas phase in the vessels provides a plausible mechanism for fibre–vessel pressure exchange.

#### 2.3.2. Sap freezing point depression

Sap contains dissolved sugars and hence experiences a reduced freezing point compared to pure water according to Blagden's law [34], where *K*_{b} is the cryoscopic constant, *C*_{s} is sugar concentration and *ρ*_{w} is water density. For example, sap containing 3% sucrose by mass experiences an FPD of . Although this temperature difference may appear insignificant, we will see that it is actually large when considered on the scale of individual cells, and indeed is sufficient to account for the existence of ice in fibres while sap in adjacent vessels remains in liquid form. This partitioning of ice and liquid in neighbouring fibre–vessel pairs induces cryostatic suction that draws liquid out of the vessel to form ice on the inner fibre wall.

#### 2.3.3. Root water uptake during freezing

No previous hypothesis for sap exudation explicitly considers the role of root water uptake. Furthermore, several studies suggest that root pressure in maple has a negligible effect on exudation [20,35]. Nonetheless, it is well known that during winter months trees can draw water from the roots if soil temperatures are high enough [1,26], which can be caused by an insulating snow cover [36]. Recent experiments on maple saplings [37,38] have provided the first direct evidence that root water uptake occurs in maple during winter while exudation is underway.

Our aim is now to incorporate these three modifications into a model of the freeze–thaw process outlined previously, and then demonstrate that the resulting equations are capable of reproducing observed behaviours.

## 3. Mathematical formulation

### 3.1. Outline of the modelling approach

The freeze–thaw mechanism outlined in the previous section involves processes operating on two distinct spatial scales: the *microscale* corresponding to individual wood cells with dimensions ranging from 10 to 100 µm; and the *macroscale* corresponding to the tree stem with a diameter range of tens of centimetres. The derivation of our mathematical model for sap exudation therefore divides naturally over these two scales. Firstly, we develop microscale equations that capture cell-level processes within libriform fibres and vessels, combining the dynamics of freezing, thawing, gas dissolution, osmotic pressure, heat transport and porous flow through the fibre–vessel wall. Secondly, we consider heat transport in the entire tree stem and apply periodic homogenization to derive an equation for the macroscale temperature that incorporates microscale cellular processes via appropriately defined transport coefficients and source terms. We proceed as follows:

— Start from an existing microscale model for the thawing half of the freeze–thaw process [39,40] in which a two-dimensional periodic microstructure is assembled from copies of a

*reference cell Y*containing a single fibre and vessel (figure 2*a*). The fibre is placed at the centre of*Y*(where the dashed line denotes the fibre–vessel wall) and the remainder of the reference cell corresponds to the vessel. This choice of geometry is a mathematical idealization that captures the volumes of the fibre and vessel compartments, but is not intended to accurately represent the actual layout of wood cells. The governing equations consist of a partial differential equation (PDE) for the microscale temperature along with five ordinary differential equations (ODEs) for phase interface locations and root water volume, coupled nonlinearly through source terms and algebraic constitutive relations.— Supplement the thawing model with analogous equations for the freezing process, which have similar structure but differ slightly depending on the precise state of freezing or thawing in the fibres and vessels.

— Apply periodic homogenization [40,41] to derive a macroscopic equation for temperature that is coupled to the microscale (reference cell) problem at each point within the tree stem (figure 3

*b*). The macroscopic heat diffusion equation contains an integral source term depending on the microscale temperature and capturing all processes on the cellular level. A similar approach has been applied to studying protein-mediated transport of water and solutes in non-woody plant tissues [42].— Exploit radial symmetry on the micro- and macroscales to reduce both PDEs for temperature to a single spatial dimension. We will see later in §4.1 that the microscale equations need only be solved on the circular sub-region

*Y*^{2}in figure 3*a*(consisting of the fibre and surrounding vessel overlap region) which is clearly radially symmetric.

### 3.2. Microscale equations for cell-level thawing process

The cell-level model is based on equations already developed for the thawing half of the freeze–thaw process by Ceseri & Stockie [39], which were subsequently homogenized by Graf & Stockie [40]. We therefore begin by considering an intermediate state in the thawing process corresponding to stage 4 in figure 2, during which the vessel sap is completely thawed while the fibre contains both liquid and ice. We extend the Ceseri–Stockie model by incorporating additional physical effects that capture the influence of ice–water surface tension, root water uptake and volume change due to ice/water phase transitions. We discuss some of the most important assumptions and modifications next, leaving the reader to consult the references [39,40] for a complete derivation and discussion of assumptions.

Our model is based on the conceptual diagram in figure 1*b* that depicts a single vessel–fibre pair. Tracheids are not treated separately but instead ‘lumped together’ with vessels because, although they are connected hydraulically to vessels via paired pits, they have a much smaller diameter and correspondingly lesser influence on sap transport than vessels. Because multiple fibres adjoin and interact hydraulically with each vessel, we introduce the parameter *N*^{f} representing an average number of fibres per vessel, which is estimated from SEM images [7] as *N*^{f} ≈ 16. Our model captures the dynamics of a single fibre and then scales all fibre–vessel flux terms by an appropriate factor of *N*^{f}.

We assume that sapwood can be represented as a doubly periodic array of idealized reference cells *Y* as shown in figure 3, where each reference cell contains a circular fibre embedded within a surrounding square liquid region representing the adjoining vessel. This choice of geometry is made for mathematical convenience in the homogenization step, and can be justified because our aim is to derive a system of equations that captures the net effect of sap flow and heat transport on the microscale, keeping in mind that any specific geometric details will ultimately be ‘averaged out’ during the homogenization process anyways.

Our two-dimensional geometry comes with the built-in assumption that axial (vertical) variations are neglected. In the absence of root water uptake, the model tree behaves as a closed system that is essentially in equilibrium. Any pressure differences initiated by phase change engender primarily horizontal flow between neighbouring cells, and negligible axial flow. Furthermore, we have already shown [39] that phase change on the microscale dominates the pressure exchange process and occurs very rapidly (of the order of milliseconds). Root water uptake induces an axial flow but this is a much slower process; therefore, over the timescales that dominate the microscale problem, axial transients may be neglected.

The fibre is a circular cylinder with length *L*^{f} and cross-sectional radius *R*^{f} as shown in figure 3*a*. Situated at the centre of the fibre is a cylindrical gas bubble with time-varying radius *s*_{g}(*t*), outside of which lies an annular ice layer with outer radius *s*_{iw}(*t*). The remaining volume extending to the fibre radius *R*^{f} contains melt-water from thawed ice. We note that this configuration is specific to the thawing process, and the ordering of ice and water layers would be reversed during freezing.

The vessel is represented by the portion of the reference cell lying outside the fibre–vessel wall (denoted by a dashed line) and the side length *ℓ* of the reference cell is chosen so that the vessel cross-sectional area equals that of a cylinder of radius *R*^{v}. Keeping in mind that there are actually *N*^{f} fibres connected to each vessel, we require that *ℓ* satisfy the area constraint
3.1

Within the vessel is a gas bubble of radius *r*(*t*), which is surrounded by liquid sap owing to the FPD effect that lowers the freezing temperature below that in the fibre. The cumulative volume of melt-water flowing through the porous fibre–vessel wall is denoted by *U*(*t*) and is measured positive from fibre to vessel. The final variable that determines the local state of the fibre–vessel system is the volume of root water uptake, denoted *U*_{root}(*t*).

We may now formulate a first-order system of five ODEs describing the time evolution of *s*_{iw}, *s*_{g}, *r*, *U* and *U*_{root}. The fibre ice–water interface is governed by the Stefan condition [43,44]
3.2where represents the normal temperature derivative on the interface (i.e. the curve along which temperature equals the melting (or freezing) point *T*_{m}) and the final term accounts for the volume of water transferred between fibre and vessel. This form of the Stefan condition assumes that liquid motion induced by phase density differences is negligible [43]. The microscale temperature *T*_{2}(*y, t*) is obtained as the solution of a heat diffusion equation that will be stated in the next section, where the microscale spatial coordinate is *y*. The parameters *ρ*_{w} and *k*_{w} denote the density and thermal conductivity of liquid water, while (*E*_{w}−*E*_{i}) is the enthalpy difference between water and ice (also called the latent heat or enthalpy of fusion) at locations where *T*_{2} = *T*_{m}. The effects of thermal expansion are known to be relatively small [23] and so have been neglected here.

Imposing mass conservation yields an equation for the fibre gas bubble radius (which in this thawing scenario is a gas–ice interface)
3.3where *ρ*_{i} is the density of ice. An equation for the vessel gas bubble radius follows from a similar mass conservation argument:
3.4where *L*^{v} denotes the length of a vessel. This last equation expresses the balance between water flux from neighbouring fibres and the slight volume change stemming from the water/ice density difference. The effect of gas dissolution has been omitted here but will be incorporated below in the gas density; this approximation was already justified in [45], which showed that incorporating dissolution in these equations has negligible impact on the bubble radii *s*_{g} and *r*.

Darcy's law governs liquid water flux through the porous fibre–vessel wall
3.5where the wall is characterized by hydraulic conductivity and surface area *A*. The pressure term in square parentheses derives from four contributions: liquid pressure in the vessel () and fibre (), osmotic pressure (*p*_{osm}) and capillary pressure () due to ice–water surface tension [46]. This latter contribution, also known as cryostatic suction, follows hand-in-hand with FPD and arises whenever ice lies on the inside surface of the wall and liquid sap is present on the vessel side, since then the small capillary pores in the adjoining wall (with radius *r*_{cap}) contain both ice and liquid. For the thawing scenario under consideration here, water lies on both sides of the fibre–vessel wall and so ; however, other stages in the freeze–thaw process can give rise to non-zero as detailed in the next section (see also figure 4).

The final ODE comes from another application of Darcy's law to root flux
3.6where is the root hydraulic conductivity and *A*_{r} denotes the portion of root surface area corresponding to a single vessel. The cut-off function ‘’ ensures that water only flows inward from soil to roots and not outward, which is consistent with experiments that demonstrate root outflow can be a factor of five smaller than that for inflow [47]. Indeed, studies of root water transport in a variety of tree species show that root conductivity can vary with factors such as temperature [48], root age [49], and time of day [47] or season [50]. Many authors attribute this selective control of water transport to membrane proteins known as aquaporins [51].

In the preceding discussion, we introduced a number of constant parameters whose values are listed in table 1. The remaining symbols correspond to intermediate variables whose definitions we provide next. First, the density of gas in the fibre and vessel bubbles depends on initial values of density and volume, modified to account for dissolved gas according to
3.7where *H* is the dimensionless Henry's constant for air in water. The various phase volumes are determined from the cylindrical cell geometry as
3.8
3.9

The corresponding gas pressures are given by the ideal gas law as
3.10where is the universal gas constant and *M*_{g} is the molar mass of air. The water and gas pressures in both fibre and vessel differ by an amount equal to the capillary pressure, which is determined by the Young–Laplace equation as
3.11where *σ*_{gw} is the air–water surface tension. The osmotic pressure across the fibre–vessel wall depends on the sap sugar concentration according to
3.12

Finally, the sap sugar content induces a reduction in freezing temperature that obeys 3.13

### 3.3. Equations for other phase transitions

In the previous section, we developed equations specific to the thawing process, during which the vessel is completely thawed and the fibre contains a mix of gas, water and ice (see stage 4 in figure 2). We describe next how these equations should be modified to capture other freeze–thaw states in the fibre and vessel. In particular, we account for the fact that phase interfaces can appear or disappear whenever ice completely thaws (or liquid completely freezes), as well as the reversal of the ice and water layers in the fibre during freezing and thawing.

In fact, many equations remain unchanged throughout the entire freeze–thaw cycle, with the exception being those for *s*_{g}, *s*_{iw} and . The required modifications for each case are listed in figure 4, referenced by the numbered stages in figure 2. We emphasize that the ice–water interface lies within pores in the fibre–vessel wall and forms a *mushy layer* wherein both solid and liquid phases coexist in the pore space. This type of phase interface (called a *frozen fringe* in the context of ice lensing in soils [46]) differs from an idealized gas–water interface in that the interfacial pressure jump increases with ice volume fraction according to
where *σ*_{iw} represents the ice–water surface tension and are the corresponding volume fractions. Finally, we note that when both fibre and vessel are completely frozen (stage 3) the equations for *r*, *U* and *U*_{root} also drop out of the system.

### 3.4. Homogenized equation for temperature

The equations derived in the preceding two sections govern microscale processes at the cellular level, whereas on the macroscale the temperature is of primary interest, and it is transport of heat between the external (ambient) environment and the interior of the tree stem that drives the freeze–thaw process. Clearly, there exists a two-way interaction between the global temperature and the local fibre–vessel state, wherein temperature governs phase change dynamics in fibres and vessels, while cellular processes in turn influence heat transport through the Stefan condition and local phase volume fractions. To simplify this complex multiscale problem, we exploit a separation in spatial scales reflected in the fact that state variables describing the fibre–vessel configuration are essentially ‘invisible’ on the macroscale except through their effect on heat transport properties of the sap- and gas-filled wood.

Because of the repeating microstructure of wood, this problem is ideally suited to the application of *periodic homogenization*. The philosophy behind this approach is to solve at each point in space a local problem on a reference cell *Y* that determines the solution state on the microscale. By using an appropriate homogenization or averaging procedure, the effect of microscale variables on the macroscale may then be incorporated into equations for the global solution variables. One technical requirement is that the reference cell must divide into two sub-regions, , separated according to whether heat diffusion is fast (in *Y*^{1}, the outer portion of the vessel) or slow (in *Y*^{2}, an overlap region covering the fibre and the remainder of the vessel). The result is two heat equations: one governing temperature on the macroscopic domain *Ω* and the second on *Y*^{2} × *Ω*. When these two equations are coupled together, we obtain a two-scale temperature solution on the domain *Y* × *Ω*. Instead of fully coupling the micro- and macroscale equations, this homogenization approach leads naturally to a simpler system of equations that captures the essential aspects of coupling between scales. A similar homogenization approach has been applied by Chavarra-Krauser and Ptashnyk to a model of water and solute transport in plants [42].

The dynamics of heat transport are best described using a mixed formulation written in terms of both temperature and specific enthalpy, which are denoted, respectively, by *T*_{2}(*x*, *y*, *t*) and *E*_{2}(*x*, *y*, *t*) on the reference cell region *Y*^{2}, and *T*_{1}(*x*, *t*) and *E*_{1}(*x*, *t*) on the macroscale. The variables *T*_{1} and *E*_{1} depend on time *t* and the global spatial coordinate *x*, whereas microscale quantities have an additional dependence on the reference cell *Y* through a local spatial coordinate *y*. Temperature and enthalpy are not independent variables but instead are related via the piecewise linear function
3.32

We introduce the large parameter (taking in practice) to impose a small but non-zero slope on the central plateau region where *T* ≈ *T*_{m}. We also make use of the fact that and choose
3.33so that the function *T*(*E*) is continuous. This form of *T*(*E*) is a regularization of the exact temperature–enthalpy relationship [56] that avoids numerical instabilities in the calculation of temperature and also recovers the exact (piecewise linear) result in the limit as and *δ*_{i}, *δ*_{w} → 0.

During the homogenization procedure [40], we find that heat transport in the reference cell must only be treated on the sub-region *Y*^{2} where temperature obeys
3.34and *D*(*E*_{2}) is a thermal diffusion coefficient that is a piecewise linear and continuous function of enthalpy [56]
3.35

We employ this non-standard definition of *D* (instead of the usual thermal diffusivity having units m^{2} s^{−1}) so that we can factor out the specific heat, thereby allowing the same coefficient to be used in both this microscale heat equation and the mixed temperature–enthalpy form we develop below for the macroscale. We include an explicit time- and global space-dependence in *Y*^{2}(*x*, *t*) to emphasize the fact that the ice region within the fibre is bounded by a moving water–ice interface, and that the fibre configuration varies from point to point throughout the tree stem. On the water–ice interface (corresponding to the inner boundary of *Y*^{2}) the temperature equals the melting point value
3.36

We thereby obtain the macroscale temperature equation
3.37where the coupling with microscale variables is embodied in a surface integral term. The factor multiplying the diffusion coefficient is a 2 × 2 matrix whose entries depend on the reference cell geometry according to
3.38for *i*, *j* = 1, 2. Here, *δ _{ij}* is the Kronecker delta symbol and

*μ*(

_{i}*y*) are solutions of a standard reference cell problem on

*Y*

^{1}[41]. The temperature on the outer surface of the tree is held at the ambient value 3.39

Finally, the micro- and macroscale solutions are coupled by matching temperature on the interior boundary 3.40

In summary, the governing equations consist of a system of differential–algebraic equations (3.2)–(3.13) and (3.34)–(3.36) for the microscale temperature and fibre–vessel state variables within each local reference cell. These are supplemented by equations (3.37)–(3.40) for the macroscale temperature on *Ω*. Both problems are solved at each spatial point and the two solutions are coupled by means of the integral source term in (3.37) and the boundary condition (3.40). The geometry of the local reference cell is also incorporated into the macroscale problem via the (constant) pre-factors multiplying the diffusion coefficient in (3.37).

## 4. Simulating daily freeze–thaw cycles

### 4.1. Numerical solution algorithm

The radial symmetry of both micro- and macroscale domains implies that all solution variables can be written as functions of a single radial coordinate and time. We use a *method of lines* approach and discretize the temperature variables in space using finite elements, yielding a large system of time-dependent ODEs. When combined with the ODEs and algebraic equations governing microscale fibre–vessel dynamics, the resulting coupled system is integrated in time using a standard ODE solver. The spatial discretization on the two scales proceeds as follows:

—

*Microscale*(*cell-level*)*equations*. The fibre ice temperature is assumed to be a uniform 0°C, and gas temperature is also taken as constant since the thermal diffusivity of gas is so much larger than that for either ice or water. Therefore, the PDE (3.34) for temperature on*Y*^{2}must only be solved on the annular region between*Γ*and the phase interface*s*_{iw}(figure 3*a*). We find that sufficient accuracy is obtained for*T*_{2}by using only four radial grid points within the annulus. Because the phase interface evolves in time, we use a moving mesh approach wherein the motion of grid points introduces an additional ‘grid advection’ term that is proportional to the mesh point velocity [57].—

*Macroscale*(*tree-level*)*equation*. The tree stem is similarly divided into equally spaced radial points, and here we find that taking 20 grid points yields sufficient accuracy in*T*_{1}. Owing to radial symmetry, the integral source term in (3.37) reduces to multiplication by the curve length . The factors depend only on the reference cell geometry and so can be pre-computed at the beginning of a simulation.

We employ an efficient split-step approach where in each time step the reference cell problem is solved for the microscale variables, and then the macroscale temperature equation is solved by holding the microscale variables constant.

The algorithm described above has been implemented in Matlab using the built-in stiff ODE solver `ode15s` to integrate the equations in time. The only algorithmic detail remaining to be described is the switching between equations required as phase interfaces appear or disappear. We can capture this switching simply and robustly using the `Events` option provided in the ODE solver suite, which signals an event based on zero-crossings of an ‘indicator function’. During any portion of the freeze–thaw cycle, the indicator function is set equal to either the thickness of a phase interface or the difference between the phase temperature and the melting temperature. When the indicator crosses zero, the time integration halts, equations are modified appropriately, and the ODE solver is restarted using the new set of equations and taking the current solution as the new initial state. The time integration then proceeds until the next phase change event is signalled. A typical simulation covering four daily temperature cycles requires between 30 and 45 min of clock time on an Apple MacBook Pro with 2.3 GHz quad-core Intel i7 processor.

### 4.2. Choice of parameters

The algorithm just described is used to simulate freeze–thaw dynamics in a typical *base case* scenario for which all parameters are listed in table 1. We take a ‘sapling’ of diameter 0.07 m consisting entirely of sapwood. The sugar content of maple sap ranges from 1 to 5% by mass [14] and so we choose a representative value of 3% that induces a vessel FPD of . To mimic temperature variations during late winter, we let ambient temperature vary sinusoidally between −10°C and +20°C over a 24-h period (this range is somewhat extreme but is chosen to correspond with the experiments of Améglio *et al.* [20] that we will describe shortly). We begin with a freezing event and initialize the tree in a thawed state with uniform temperature 0.35°C, just slightly above the freezing point. Each fibre initially contains gas and water with 75% gas by volume, whereas the vessel has a much smaller initial gas content of 8%.

There remain two parameters whose values we have not been able to obtain reasonable estimates from the literature—root hydraulic conductivity and capillary pore radius *r*_{cap}—and so we have had to adjust their values in order to match numerical results with experimental data. First, we choose so that pressure and root uptake vary over timescales similar to those observed in experiments [14,20]. Then, we take *r*_{cap} = 2.8 × 10^{−7} m so that the exudation pressure build-up is within the observed range of 80–150 kPa [7,21]. This pore size is also consistent with that measured in other membranes that hinder transport of sucrose molecules [55].

### 4.3. Base case: pressure build-up during temperature cycling

Using these base case parameters and initial conditions, we perform two numerical simulations: one with root water uptake corresponding to a soil pressure of *p*_{soil} = 203 kPa, and a second with no root uptake (e.g. consistent with a completely frozen soil). Vessel sap pressures are compared in figure 5*a*, and in both cases we observe a periodic variation in pressure synchronized with daily temperature fluctuations. Without root uptake, the vessel pressure simply oscillates between two fixed values of 20 and 200 kPa, and there is no pressure build-up over multiple freeze–thaw cycles. However, when root uptake is included there is a gradual pressure increase superimposed on the background oscillations, with a total increase (measured from the local maximum in each cycle) of roughly 80 kPa over the 4 days. The accompanying plot of total root uptake in figure 5*a* shows that the majority of root water is absorbed during the first freeze–thaw cycle, followed by a more gradual uptake that is essentially complete after 3 days.

We next draw a direct comparison with the experiments of Améglio *et al.* [20] who studied black walnut trees (*Juglans nigra*) in a controlled laboratory setting where the living stump of an excised tree branch was connected via a sealed pipe to a pressure transducer. We calculate vessel sap pressure in our simulations as an average pressure across the stem cross section to be as close as possible to such a transducer measurement. We are unaware of any comparable data for sugar maple, but we claim that a meaningful comparison may still be drawn with Améglio's results since black walnut is closely related to maple and undergoes exudation under similar conditions [7,20,58]. The curves to focus on from Améglio's study (figure 5*b*) are the air temperature (labelled ‘*T* trunk’) and stem pressure (labelled ‘*P* control’).

The qualitative agreement between simulated and experimental pressures is remarkable considering the complexity of the processes involved and the minimal parameter fitting required. The overall shape of pressure curves is similar, with each freeze–thaw cycle exhibiting a rapid increase whenever temperature exceeds the freezing point. The pressure then attains a maximum, after which there is a slight decrease over roughly 6–8 h, followed by a rapid drop as ambient temperature crosses the freezing point again. We remark that there is also a rough quantitative agreement between simulations and experiments in that pressure oscillations have an amplitude of 80–100 kPa, and the total pressure build-up over 4 days also is similar. On the other hand, the maximum value of our simulated pressure is 290 kPa, which is almost double the 160 kPa observed in the black walnut experiments; however, it is possible that more time is needed for the experiment to reach steady state, not to mention that there are species-specific differences that could influence pressure.

A more quantitative comparison can be drawn based on two characteristic features of the pressure in figure 5*a* labelled as Δ*p*_{1} and Δ*p*_{2}. The first corresponds to the amplitude of oscillations in the absence of root uptake, which derives mainly from cryostatic suction and so can be estimated using the formula This value is close to the computed amplitude of the ‘no root uptake’ curve, as well as to the rise in vessel pressure during the initial thawing event for the ‘base case’. The second feature Δ*p*_{2} captures the exudation pressure build-up during the first freeze–thaw cycle which arises mainly from root water uptake. Because this additional water acts to compress the gas in fibre and vessel, we apply the differential form of the ideal gas law at constant temperature, , during the first freezing event. Substituting values of *p* ≈ 200 kPa for the initial vessel pressure, Δ*V* *≈* 0.4 cm^{3} for the root water volume uptake (taken from figure 5*a*) and *V* ≈ 1.15 cm^{3} for the initial gas volume in a slice through the tree cross section (with thickness equal to that of the reference cell, *L*^{f}), we obtain . The correspondence between this estimate and the computed value of 50 kPa is reasonable, considering that it ignores effects such as gas dissolution.

Despite the abundance of experimental data available for sugar maple [7,14,21,59], most experiments measure sap outflux from tapped trees [21] rather than the ‘closed system’ corresponding to an untapped tree that we consider here. Other measurements have been taken of excised wood samples rather than living trees, while yet others were taken in uncontrolled external conditions with irregular variations in pressure and ambient temperature. Consequently, we hesitate to attempt a detailed comparison between any of these experiments and our simulations; nonetheless, we can still draw a few quantitative comparisons. For instance, Tyree [14] performed experiments on excised maple branches that absorbed water at a maximum rate of 12 cm^{3} h^{−1}; for similar sized branches, our model yields a comparable maximum absorption rate of roughly 13 cm^{3} h^{−1} as well as qualitatively similar solution profiles. Another experiment by Johnson *et al.* [59] yielded a total root uptake of 2.0 cm^{3} during freezing, followed by a much smaller uptake of 0.1 cm^{3} during a subsequent freezing event. We see similar qualitative behaviour in our simulations, as well as measuring 2.2 and 0.2 cm^{3} of water absorbed during the first and second freeze, respectively.

### 4.4. Two crucial mechanisms: root water uptake and freezing point depression

To evaluate the relative importance of the various physical mechanisms, we present in figure 6 the base case pressure and root water uptake alongside simulations with each of the following mechanisms ‘turned off’: FPD, root uptake (repeated from figure 5*a*), osmosis and gas dissolution. The first two mechanisms clearly have the greatest impact on the build-up of exudation pressure. We have already discussed the crucial role of root uptake in facilitating pressure accumulation over multiple freeze–thaw cycles. This effect is underscored by the plots in figure 7*a* depicting the pressure response when root conductivity varies between zero (no uptake) and nearly 10 times the base value.

Without FPD, the vessel pressure remains nearly constant and there is minimal root uptake, whereas without osmosis the vessel pressure increases. We therefore conclude that the predominant impact of sugar on exudation is through FPD rather than osmosis, and even though is small it nonetheless plays a critical role in facilitating pressure transfer between fibre and vessel. This dependence is illustrated further in figure 7*b*, where the sugar content is varied between 0 and 7%, and we observe that both net pressure build-up and oscillation amplitude increase with sugar content.

One assumption requiring further investigation is that of zero conductivity to root outflow in (3.6), which we motivated by citing experimental results that exhibit a small but still non-zero root outflow [47]. To study this outflow effect, we take four different outflow conductivities equal to the inflow value scaled by a factor between 0 and 1 (where 0 corresponds to the base case). The results in figure 8 clearly show that allowing even a small outflow has a major influence on the root water uptake by preventing accumulation of water over multiple freeze–thaw cycles and thereby reducing build-up of exudation pressure. Because of the obvious sensitivity of these results to root outflow, a more extensive experimental study of root conductivity in maple is warranted.

We end this section by addressing the seemingly counterintuitive result in figure 6*a* that introducing osmosis decreases vessel sap pressure. This result can be most easily explained by considering the water flux equation (3.5) over a long enough time that the fibre and vessel have reached a quasi-steady state and . Then (3.5) reduces to the simple pressure balance

The ice–water capillary pressure is a constant, and our simulations show that osmosis has relatively small impact on fibre bubble size and pressure (the latter effect was discussed in [39]). Therefore, the primary influence of osmosis is within the vessel: osmotically driven flow from fibre to vessel compresses the vessel bubble, which not only increases the vessel gas pressure , but also increases the capillary pressure term (via a reduction in bubble radius *r*). The contribution from surface tension dominates and so the net effect is actually a *decrease* in vessel sap pressure , which is consistent with figure 6*a* and the results reported in [39].

### 4.5. Phase change dynamics on the microscale

When a completely frozen tree warms above 0°C during the day, a thawing front develops near the bark (wherein water and sap are frozen ahead of the front and thawed behind) and advances into the stem; an analogous scenario occurs upon freezing. Clearly, the ‘interesting’ solution dynamics will occur in the vicinity of this front, and hence knowledge of phase change on the microscale is desirable for understanding solution behaviour. Over a century ago, Wiegand [23] recognized the existence of freezing and thawing fronts that ‘penetrate the wood in a wave-like manner’ and in which ‘but few cells would actually take part in the production of pressure at any one time’; however, there has so far been no attempt to develop a mathematical model for this phenomenon. In particular, the role of FPD in governing the progress of these phase transitions throughout the sapwood has not been investigated before.

Phase change dynamics are most easily studied by means of a temperature–enthalpy diagram as depicted in figure 9*a,c*, which are each taken at a fixed time during a freeze or thaw event. Both plots feature a plateau region at the melting temperature, which has a horizontal extent equal to the *enthalpy of fusion*. Note that there are two distinct melting temperatures in fibre and vessel equal to *T*_{m} and respectively. The corresponding plots of temperature versus radius are shown in figure 9*b,d* which depict the local state of each point within the tree stem. For example, in the freezing case (top) the three grid points closest to the stem centre are completely thawed (red), the outermost point is frozen (blue) and the intervening points are undergoing freezing (purple). Owing to FPD, water in the fibre freezes before the vessel sap, thereby introducing a time delay in formation of ice between the fibre and vessel.

For both freezing and thawing, the bulk of the stem is in a state located at the leading edge of the enthalpy plateau (right edge for freezing, left edge for thawing). This behaviour can be explained by considering the local rate of phase change: conservation of energy at a phase interface is expressed mathematically using the well-known Stefan condition, which states that the rate of freezing (or thawing) is proportional to the temperature gradient. Referring to figure 9*b,d*, the temperature difference between adjacent points is smaller near the tree centre and larger near the bark. Consequently, at any location in the tree a freezing event begins within the fibre as a slow process, followed at a later time in the vessel which freezes relatively quickly. By contrast, a thawing event begins with a slow thawing of the vessel sap, followed by rapid thawing in the fibre.

## 5. Concluding remarks

We have developed the first complete mathematical model for the tree sap exudation process based on a prevailing freeze–thaw hypothesis. We introduced a number of additions to this hypothesis, and identified root water uptake and FPD as the two main driving mechanisms for sap exudation. In particular, we showed that the primary mechanism whereby sugar induces exudation pressure is via the FPD and not osmosis as was previously believed. Numerical simulations of the governing equations demonstrate qualitative and quantitative agreement with experimental data on sugar maple and the related species black walnut. The quality of agreement is striking considering that the model parameters were determined using a minimum of parameter fitting. Our work clearly demonstrates the need for further experiments on sugar maple that parallel the work of Améglio *et al.* on walnut [20]. Our model results lead to the important conclusion that FPD is a primary driver of sap exudation, which also requires experimental validation. Furthermore, because we have only rough estimates at present for two of the model inputs—capillary pore size *r*_{cap} and root conductivity (especially differences between conductivity to inflow and outflow)—more accurate measurements of these parameters are also required.

Our model provides an ideal platform from which to investigate other problems related to sap flow in maple and related species. First of all, we aim to extend our current model of a two-dimensional stem cross section to three dimensions. This will permit us to incorporate variations in gravitational pressure head and sugar concentration with height [23] and to study problems of practical importance to the maple syrup industry such as optimizing tap-hole placement or determining sensitivity to changes in soil or climatic conditions. Finally, there are a number of intriguing parallels between exudation and the phenomenon of freeze-induced winter embolism [16,17] that are also worthy of future investigation.

## Authors' contributions

J.M.S. designed the study. I.G. and M.C. designed the numerical algorithm. I.G. carried out the computational studies. All authors derived the mathematical model, analysed the results and wrote the manuscript. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by research grants from the Natural Sciences and Engineering Research Council of Canada and the North American Maple Syrup Council (to J.M.S.), a Postdoctoral Fellowship from Mitacs (to M.C.) and a Feodor Lynen Fellowship from the Alexander von Humboldt Stiftung (to I.G.).

## Acknowledgements

We are indebted to Chris Budd (University of Bath) for his helpful comments on an earlier version of this manuscript.

- Received July 26, 2015.
- Accepted August 28, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.