Freeze dehydration vs. supercooling in tree stems: physical and physiological modelling

Frost resistance is the major factor affecting the distribution of plant species at high latitude and elevation. The main effects of freeze-thaw cycles are damage to living cells and formation of gas embolism in tree xylem vessels. Lethal intracellular freezing can be prevented in living cells by two mechanisms: dehydration and deep supercooling. We developed a multiphysics numerical model coupling water flow, heat transfer, and phase change, considering different cell types in plant tissues, to study the dynamics and extent of cell dehydration, xylem pressure changes, and stem diameter changes in response to freezing and thawing. Results were validated using experimental data for stem diameter changes of walnut trees. The effect of cell mechanical properties was found to be negligible as long as the intracellular tension developed during dehydration was sufficiently low compared to the ice induced cryostatic suction. The model was finally used to explore the coupled effects of relevant physiological parameters (initial water and sugar content) and environmental conditions (air temperature variations) on the dynamics and extent of dehydration. It revealed configurations where cell dehydration could be sufficient to protect cells from intracellular freezing, and situations where supercooling was necessary. This model, freely available with this paper, could easily be extended to explore different anatomical structures, different species and more complex physical processes.


23
Plant species distribution at high latitude and elevation is mainly driven by their resistance to freezing stress 24 (Charrier et al., 2013a). Freeze-thaw cycles are known to damage living cells (Siminovitch and Scarth, 1938; water or ice depending on local temperature, and gas. This gas compresses or expands in response to water 147 flows entering or leaving the vessels, thus creating pressure variation according to the ideal gas law, as done in 148 Ceseri and Stockie (2013); Graf et al. (2015). 149 Heat transfer and phase changes are calculated at the tissue scale and driven by external temperature 150 variations. Temperature (and phase for xylem sap) is computed at the cellular scale from the tissue-scale values. 151 Vessel sugar content impacts tissue-scale phase change through freezing point depression (FPD). 152 Water fluxes occur between the different elements, see the blue arrows in figure 1c. These water fluxes are 153 determined by the differences in water potential (hydrostatic/turgor, osmotic, cryostatic) across cell membranes. 154 For each elastic compartment (VACs, bark cells), the balance of water fluxes results in volume changes, which 155 are then used to calculate changes in tissue diameter, as well as changes in turgor and osmotic potential. The anatomical description used in the model is shown in figure 1c. Vessels are arranged regularly along the 159 ray, with the vessel number, N vessels , computed using a linear vessel density and the size of the xylem tissue: where lvd, R xylem and R pith are the linear vessel density, the xylem and the pith radius, respectively. Each 161 vessel has a given number of VACs associated with it, N vac , calculated using a ratio of the vessel-VAC exchange 162 area to the projected VAC area: where A vac−v , R vac and l vac are the vessel-VAC exchange area, the VAC radius and VAC length, respectively.

164
The number of parenchyma rays is computed using a tangential ray density and the branch diameter: where trd is the tangential ray density. The number of bark cells connected to the parenchyma rays is with V 0 w,bark the inital volume of water in the bark accesible from the rays, equal to BW F × V 0 bark , with BW F 167 the bark water fraction accessible from the rays, and V 0 bark the initial bark volume. V 0 w,Bark cell is the initial 168 bark cell water volume.

169
Heat transfer and phase change 170 Heat transfer and phase change are calculated at the tissue scale through the heat equation in a 1D axi-symmetric 171 model in cylindrical coordinates: where H is the enthalpy, T the temperature, k th the thermal conductivity and ρ the density. This equation is temperature at the tissue scale is computed locally based on C v s , the local vessel sugar content: with T c = 273.15K. Note that Eq. (7) is only valid for water. Heat equation, Eq. (5), is solved on a grid finer 181 than that formed by the regular arrangement of vessels along the ray, but with grid points that correspond to 182 the vessel locations, so that temperature and enthalpy values do not need to be interpolated to be obtained at 183 the cell scale. The cell size of the thermal grid, ∆ r , is small enough to reach spatial convergence of the results.

184
The vessel sugar concentration is linearly extrapolated from the vessel locations to the heat equation grid so 185 that Eq. (7) can be applied at any thermal grid point. A vessel ice fraction, δ iv , is computed for each vessel as Water fluxes between elements are computed using Darcy's law. Along the parenchyma ray: where k ray is the ray porosity, ∆l v is the distance along the ray between two vessels, µ is the dynamic viscosity where k vac is the vessel-VAC wall porosity, W the vessel-VAC wall thickness, and p v w the vessel water pressure.

196
Q vac−v is positive for water fluxes going towards the vessels. The cryo-suction pressure induced by vessel freezing 197 is computed at each vessel location as (Loch, 1978;Beck et al., 1984) 198 where ρ w and L are the water density, and water latent heat of fusion, respectively.

199
Eq. (9) implies that cryo-suction will draw water in a vessel from its VACs once this vessel is frozen. It is 200 assumed that this water turns into ice when entering a vessel. We furthermore assume that, if a vessel is frozen, 201 i.e., if δ iv = 1, Q vac−v cannot be negative.

202
Pressure-volume relationships in living cells 203 In living cells, the balance of water fluxes results in volume changes. For the VACs: where the second term on the right hand side represents the balance of fluxes entering/leaving the VAC from/to 205 the ray. Between the xylem tisue and the bark, the total water flux is where Q R xylem ray is the water flux computed between one bark cell and the VAC closest to the bark. This water 207 flux is rescaled by the number of bark cells to obtain the volume change at the bark cell scale: These living cells volume changes are related to turgor pressure variation through (Steudle et al., 1977) 209 dp vac for the VACs, and 210 dp bark cell for simplicity and for the validity of Eq. (16), the elastic modulus is assumed to be constant, even though 221 experimental results suggest it could vary with turgor pressure and/or cell volume (Steudle et al., 1977).

222
Volume changes also result in osmotic pressure changes through changes in sugar concentration: and for the bark cells, with T vac m , the freezing temperature computed for each VAC, and T bark cell m , the freezing temperature for the 229 bark cells. Note again that these last two equations are only valid for water.

230
Pressure-volume relationships in vessels 231 Vessels contain gas that compresses or expands depending on water fluxes leaving or entering vessels. Following , and as shown in figure 2d, we assume that this gas is contained 233 in one cylindrical bubble located at the center of each vessel. Applying flow rate conservation between the 234 gas/water (or gas/ice) interface and the vessel/VAC wall, the bubble radius, r v g , varies as where L z is a vertical dimension that is introduced for unit consistency but that has no influence on model 236 results. Changes in r v g induce changes in gas pressure through the ideal gas law: In previous equation, n v g is the gas quantity inside the gas bubble and R g the ideal gas constant. Finally, the 238 pressure in the liquid water/ice phase is obtained using Laplace equation: where σ gw is the liquid water/gas interface surface tension. Note that we assume that no change occurs in Eq.

240
(22) when the vessel is entirely frozen. Similarly to living cells, volume changes in vessels also result in sugar 241 concentration changes: where n v s and R v are the vessel sugar quantity, taken as a constant, and vessel radius, respectively. These 243 changes in vessel sugar concentration impact phase change at the tissue scale through Eq. (7).

244
Diameter changes 245 Finally, diameter changes are obtained from living cell volume changes, Eqs. (11) and (13). The total volume 246 of water in VACs is computed at each instant: with the volume variation equals to The xylem diameter, considered as a cylinder, is computed as At each instant, the volume of the bark tissue is equal to the initial volume minus the volume lost by dehydration: from which the stem diameter, considered as a cylinder, can be computed: VACs and bark cells.

280
In the present section we present the results obtained using the model described previously. Unless stated 281 otherwise, all parameter values are presented in table 2.

282
Freeze-thaw cycle 283 In the present section we introduce the model result for a 24 hours long freeze-thaw cycle (base case 1 in table   284 2). The case simulated corresponds to a walnut stem with a diameter of 2 cm, initially at +5 • C and subjected to variation, and dehydration is fast enough to limit intracellular freezing probability, i.e., melting temperature 305 based on living cell sugar concentration is almost always lower than the external temperature.

306
Cell dehydration from freezing also induces tissue diameter shrinkage, see figure 3e. Once temperature increases, cryostatic suction starts decreasing (incrasing ice potential) but as we hypo-315 thetized that ice does not flow out of vessels, the osmotic potential diverges from the decrease in cryosuction. It 316 does rather slightly increase with temperature, as vessel water pressure does, in relation to temperature effect 317 in ideal gas law. As the dehydration of living cells ceased, tissue diameter and FPD remained stable until the 318 temperature reaches the vessel melting point.

319
When thawing starts, the difference in osmotic potential between living cells and vessels is so high that the 320 swelling of the diameter is extremely rapid. As no irreversible process takes place, e.g., cell mortality due to 321 intracellular freezing, diameter and pressure changes are fully reversible.  11.8MPa. On the other hand, using Eq. (10) for T = −10 • C, the cryostatic suction is p v ice = −12.4MPa. Thus, 406 a sufficient condition for protection against intracelluler freezing becomes where Π vac is the living cell osmotic potential. At equilibrium, apoplast and symplast potentials are equal.

408
Using Eq. (9) and assuming the extracellular ice has zero osmotic potential leads to Hence combining the last two equations, we observe that negative p t (wall tension generation) and positive p v w 410 can induce Π vac < −p v ice . Note that, as shown with the previous example for T = −10 • C, Eq. (29) is often 411 largely sufficient to protect against intracellular freezing, but Eq. (30) is useful to evidence some scenarios that 412 can lead to non-equilibrium freezing, i.e., to the need of supercooling for protection against intracellular freezing.

413
We thus find the following scenarios: (1) if equilibrium is not attained, i.e., if water flow is too slow or impossible,

415
The model that we developed and validated has been used to explore the effect of environmental, physiological 416 and mechanical parameters that could favour one or several of these different scenarios.

417
The first scenario, i.e., transient supercooling, is favoured by parameters like cooling rate and initial sugar 418 concentration as seen in figure 6. We limited our study to cooling rates up to 8K/h, which is already high 419 for field conditions. Increasing cooling rate leads to an increased supercooling degree, which is coherent with is likely to be a resistance mechanism mainly through its effect on nucleation temperature of supercooled cell 429 water rather than through its effect on freeze-induced dehydration.

430
The second scenario, i.e., steady supercooling due to the generation of cell wall tension has been explored 431 in figure 5. We demonstrated that cell wall elasticity has no effect on dehydration dynamic, xylem pressure

460
Lowering the minimal temperature systematically increases the maximal supercooling degree ( figure 7).

461
Whether it is transient or steady depends on cell wall tension or on the water content, as explained previously.

462
But even when it is transient, at the lowest tempererature (−15 • C) the viscosity is so high, i.e., more than 463 1000 Pa.s, that the time needed to reach equilibrium can be up to more than 60 hours! In practice, for 464 temperatures lower than −10 • C dehydration is never sufficient to protect cells against intracellular freezing.

465
For such temperature, at the lowest water content, the predicted living cell sugar concentration exceeds 8000 466 mol/m 3 , which is higher than the supersaturation limit (7152 mol/m 3 , Hartel and Shastry (1991)). Thus, at such

489
We hypothesize that ice appears in previously gas filled compartments due to the freezing of a thin water 490 layer covering the walls of these compartments. The resulting ice then draws additional water by cryosuction.

491
This thin water layer could appear due to the condensation of water vapour under decreasing temperature.

492
Water vapour may also condensate at a decreased humidity level in nanopores (Kelvin effect). Ice propagation 493 may also occurs from vessels to intercellular spaces (Ball et al., 2002). We note also that, following Konrad i.e., as C s = n s /V , the same concentration is obtained for a higher V when n s is higher. Decreasing bark water 511 content leads to reduced stem diameter shrinkage and pressure ( figure 7). We have explained in the previous 512 section that decreasing bark water content indeed lead to a better protection by freeze induced dehydration, 513 while increasing sugar content has an effect that depends on water content. The model therefore predicts that 514 non-acclimated trees have a much higher shrinkage than acclimated trees.
Water/Ice Gas Figure 2: a) Enthalpy-Temperature relationship: the regularized relationship spreads latent heat release/storage over a very small temperature interval of thickness ∆T m = L/c ∞ . In the ice or liquid state the enthalpy is a linear function of temperature, i.e., H = c i T or H = c w (T − (T m + ∆T m /2)) + L, with c i and c w the ice and liquid water specific heat capacities, respectively. b) Changes of physical properties and vessel ice fraction with local enthalpy value. c) Changes of cell elastic modulus with cell volume, from B t the value for a turgid cell, to B f the value for a flacid cell. Around V collapse , the elastic modulus is regularized in the same way as the H-T relationship around the melting temperature. d) Geometry of a vessel of radius R v , wall thickness W/2, containing a gas bubble of radius r v g . Figure 3: Model results for one freeze-thaw cycle. The stem (2 cm in radius), initially at +5 • C, is subjected to a sinusoidal external temperature variation, down to −10 • C, over 24 hours. The parameters used for this case correspond to base case 1 in table 2. a) FPD dynamic for two VACs and one bark cell, along with external temperature changes. b) Time course of osmotic pressure for two VACs and one bark cell, and cryostatic pressure for two vessels. c) Time course of p v w for two vessels. d) Time course of living cell turgor pressure for two VACs and one bark cell. Note that for the bark cell, only one data point every 50 is shown to enhance readability. e) Tissue diameter changes (total stem and xylem) and value given by the correlation in (Améglio et al., 2001b). f) Mean vessel pressure.