A two-scale Stefan problem arising in a model for tree sap exudation
The study of tree sap exudation, in which a (leafless) tree generates
elevated stem pressure in response to repeated daily freeze–thaw
cycles, gives rise to an interesting multiscale problem involving
heat and multiphase liquid/gas transport. The pressure generation
mechanism is a cellular-level process that is governed by
differential equations for sap transport through porous cell
membranes, phase change, heat transport, and generation of osmotic
pressure. By assuming a periodic cellular structure based on an
appropriate reference cell, we derive an homogenized heat equation
governing the global temperature on the scale of the tree stem, with
all the remaining physics relegated to equations defined on the
reference cell. We derive a corresponding strong formulation of the
limit problem and use it to design an efficient numerical solution
algorithm. Numerical simulations are then performed to validate the
results and draw conclusions regarding the phenomenon of sap
exudation, which is of great importance in trees such as sugar maple
and a few other related species. The particular form of our
homogenized temperature equation is obtained using periodic
homogenization techniques with two-scale convergence, which we
investigate theoretically in the context of a simpler two-phase
Stefan-type problem corresponding to a periodic array of melting
cylindrical ice bars with a constant thermal diffusion coefficient.
For this reduced model, we prove results on existence, uniqueness
and convergence of the two-scale limit solution in the weak form,
clearly identifying the missing pieces required to extend the proofs
to the fully nonlinear sap exudation model. Numerical simulations
of the reduced equations are then compared with results from the
complete sap exudation model.periodic homogenization; two-scale convergence; Stefan problem;
multiphase flow; phase change.
MSC(2010): 35B27, 35R37, 76T30, 80A22, 92C80.
This paper is motivated by the study of sap flow in sugar maple trees that are subject to repeated cycles of thawing and freezing during the sap harvest season in late winter (Ceseri & Stockie, 2013). We seek insight into the phenomenon of sap exudation, which refers to the generation of elevated sap pressure within the maple stem when the tree is in a leafless state and no transpiration occurs to drive the sap flow. Our work is based on the model derived in Ceseri & Stockie (2013) that captures the physical processes at the microscale (i.e., at the level of individual wood cells) and includes multiphase flow of ice/water/gas, heat transport, porous flow through cell walls, and osmosis. There is an inherent repeating structure in sapwood at the cellular scale that lends itself naturally to the use of homogenization ideas that we exploited in Graf et al. (2015) to obtain a multiscale model for the macroscale temperature that is coupled to a corresponding system of equations governing the microscale cellular processes. Our main objective in this paper is to provide a more rigorous theoretical justification for this multiscale model by working through the details of the homogenization process and proving results regarding existence, uniqueness and two-scale convergence.
Multiscale problems such as the one just described are characterized by geometric, material or other features that exhibit variations on widely differing spatial scales. Many mathematical and numerical methods have been developed to capture such scale separation as well as the interactions between physical phenomena operating on disparate scales (Engquist et al., 2005; Hornung, 1997). For problems having a periodic microstructure, a mathematical technique that has proven to be very effective is known as periodic homogenization (Cioranescu & Donato, 1999), and more specifically the method of two-scale convergence (Allaire, 1992; Nguetseng, 1989), which has also been extended to capture non-periodically evolving microstructures (Peter, 2007a, b; Peter & Böhm, 2009). We are interested here in applying two-scale convergence to analyze solutions of a Stefan-type problem that governs the dynamics of the ice/water interface within individual tree cells. Locally, temperature obeys the heat equation and is coupled with a Stefan condition that governs solid–liquid phase transitions at the interface. Many different approaches have been developed to analyze such phase transitions, which are well-described in Visintin (1996). With the exception of a few studies of (single-phase) water and solute transport in plant tissues (Chavarría-Krauser & Ptashnyk, 2010, 2013), periodic homogenization techniques have not been applied in the context of heat or sap flow in trees.
The approach we employ in this paper has the advantage that it applies homogenization techniques in a straightforward manner in order to obtain an uncomplicated limit model, the simplicity of which ensures that numerical simulations are relatively easy to perform. In particular, we define a reference cell that is divided into two sub-regions: , where the temperature diffuses rapidly; and , on which we define a second temperature field that diffuses slowly. Refer to Arbogast et al. (1990) and Peter & Böhm (2008) for similar homogenization approaches involving slow and fast transport. One particular challenge arising in the study of Stefan problems is that the diffusion coefficient depends on the underlying phases, so that heat diffuses differently in water or ice. Consequently, the diffusion coefficient depends on temperature (or equivalently on enthalpy) so that the governing differential equation is only quasi-linear.
Rather than attempting to analyze the sap exudation problem in its full complexity, we find it more convenient to develop our homogenization results in the context of a simpler “reduced model” defined on a similarly fine-structured domain wherein the cell-level processes are governed by a Stefan problem that involves only heat transport and ice/water phase change. In particular, we consider a domain consisting of a periodic array of cylindrical ice inclusions immersed in water. To handle the multiplicity of the ice bars, we apply the technique of periodic homogenization with two-scale convergence established in Allaire (1992) and Nguetseng (1989). Several authors have previously applied homogenization to Stefan problems, such as Bossavit & Damlamian (1981), Damlamian (1981) and Visintin (2007), where the phase change boundary is handled by separately homogenizing an auxiliary problem. In Eck (2004) on the other hand, an additional function is introduced for an aggregate state that diffuses on a slow time scale and with which all microscopic phase changes are properly captured. When we show existence for the heat equation with phase transitions, we deduce a general existence result for quasi-linear parabolic differential equations having a non-monotone nonlinearity in the diffusion operator, which is of general interest in the context of heat transport and Stefan problems (even in a single-scale setting).
This paper is organized as follows. We begin in Section A two-scale Stefan problem arising in a model for tree sap exudation by providing background material on the physics of maple sap exudation, along with a description of the governing equations at the cellular level. A reduced model involving only melting of ice is introduced in Section A two-scale Stefan problem arising in a model for tree sap exudation for the purposes of more easily deriving the two-scale convergence results. The main analytical results on existence, a priori estimates, two-scale convergence and uniqueness are presented in Section A two-scale Stefan problem arising in a model for tree sap exudation, and detailed proofs of the key results are relegated to the Appendix. Following that, we state in Section A two-scale Stefan problem arising in a model for tree sap exudation the strong form of the limit problem for the reduced model, which in turn suggests a corresponding strong form of the original sap exudation model in Section A two-scale Stefan problem arising in a model for tree sap exudation. These limit problems lead naturally to a multiscale numerical algorithm that is described in Section A two-scale Stefan problem arising in a model for tree sap exudation, after which numerical simulations of both problems are presented and compared.
Before presenting the details of the mathematical model, it is necessary to introduce some background material on the phenomenon of sap exudation. Sugar maple trees (along with a few related species such as red or black maple, black walnut, and birch) have a unique ability compared to other deciduous tree species in that they exude large quantities of sap during the winter when they are in a leafless state. Sap exudation originates from an elevated pressure in the tree stem that is generated over a period of several days during which the air temperature oscillates above and below the freezing point. The ability of maple to exude sap has intrigued tree physiologists for over a century, and various physical and biological processes have been proposed to explain this behaviour (Johnson et al., 1987; Milburn & Kallarackal, 1991; Tyree, 1995). Until recently, a significant degree of controversy existed over the root causes of sap exudation, and the most plausible and widely-accepted explanation has been a freeze–thaw hypothesis proposed by Milburn & O’Malley (1984). This hypothesis forms the basis of the mathematical model for the cellular processes underlying exudation during a thawing event that was developed by Ceseri & Stockie (2013), which was subsequently extended to capture a complete freeze–thaw cycle by Graf et al. (2015).
The Milburn–O’Malley process depends crucially on the distinctive microstructure of sapwood (or xylem) in sugar maple trees (Acer saccharum). Wood in most deciduous tree species consists of roughly cylindrical cells that are on the order of in length. These cells can be classified into two main types: vessels having an average radius of μ, which are surrounded by the much more numerous (libriform) fibers with a radius of approximately 3–4 μ. The repeating structure of vessels and fibers is illustrated in Figure 1a. The vessels have a significantly larger diameter and therefore comprise the main route for sap transport between roots to leaves during the growing season, whereas the fibers are understood to play a largely passive and more structural role. Under normal conditions the vessels are filled with sap, which is composed primarily of water but also contains as much as 2–5% sugar by weight in species like Acer. On the other hand, the fibers are thought to be primarily filled with gas (i.e., air). We note that experiments exhibit small but measurable amounts of gas also being present within the vessel sap, either as bubbles or in dissolved form.
|(a) Sapwood microstructure||(b) A single fiber/vessel pair|
Milburn and O’Malley hypothesized that during late winter when daily high temperatures peak above the freezing point, and just as evening temperatures begin to drop below zero, sap is drawn through tiny pores in the fiber/vessel walls by capillary and adsorption forces into the gas-filled fibers where it forms ice crystals on the inner surface of the fiber wall (top “cooling sequence” in Figure 2). As temperatures drop further, the ice layer grows and the gas trapped inside the fiber is compressed, forming a pressure reservoir of sorts. When temperatures rise above freezing again the next day, the process reverses, with the ice layer melting and the pressurized gas driving liquid melt-water back into the vessel where it then (re-)pressurizes the vessel compartment (bottom “warming sequence” in Figure 2). Milburn and O’Malley also stressed the importance of osmotic pressure in terms of maintaining the high stem pressures actually observed in sugar maple trees. This essential role of osmosis has since been verified experimentally by Cirelli et al. (2008) who confirmed the existence of osmotic pressure arising from a selectively permeable membrane within the fiber/vessel wall. They showed that the cell wall permits water to pass but prevents larger sugar molecules contained in the vessel sap from entering the fiber, thereby introducing a significant osmotic pressure difference between the sugar-rich vessel sap and the pure water contained in the fiber.
There are two additional physical effects not explicitly addressed by Milburn & O’Malley (1984) that are essential in order to obtain physically-consistent results for the sap thawing process. First of all, Ceseri & Stockie (2013) demonstrated the necessity for including gas bubbles suspended within the vessel sap that permit an exchange of pressure between the vessel and fiber compartments, which would otherwise not be possible owing to the incompressibility of water. Secondly, despite the pervading belief that there is no significant root pressure in maple during winter (Kramer & Boyer, 1995; Wilmot, 2011), we found to the contrary that including uptake of root water during the freezing process is absolutely essential in order that pressure can accumulate over multiple freeze–thaw cycles (Graf et al., 2015). Indeed, the need to include root pressure is confirmed by recent experiments (Brown, 2015) that demonstrate the existence of root pressure in maple trees during the sap harvest season.
The modified Milburn–O’Malley description just presented (with the exception of root pressure) was employed by Ceseri & Stockie (2013) and Graf et al. (2015) to derive a mathematical model for cell-level processes governing sap exudation during a thawing cycle. In this study, we study the same problem, including the effect of the gas phase in both cell chambers (fiber and vessel), but we will assume for the sake of simplicity that the effects of gas dissolution and nucleation are negligible. This is the primary difference between our microscale model and that in Ceseri & Stockie (2013) and Graf et al. (2015), on which it is based. Neglecting root pressure is a reasonable simplification because we are only interested here in studying a single thawing event and not capturing repeated freeze–thaw cycles.
With the above assumptions in mind, we approximate the sapwood as a periodic array of square reference cells pictured in Figure 3a. Each reference cell contains a circular fiber of radius located at the centre, surrounded by a vessel compartment that makes up the remainder of the cell. Because the vessels have considerably larger diameter, we assume that on the scale of a fiber the cylindrical geometry of the vessel can be neglected as long as we ensure that appropriate conservation principles (for mass and energy) are maintained within the vessel. This choice of reference cell is obviously a mathematical idealization that may influence fine details of vessel transport on the microscale but ultimately has minimal impact on the homogenized solution.
The fiber compartment is sub-divided into nested annular regions containing gas, ice and liquid, and the outer radii of the phase interfaces are denoted (for gas/ice) and (for ice/water). The vessel contains a circular gas bubble of radius which has no specified location but rather is included simply to track the amount of gas for mass-conservation purposes. One additional variable is introduced to measure the total volume of water transferred from fiber to vessel. The region lying outside the fiber and inside the boundary of represents the sugary sap-filled vessel. Note that during a thawing cycle, we are only concerned with a vessel containing liquid sap (no ice) because of the effect of freezing point depression, which ensures that any given vessel thaws before the adjacent fiber(s). This reference cell geometry should be contrasted with that depicted in Ceseri & Stockie (2013, Fig. 3.1).
|(a) Reference cell, with ice layer||(b) Reference cell, completely melted|
For the moment, we will consider the four solution variables , , , as depending on time only, with an additional dependence of temperature on the microscale spatial variable; however, beginning in Section A two-scale Stefan problem arising in a model for tree sap exudation when we derive macroscale equations for the homogenized problem, these variables will also depend on the global spatial variable that denotes the location of the reference cell within the tree stem. Within a reference cell, the dynamics for , , and are governed by four differential equations whose derivation can be found in Ceseri & Stockie (2013). The first is the Stefan condition for the ice/water interface in the fiber
|where denotes the microscale temperature variable that depends on both time and the local spatial variable (which needs to be distinguished from the macroscale temperature variable introduced later) and is the normal derivative at the interface. Here, the enthalpies of water and ice ( and , resp.) are evaluated at the freezing point, ; consequently, the difference represents the latent heat of fusion. We describe heat transport using a mixed temperature–enthalpy formulation, in which the thermal diffusion coefficient is written as a function of enthalpy . Following Visintin (1996), we take to have the piecewise affine linear form|
|where , are the densities of water and ice respectively, and , are the thermal conductivities. Note that in this temperature–enthalpy formulation has units of and is referred to as a thermal diffusion coefficient, to distinguish it from the more usual “thermal diffusivity” (which is defined as the ratio and has units of ). The governing equations for and are discussed later in Sections A two-scale Stefan problem arising in a model for tree sap exudation–A two-scale Stefan problem arising in a model for tree sap exudation as a result of the two-scale convergence analysis and are the solutions of the system (4.20a–e). Note that the final term in the Stefan condition (2.1a) was neglected in Ceseri & Stockie (2013) and serves to capture the effect on the phase interface of fiber–water volume changes due to porous flow through the fiber/vessel wall.|
|The next two differential equations embody conservation of mass in the fiber|
|and the vessel|
|Note that within the sapwood there are many more fibers than vessels (as depicted in Figure 1a), so that the effect of fiber–vessel flux terms should be increased to account for the multiplicity of fibers. With this in mind, we have multiplied appropriate fluxes by the parameter in (2.1d) that represents the average number of fibers per vessel and has a typical value of . The final differential equation describes water transport through the porous fiber/vessel wall in response to both hydraulic and osmotic pressure|
Here, we denote the pressure variable by , where superscripts / refer to fiber/vessel and subscript denotes the liquid water phase. The constant parameter is the fiber/vessel wall conductivity, is the wall surface area, is the vessel sugar concentration, and is the universal gas constant. Note that because is defined inside , we should strictly be using the microscale temperature in the osmotic term in (2.1e), but this would lead to a significant complication in any numerical algorithm due an additional nonlinear coupling between scales. Therefore, we have used instead, which is a reasonable approximation because temperature variations throughout the reference cell are small.
Several intermediate variables have been introduced into the above equations. They are determined by the following algebraic relations:
All constant parameters appearing in the above equations are listed in Table 1 along with typical values.
|Side length of reference cell|
|, Radius of|
|Area of fiber/vessel wall|
|Thickness of fiber/vessel wall|
|Number of fibers per vessel||16||–|
|Tree stem radius||0.25|
|Specific heat of water|
|Specific heat of ice|
|Enthalpy of water at|
|Enthalpy of ice at|
|Thermal conductivity of water|
|Thermal conductivity of ice|
|Density of water|
|Density of ice|
|Freezing temperature for water||273.15|
|Heat transfer coefficient|
|Molar mass of air||0.029|
|Universal gas constant||8.314|
|Gas/liquid surface tension||0.076|
|Vessel sugar concentration (2%)||58.4|
|Hydraulic conductivity of fiber/vessel wall|
There is one special case to consider, namely when a fiber initially containing an ice layer is above the freezing point for long enough time that the ice melts completely. In the moment the ice layer disappears, the reference cell geometry appears as in Figure 3b and the cell-level equations must be modified as follows. First of all, must change to account for the fact that there are two possible values of thermal diffusivity, one in the region containing the gas and another in the liquid. Furthermore, the gas/ice and ice/water interfaces merge so that Eq. (2.1a) drops out and we identify a new fiber gas/water interface as . This leads to the following simplified version of (2.1c)
The equations for the temperature and enthalpy variables appearing in the microscale model above are derived in the next section in the context of a simpler problem involving only melting ice. Despite the fact that this reduced model involves only a single microscale variable for the dynamics of the ice–water interface (in addition to the temperature), the equations for temperature and enthalpy remain the same, and we will show that the microscale model above is completed by Eqs. (4.20a)–(4.20e).
We now shift our attention to the macroscale problem, which captures the dynamics of thawing sap within a cylindrical tree stem having a circular cross-section . There is a clear separation of scales in that the tree has radius on the order of tens of centimetres whereas the cell-level processes occur over distances on the order of microns. Let represent the macroscale spatial variable and the microscale variable on the reference cell. Then, our main aim in this section is to determine equations for the temperature and enthalpy variables not only in the reference cell, and , but also on the macroscale, and .
The derivation of these equations may be simplified significantly by considering a reduced problem that involves only ice/water phase change and leaves out all other physical processes (porous flow, gas bubbles, surface tension, etc.). To this end, we consider a periodic array of melting “ice bars” as pictured in Figure 4a, situated inside a slightly more general domain having Lipschitz boundary that contains both water and ice in the form of circular inclusions. Let be a reference cell that captures the configuration of the periodic microstructure, and for which represents its actual physical size with (although we focus on dimension , the theoretical results proven here apply to any dimension). The reference cell is divided into two sub-domains and that are separated by a Lipschitz boundary as shown in Figure 4b. For simplicity, we take to be a circle of radius satisfying . The primary feature that we exploit in our homogenization approach is that within heat must diffuse rapidly, whereas in there is a relatively slow diffusion of heat.
|(a) Periodically-tiled domain||(b) Reference cell for reduced model|
We next introduce a small parameter that corresponds to the size of the periodic microstructure (and must be distinguished from the physical size because we will eventually take the limit as ). The domain may then be decomposed into three -dependent sub-domains: (which is connected), and two disconnected components consisting of the region and the boundary curves . This decomposition is illustrated in Figure 4. To avoid technical difficulties, we assume that does not touch the outer boundary of , so that and .
The major advantage of this reduced model is that the reference cell problem simplifies significantly, with the only unknowns being and temperature. We proceed with the temperature and enthalpy equations.
Throughout the analytical developments of this paper, we employ what is known as the two-phase formulation of the Stefan problem, in which the heat diffusion equation is posed in a mixed form involving both temperature and enthalpy. Assuming that material properties of water and ice remain constant, the temperature can be written as a piecewise linear function of enthalpy as follows (Visintin, 1996)
where and denote specific heats of water and ice respectively, and is the freezing point of water (parameter values are listed in Table 2). A distinguishing feature of this temperature–enthalpy relationship is that when temperature is equal to the freezing point, the enthalpy varies while temperature remains constant – this behavior derives from the fact that a certain amount of energy (called latent heat) is required to effect a change in phase from solid to liquid at the phase interface.
Because the function is neither differentiable nor invertible, we instead employ in our model a regularized version defined as
which has “rounded corners” that are smoothed over the short intervals and . Note that we have also introduced a small positive slope within the central plateau region near (refer to Figure 5). These modifications ensure that is a continuously differentiable, invertible and monotone increasing function of enthalpy. Incidentally, such a regularized function is most likely a more accurate representation of what one would actually observe in a real physical system.
We now describe the solution decomposition into slow and fast diffusing variables. Let functions and denote the fast-diffusing temperature and enthalpy components respectively, with both defined on the sub-region . Similarly, let and denote the slowly-diffusing temperature and enthalpy on . We may then state the strong formulation of the two-phase Stefan problem as
where is given by (2.1b) and is the ambient temperature imposed on the outer domain boundary. Note that these equations capture the phase interface location implicitly through the relationship and no explicit Stefan-type phase interface condition is imposed.
It is not feasible to solve the system of equations (3.2) derived in the previous section using direct numerical simulations owing to the presence of the microstruture with being very small. This section contains the primary theoretical results pertaining to an upscaling of the reduced problem which is achieved by characterising the limit as . Lemmas and theorems are stated here, and the proofs are relegated to the Appendix.
In order to make the problem ameanable to periodic-homogenization techniques, we begin by transforming Eqs. (3.2) into a weak formulation. To avoid technical difficulties, we replace the Robin boundary condition (3.2c) by a Dirichlet condition (i.e., we consider the formal limit ) and refer the reader to Graf & Peter (2014) for a detailed discussion. This requires first defining some appropriate solution spaces:
where the “primes” denote dual spaces and represents the time interval of interest for some fixed . The corresponding test spaces are , and . We also need to introduce notation for inner products, with representing the -inner product with respect to space of two functions in for , whereas denotes that an additional time integration is performed over the interval with . Finally, we let denote the dual pairing on . Later, we will show that , so that we can interpret as , where represents the Riemann curvature tensor.
We are now prepared to state the weak form of the heat-diffusion problem. Assuming that initial values and are smooth, non-negative and bounded functions, and that a Dirichlet condition is imposed at the outer boundary , our goal is to find such that
for all . Note that represents the outward-pointing unit normal vector on , and that temperature and enthalpy are connected via , or equivalently . We assume that is positive, bounded and smooth such it can be extended to in . We note again in closing that slow diffusion is induced in the problem via the factor multiplying terms in Eqs. (4.1) that involve the diffusion coefficient .
In this section, we apply a procedure developed by Arbogast et al. (1990) to transform the model (4.1) by combining enthalpies and into a single function defined on the whole -independent domain . We use the fact that is the only boundary of to obtain
for all . After substituting this expression into (4.1a) we obtain
for all . Hence, Eqs. (4.1) have been replaced with
for all . We then define the function by
so that with conditions (3.2b) and (3.2e) the function is guaranteed to be continuous and weakly differentiable. Furthermore, we define where for are indicator functions for and respectively. Solving (4.1) is then equivalent to finding such that
for all , where we have used that and .
We perform one further transformation of (4.4) that makes the Dirichlet boundary condition homogeneous. To this end, we define where is extended continuously to . Then (4.4) is equivalent to finding such that
for all .
To prove the existence of a solution to (4.4) for every , we formulate a theorem of existence, which is strongly inspired by a proof for a related result found in a set of unpublished lecture notes by Wolff (2016). In Theorem 4.1 we introduce a general existence result for parabolic equations with non-monotone non-linearities in the diffusion operator. The proof is based on the Rothe method. We state the theorem next and provide the proof in Appendix A two-scale Stefan problem arising in a model for tree sap exudation.
Consider the initial–boundary value problem
Let the conditions (4.7) be satisfied. Then
define a family of operators and an element for which the following hold: