A thermodynamically consistent quasi-double-porosity thermo-hydro-mechanical model for cell dehydration of plant tissues at subzero temperatures

Abstract

Many plant tissues exhibit the property of frost resistance. This is mainly due to two factors: one is related to metabolic effects, while the other stems from structural properties of plants leading to dehydration of their cells. The present contribution aims at assessing the impact of ice formation on frost-resistant plant tissues with a focus on structural properties specifically applied to Equisetum hyemale. In this particular case, there is an extracellular ice formation in so-called vallecular canals and the pith cavity, which leads to a dehydration of the tissue cells to avoid intracellular ice formation, what would be fatal for the cells and subsequently for the whole plant. To address the underlying phenomena in the plant, a coupled thermo-hydro-mechanical model based on the Theory of Porous Media is introduced as the modelling framework. The dehydration of the tissue cells is referred to as of quasi-double-porosity nature, since the water is mobile within the intercellular space, but confined to the cells in the intracellular space and consequently kinematically coupled to them. However, the mass exchange of water across the cell wall is considered. The presented numerical example shows the strong coupling of the underlying processes as well as the quasi-double-porosity feature. Finally, it supports the experimental finding of the vallecular canals as the main location of ice formation.

1 Introduction

Frost hardiness is an evolutionary adaption of terrestrial plants to climates with frost occurrence. Two main strategies are essential. Whereas the (deep) supercooling process alters ice crystal growth and form [20], dehydration prevents ice formation within living cells (freezing avoidance [34]). In the latter case, extracellular ice formation occurs in specific intercellular spaces or extracellular on the plant surface, and water is drawn from adjacent cells by the lower water potential of the ice which accumulates continuously. Cell dehydration increases cellular solute concentration, thereby lowering the freezing point further [6, 38, 40].

Freezing avoidance is of interest in plant research for more than 150 years [50]. So far, extracellular freezing was detected in many different species including spore-bearing vascular plants and mosses [40, 44, 61] and different plant parts as petioles [43, 50], leaves [5], buds [36] and woody tissue [2, 60].

In the present contribution, the process of extracellular ice formation is studied numerically with regard to Equisetum hyemale var. robustum. Equisetum hyemale (winter scouring rush) is a spore-bearing plant consisting of nodes and internodes which prefers moist habitats and shows an extensive internal canal system [64]. Extracellular freezing takes place within the pith cavity (PC), the vallecular canals (VC), the intercellular spaces of the chlorenchyma and the substomatal chambers as described by [61]. Therein, the analysis of the freezing sites is presented and taken into account for the current model and the numerical investigation together with various other studies regarding the extracellular freezing process within Equisetum hyemale by [46, 58] and [38].

So far, a holistic treatment of biomechanical phenomena as multiphasic porous biomaterial in terms of mechanical, genetic, metabolic, chemical and potentially other influences, is already established in the field of human biomechanics, compare, for example, [28] and quotations therein. In plant biomechanics, poroelastic approaches have been used occasionally, compare [4, 54] or [56, 57], where plants are treated as a multiphasic material of solid and fluid based on a simple binary poroelastic model according to Biot’s theory [8]. For more sophisticated models with multiple constituents, the first-principle-based framework of the Theory of Porous Media (TPM) is particularly suited to derive thermodynamically consistent constitutive models, while it also accounts for classical biological approaches. Regarding the fluid dynamics of macro(bulk)-flow in plants, compare [12] or [45], as well as micro-porosity-based flow such as cell dehydration, compare [16] or [39], one usually proceeds from a water potential as the flow-driving force. In the present contribution, this somehow fundamental idea is combined with the continuum-thermodynamical TPM, compare [9, 19, 23, 24] and citations therein, which allows to couple the mechanical, hydraulic and thermal effects of plants in the framework of finite deformations.

Generally, when frost-resistant plants are cooled down, the temperature in the tissue is decreasing faster than in air/vapour-filled compartments due to the higher thermal conductivity of the tissue. This leads to a condensation of water in these compartments at the surface from the tissue. Note that the condensed water freezes at $273.15\,$K under atmospheric pressure starting at the surface from the biggest compartment to smaller ones. The ice formation leads to a so-called cryo-suction of liquid water to the freezing site, when the temperature is decreased further [17]. Since the focus is on the strategy of dehydration as key mechanism of frost hardiness, the modelling approach is chosen accordingly. Regarding the impact of ice formation on plant tissues, distinction is made between localised and rather dispersed ice. From the modelling point of view, dispersed ice formation within the plant tissue, as for example in Betula nana (dwarf birch) [60], requires the explicit inclusion of ice as a model constituent. As a result, phase transformation processes between liquid water and solid ice have to be included into the overall approach. This results in a quaternary TPM model with a solid skeleton as load-bearing structure, water in liquid and frozen states, and air, as emphasised in [30, 31]. However, for a local accumulation of larger ice bodies within extracellular spaces, as for example predominantly in Equisetum hyemale [61], the impact of ice formation on the plant tissue can be modelled by a reduced ternary TPM model with the constituents solid skeleton, liquid water and air, as pointed out in [32]. In this particular case, the main locations of ice formation are in the pith cavity (open space in the centre of the stem) as well as in the vallecular canals (open spaces in the wall of the stem). These compartments are not part of the porous medium. Note in passing that ice formation may also occur at other locations within the porous medium, but is there of minor importance. Following this, ice build-up appears at the external side of the body’s boundary causing a drop in water potential [46], which attracts even more water to the freezing site. Imposing Dirichlet pressure boundary conditions at these surfaces, where ice formation can occur, exhibits ice build-up by the ice-induced drop of the water potential. This approach reduces the computational effort dramatically compared to a quaternary model that additionally has to numerically account for phase transformations. In this article, the case of localised ice in Equisetum hyemale is addressed.

On a conceptual basis, the distinction of ice-formation pattern and the corresponding modelling approaches represent a novelty with respect to an understanding of the underlying thermo-hydro-mechanical (THM) processes and the influence of the microstructure upon ice formation in plant tissues. In particular, the TPM approach allows for a coupled assessment of ice build-up and the deformation of Equisetum hyemale, which is in good agreement with experimental investigations [61]. The location of maximum ice accumulation in Equisetum hyemale has been determined and is also confirmed by the same experimental investigations [61]. Generally, the key is the understanding of the availability of water with focus on a microstructurally based process of plants, which is essential for their frost hardiness. This is the dehydration of their cells occurring via water percolation through pores in the cell walls.

The importance of the porosity at the cell scale has been emphasised by [39, 66] and many others. This feature indicates the existence of pore spaces at two scales, at the microscale within the cells and at the macroscale in terms of the vascular system that transports water and nutrients at the scale of the whole plant. Within porous-media theories, a multiscale porosity feature has been set into a sound continuum-mechanical framework by [11, 15, 33] as well as [70]. These articles proceed from the assumption of individual fluid constituents at two porosity scales and, consequently, individual states of motion for the fluids in the micro- and macro-pore spaces with appropriate mutual interactions. This is usually referred to as been of double-porosity nature. In contrast, the present approach does not require an independent state of motion for water at two porosity scales, since the water is either confined to the solid skeleton as cell water in the micro-pore space or it is mobile within the macro-pore space. Hence, there is no independent state of motion for the water at the microscale, although the exchange of water of the cells with the macro-pore space is considered via cell dehydration. This introduces the so-called quasi-double-porosity feature. Therefore, beyond establishing the framework for the description of the involved THM processes in plant tissues, this article specifically aims at assessing the water management of plant tissues in a freezing environment with a focus on the availability of water at two porosity scales in Equisetum hyemale.

In the following sections, a macroscopic ternary model of solid skeleton, liquid water and air is utilised for the description of stem-based plant tissues like that of Equisetum hyemale. Basically, triphasic materials of solid, liquid and gas need to be described by their mass, momentum and energy balances of all constituents resulting in a coupled set of fifteen scalar equations as far as fully three-dimensional (3-d) problems with different constituent temperatures have to be considered. While the solid mass balance can usually be integrated to yield a relation between the initial and the current partial solid density, this does not apply in case of the existence of mass production terms. In this case, the integration of the mass production results in an implicit equation that has to be solved numerically. As a result, the balance equations are governed by primary variables given by the solid density or, as in the present contribution, by the solid volume fraction and constituent-individual fluid pressures for the mass balances, the solid displacement and the fluid or seepage velocities for the momentum balances and by constituent-individual temperatures for the energy balances. Beyond the primary variables, a lot of secondary variables have to be determined as functions of the primary variables, such that the set of equations can be solved. In this context, it should be noted that except for some very specific simple problems, there is no chance to obtain exact solutions from the so-called strong formulation of the initial-boundary-value problem (IBVP). Instead, weak solutions have to be preferred that can be based on the finite element method (FEM), where a weak solution ${\varvec{y}}$ has to be found for the function ${{\mathcal {G}}}({\varvec{y}})$, such that ${{\mathcal {G}}}({\varvec{y}})=\mathbf{0}$. Here, ${{\mathcal {G}}}({\varvec{y}})$ is the vector of all governing equations depending on all primary variables included in the vector ${{\varvec{y}}}$.

These general statements will be reduced to the specific model for Equisetum hyemale in the following sections.

2 Macroscopic TPM model

2.1 Theoretical modelling and basic quantities

The theoretical foundation for the description of extracellular ice formation in plants is ideally given by the TPM. Regarding Equisetum hyemale, the plant geometry is schematically depicted in Fig. 1 with its constituents and involved pore spaces. The THM model applied to the plant tissue comprises the porous material with its micro- and macro-pore spaces but excludes the domains beyond the surfaces towards the pith cavity and the vallecular canals. In detail, the cross section of Equisetum hyemale is shown in Fig. 1a in its natural state with the pith cavity in the centre and the vallecular canals within the tissue, both pointing with their longitudinal axes towards the out-of-plane direction, thus defining a sort of preferred direction. This indicates that the mechanical, hydraulic and potentially also the thermal behaviour of the plant are anisotropic. A representative elementary volume (REV) of the porous tissue material excluding the pith cavity and the vallecular canals, as shown at the top of Fig. 1b, is composed of a solid, water and air. Note again that the solid skeleton $\varphi ^{S}$ is itself a multiphasic material made of lignified elements as well as of tissue cells that contain a significant amount of water up to 90 %, while the macro-pore space is filled with two immiscible fluids, liquid water $\varphi ^{L}$ and gaseous air $\varphi ^{G}$, such that the overall aggregate $\varphi$ is given as:

$$\begin{aligned} \varphi = \underset{\alpha }{\bigcup } \, \varphi ^{\alpha } = \varphi ^{S}\,\cup \,\varphi ^{L}\,\cup \,\varphi ^G, \quad \text {where} \quad \alpha = \{S, L, G\}. \end{aligned}$$

(1)

Furthermore, the detailed view at the bottom of Fig. 1b indicates the macro-pore space (or in plant-biomechanics terminology also referred to as the intercellular or extracellular space). The remaining part of the detailed view is composing the solid skeleton, which includes the intracellular space as well as the cell wall with its micro-pores. All constituents are smeared out (volumetrically averaged) over the entire control space, as depicted in Fig. 1c, with three emerging substitute continua.

Fig. 1

In a, the cross section of Equisetum hyemale is shown in its natural state with the pith cavity (PC) and the vallecular canals (VC). At the top of b, the basic constituents are displayed schematically including the multiphasic solid skeleton together with liquid and gas in the macro-pore space. The detail at the bottom of b reveals the two pore spaces as well as the cell wall and the intracellular space. Note that the intracellular space and the cell wall including the micro-pore space constitute the solid skeleton. In c, the idealised model yields three substitute continua (solid, gas and liquid), which are smeared over the entire control space

Note that the impact of extracellular ice formation in the pith cavity and the vallecular canals is addressed by introducing pressure boundary conditions representing a drop in the water potential due to the formation of ice at the surface of the porous tissue material. As a result, the inclusion of an ice phase in the modelling approach can be omitted. Moreover, note that air is exclusively present within the intercellular (macro-pore) space, such that the double-porosity feature is considered for water only.

The local amount of the individual constituents introduced in (1) is described by their volume fraction

$$\begin{aligned} n^\alpha = \frac{{\mathrm d} v^\alpha }{{\mathrm d} v}, \end{aligned}$$

(2)

where ${\mathrm d} v^\alpha$ is the volume element of $\varphi ^\alpha$ and ${\mathrm d} v$ the volume element of the overall aggregate $\varphi$. As there is no vacant space within the plant tissue, the saturation condition holds at any instant of time yielding

$$\begin{aligned} \sum \limits _\alpha n^\alpha = n^S+n^L+n^G\,=\,1. \end{aligned}$$

(3)

The volume fraction $n^\alpha$ of a constituent $\varphi ^\alpha$ also relates the effective (intrinsic or realistic) density $\rho ^{\alpha R}$ to the partial density $\rho ^\alpha$, such that

$$\begin{aligned} \rho ^{\alpha R} = \frac{{\mathrm d} m^\alpha }{{\mathrm d} v^\alpha }, \rho ^\alpha = \frac{{\mathrm d} m^\alpha }{{\mathrm d} v} \longrightarrow \rho ^\alpha = n^\alpha \rho ^{\alpha R}. \end{aligned}$$

(4)

Therein, ${\mathrm d} m^\alpha$ is the mass element of $\varphi ^\alpha$. Note that the intrinsic density is constant under isothermal conditions for a materially incompressible constituent; however, the intrinsic density may vary due to temperature changes. Furthermore, the immiscible fluids within the macro-pore space are usually addressed by their saturation $s^\beta$, where $\beta =\{L,G\}$. The respective saturation $s^\beta$ is defined by relating the volume fraction $n^\beta$ to the porosity $n^F$ via

$$\begin{aligned} s^\beta = \frac{n^\beta }{n^F} \text {with} n^F = n^L+n^G. \end{aligned}$$

(5)

In terms of saturation, the constraint

$$\begin{aligned} \sum \limits _\beta s^\beta = s^L+s^G = 1 \end{aligned}$$

(6)

holds according to (3) also at any time. Due to the volume-fraction concept, all kinematical as well as physical quantities are defined in the entire control space in an averaged sense.

2.2 Kinematics and deformation measures

In the framework of the TPM, each constituent is assigned an own motion function ${\varvec{\chi }}_\alpha$ with corresponding velocities $\mathbf{v}_\alpha$ and accelerations $(\mathbf{v}_\alpha )^{\prime }\alpha$. Thus,

$$\begin{aligned} {\mathbf{x}}={\varvec{\chi }}_\alpha ({\mathbf{X}}_\alpha ,t) \rightarrow \left\{ \begin{array}{lcl} {\mathbf{v}}_\alpha &{}=&{}\overset{\prime }{\mathbf{x}}_\alpha =\displaystyle \frac{{\mathrm{d}}}{{\mathrm{d}}t}\,{\varvec{\chi }}_\alpha ({\mathbf{X}}_\alpha ,t)\, \\ (\mathbf{v}_\alpha )^\prime _\alpha &{}=&{}\overset{\prime \prime }{\mathbf{x}}_\alpha =\displaystyle \frac{{\mathrm{d}}^{2}}{{\mathrm{d}}t^2}\,{\varvec{\chi }}_\alpha ({\mathbf{X}}_\alpha ,t)\,. \end{array}\right. \end{aligned}$$

(7)

Therein, $\mathbf{x}$ is the current position of the overall material at time t with corresponding reference position $\mathbf{X}_\alpha$ of each $\varphi ^\alpha$ at time $t_0$. Furthermore, $(\,\cdot \,)^\prime _\alpha$ and $(\,\cdot \,)^{\prime \prime }_\alpha$ are the first and second material time derivatives of $(\,\cdot \,)$ following the motion of $\varphi ^{\alpha }$. As is usual in porous-media mechanics, the solid skeleton is described within a Lagrangian setting with respect to the reference configuration via the displacement field

$$\begin{aligned} {{\mathbf {u}}}_S = {{\mathbf {x}}}-{{\mathbf {X}}}_S\,, \end{aligned}$$

(8)

while the pore fluids are described within a modified Eulerian setting governed by the seepage velocities

$$\begin{aligned} {{\mathbf {w}}}_\beta = \overset{\prime }{{\mathbf {x}}}_\beta -\overset{\prime }{{\mathbf {x}}}_S. \end{aligned}$$

(9)

As not only the motion itself is of interest, but also the deformation of the tissue material, the corresponding deformation measures are introduced. Within thermoelasticity, the idea of splitting the material deformation gradient ${{\mathbf {F}}}_{\! S}$ into a thermal part ${{\mathbf {F}}}_{\! S \theta }$ and a mechanical part ${{\mathbf {F}}}_{\! S M}$ is a common practice, compare [27] or [41]:

$$\begin{aligned} {{\mathbf {F}}}_{\! S} = \frac{\partial \varvec{\chi }_S({{\mathbf {X}}}_S,t)}{\partial {{\mathbf {X}}}_S} = {{\mathbf {F}}}_{\! S M} \, {{\mathbf {F}}}_{\! S \theta }. \end{aligned}$$

(10)

Thermally induced volume changes are assumed to be isotropic. Hence, the thermal deformation is given via

$$\begin{aligned} \begin{aligned} {{\mathbf {F}}}_{\! S \theta } = (\mathrm {det}\,{{\mathbf {F}}}_{\! S \theta })^{1/3}\,{{\mathbf {I}}}, \text {where} \mathrm {det}\,{{\mathbf {F}}}_{\! S \theta } = \mathrm {exp}\,\Big [3\,\alpha ^S(\theta -\theta _{0})\Big ]. \end{aligned} \end{aligned}$$

(11)

Therein, $(\theta -\theta _0)$ is the difference of the Kelvin’s temperatures in the current ($\theta$) and the reference ($\theta _0$) configurations, while $\alpha ^S$ is the linear thermal expansion coefficient. Note that a single temperature for all constituents is considered.

As the solid constitutive equations proceed from the basic principles of rational thermodynamics [68], the deformation gradient ${{\mathbf {F}}}_{S}$ is substituted by the right and left Cauchy–Green deformation tensors ${{\mathbf {C}}}_S$ and ${{\mathbf {B}}}_S$, respectively, that will be used along with the spatial Karni-Reiner strain tensor ${{\mathbf {K}}}_S$:

$$\begin{aligned} \begin{aligned} {{\mathbf {C}}}_S={{\mathbf {F}}}_{\! S}^T {{\mathbf {F}}}_{\! S}, \quad {{\mathbf {B}}}_S={{\mathbf {F}}}_{\! S} {{\mathbf {F}}}_{\! S}^T, \quad {{\mathbf {K}}}_S= \displaystyle \frac{1}{2} ({{\mathbf {B}}}_S - {{\mathbf {I}}}). \end{aligned} \end{aligned}$$

(12)

Finally, the spatial velocity gradient ${\mathbf {L}}_\alpha$ and the rate of deformation tensor ${\mathbf {D}}_\alpha$ are defined as:

$$\begin{aligned} {\mathbf {L}}_\alpha = \mathrm {grad}\, \overset{\prime }{{\mathbf {x}}}_\alpha , {\mathbf {D}}_\alpha = \displaystyle \frac{1}{2} \Big (\mathrm {grad}\, \overset{\prime }{{\mathbf {x}}}_\alpha + \mathrm {grad}^T\, \overset{\prime }{{\mathbf {x}}}_\alpha \Big ) \end{aligned}$$

(13)

with $\hbox {grad}\,(\,\cdot \,)$ as the gradient operator applied to $(\,\cdot \,)$. Note that the above rates are solely needed for the evaluation of the entropy inequality of the overall process. Note that the time effect, as described in [64] for Equisetum hyemale, is assumed to originate from the interaction of the solid skeleton with the pore fluids.

2.3 Balance relations

2.3.1 General form

The local balance relations of mass, linear momentum and energy for each constituent $\varphi ^\alpha$ read

$$\begin{aligned} \begin{aligned}&(\rho ^\alpha )^\prime _\alpha + \rho ^\alpha \, \mathrm {div}\, \overset{\prime }{{\mathbf {x}}}_\alpha \,=\, {\hat{\rho }}^\alpha , \\&\rho ^\alpha \, \overset{\prime \prime }{{\mathbf {x}}}_\alpha \,=\, \mathrm {div}\, {{\mathbf {T}}}^\alpha +\rho ^\alpha \,{{\mathbf {b}}}^\alpha + \hat{{\mathbf {p}}}^\alpha , \\&\rho ^\alpha \, (\varepsilon ^\alpha )^\prime _\alpha \,=\, {{\mathbf {T}}}^\alpha \cdot {{\mathbf {L}}}_\alpha -\mathrm {div}\, {{\mathbf {q}}}^\alpha +\rho ^\alpha \, r^\alpha + {\hat{\varepsilon }}^\alpha , \end{aligned} \end{aligned}$$

(14)

where $\hbox {div}\,(\,\cdot \,)$ is the divergence operator corresponding to $\hbox {grad}\,(\,\cdot \,)$. Furthermore, ${\hat{\rho }}^\alpha$ is the so-called mass-production term describing mass interactions. In (14)${}_2$, ${{\mathbf {T}}}^\alpha$ is the partial Cauchy stress and ${{\mathbf {b}}}^\alpha$ the body force per unit mass that is usually substituted by the gravitation vector $\mathbf{g}$, while $\hat{{\mathbf {p}}}^\alpha$ is the direct momentum production that can be interpreted as the local interaction force acting on $\varphi ^\alpha$ through all other constituents of the overall model. In (14)${}_3$, $\varepsilon ^\alpha$ is the mass-specific internal energy that can be to the Helmholtz free energy $\psi ^\alpha$ via the Legendre transformation $\varepsilon ^\alpha = \psi ^\alpha + \theta ^\alpha \eta ^\alpha$, where $\theta ^\alpha$ is the constituent-specific temperature and $\eta ^\alpha$ the mass-specific entropy. The heat influx is denoted by ${{\mathbf {q}}}^\alpha$, $r^\alpha$ is the radiation and ${\hat{\varepsilon }}^\alpha$ is the direct energy production describing the energetic interaction among the constituents.

The balance relations (14) are accompanied by the following constraints for the production terms on the right-hand side of the respective equation:

$$\begin{aligned} \begin{aligned}&\sum _\alpha {\hat{\rho }}^\alpha = 0,\\&\sum _\alpha (\hat{{\mathbf {p}}}^\alpha + {\hat{\rho }}^\alpha \overset{\prime }{{\mathbf {x}}}_\alpha ) = {\mathbf {0}},\\&\sum _\alpha [\,{\hat{\varepsilon }}^\alpha + \hat{{\mathbf {p}}}^\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha + {\hat{\rho }}^\alpha (\varepsilon ^\alpha + \displaystyle \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha )\,] = 0. \end{aligned} \end{aligned}$$

(15)

In (15)$_2$, the total momentum production is split into a direct part $\hat{{\mathbf {p}}}^\alpha$ and a part that is associated with mass production, i. e. ${\hat{\rho }}^\alpha \overset{\prime }{{\mathbf {x}}}_\alpha$. In (15)$_3$, the total energy production is also split into a direct part ${\hat{\varepsilon }}^\alpha$, a part that is associated with the production of linear momentum, i. e. $\hat{{\mathbf {p}}}^\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha$ and a part that is coupling the mass production with the mass-specific internal and kinetic energies, i. e. ${\hat{\rho }}^\alpha (\varepsilon ^\alpha + \displaystyle \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha )$. For a derivation of the balance equations (14) and their restrictions (15), the interested reader is referred to [9] or [23, 24]. Note that the entropy inequality for the overall aggregate will be evaluated in detail when deriving admissible material laws.

For the set of balance equations (14), the restrictions (15) and further considerations within a constitutive setting, the following assumptions and simplifications will be introduced:

• Local thermal equilibrium is considered, which implies a common temperature, such that $\theta ^\alpha =\theta$.

• Materially incompressible solid skeleton, such that $\rho ^{SR}=\rho ^{SR}(\theta )$.

• Materially incompressible liquid water, such that $\rho ^{LR}=\rho ^{LR}(\theta )$.

• Mass interaction solely as cell dehydration between the solid skeleton and the macro-pore water, such that ${\hat{\rho }}^S + {\hat{\rho }}^L = 0$. Although this is a simplification as no sort of conversion of tissue material to water occurs, it is still a valid assumption, since the cells are filled by up to 90 % of water.

• Quasi-static conditions, i. e. $\rho ^\alpha \, \overset{\prime \prime }{{\mathbf {x}}}_\alpha \approx {\mathbf {0}}$, since the processes upon ice formation in Equisetum hyemale are rather slow.

• The gravitation ${{\mathbf {b}}}^\alpha = {{\mathbf {g}}}$ will be neglected, as the dynamics of interest is acting in the horizontal plane, while the gravitation naturally acts vertically along the twig of Equisetum hyemale.

• Note that for the material under study symmetric Cauchy stresses at the macro-scale are considered, i. e. ${{\mathbf {T}}}^\alpha = ({{\mathbf {T}}}^\alpha )^T$, such that the balance of angular momentum is always satisfied.

• The radiation $r^\alpha$ will be neglected.

• The mass-specific kinetic energy $\displaystyle \frac{1}{2}\overset{\prime }{{\mathbf {x}}}_\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha$ is considered small in comparison with the internal energy $\varepsilon ^\alpha$, such that $\displaystyle \frac{1}{2}\overset{\prime }{{\mathbf {x}}}_\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha \ll \varepsilon ^\alpha$.

2.3.2 Materially incompressible thermo-elastic solid

Based on (14)$_1$, an alternative form of the mass balance for a thermo-elastic, materially incompressible solid is obtained by integration. In case of material incompressibility, the effective solid density ${\rho }^{SR}$ does not change under pressure but can change as a result of temperature variations. On the other hand, mass balance equations are given on the basis of partial densities, such as ${\rho }^S=n^S{\rho }^{SR}$, cf. (4). As a result, the integration of (14)$_1$ yields

$$\begin{aligned} \rho ^S = \rho ^S_{0S} \, \left( \mathrm{det}\,{{\mathbf {F}}}_{\! S} \right) ^{-1} \, \mathrm{exp}\, \left( \int _{t_0}^{t} \frac{\hat{\rho }^S}{\rho ^S}\, \mathrm{d}{\tilde{t}} \right) . \end{aligned}$$

(16)

This equation can be further elaborated by considering the split of the solid deformation gradient, compare (10), viz.

$$\begin{aligned} n^S \rho ^{SR} = n^S_{0S} \, \rho ^{SR}_{0S} \, \left( \mathrm{det}\,{{\mathbf {F}}}_{\! SM} \right) ^{-1} \left( \mathrm{det}\,{{\mathbf {F}}}_{\! S\theta } \right) ^{-1} \, \mathrm{exp}\, \left( \int _{t_0}^{t} \frac{\hat{\rho }^S}{\rho ^S}\, \mathrm{d}{{\tilde{t}}} \right) . \end{aligned}$$

(17)

Therein, $n^S_{0S}$ and $\rho ^{SR}_{0S}$ are the initial values of $n^S$ and $\rho ^{SR}$. Furthermore, it is seen that the partial density $\rho ^S = n^S \rho ^{SR}$ can change due to various sources, such as mechanical deformation, thermal deformation or mass production, respectively. In case of materially incompressible solids, the intrinsic density $\rho ^{SR}$ is constant under isothermal conditions, while it can vary under non-isothermal conditions according to [27] yielding

$$\begin{aligned} \begin{aligned} \rho ^{SR} = \rho ^{SR}_{0S}(\mathrm {det}\,{{\mathbf {F}}}_{\! S \theta })^{-1} = \rho ^{SR}_{0S}\, \mathrm {exp}\,[{-3\,\alpha ^S(\theta -\theta _0)}]. \end{aligned} \end{aligned}$$

(18)

By use of this relation, (17) reduces to

$$\begin{aligned} n^S \rho ^{SR} = n^S_{0S} \, \rho ^{SR} \, \left( \mathrm{det}\,{{\mathbf {F}}}_{\! SM} \right) ^{-1} \, \mathrm{exp}\, \left( \int _{t_0}^{t} \frac{{{\hat{n}}}^S\rho ^{SR}}{n^S\rho ^{SR}}\, \mathrm{d}{{\tilde{t}}} \right) . \end{aligned}$$

(19)

Note that in (19), the mass production ${{\hat{\rho }}}^S$ has been substituted by $\hat{n}^S\rho ^{SR}$ meaning that the production of partial density does not mean that the density function $\rho ^{SR}$ is growing or shrinking but the amount of volume that is covered by $\varphi ^S$ with density $\rho ^{SR}$ is changing through $\hat{n}^S$. Thus, while $\rho ^{SR}$ is governed by (18), the solid volume fraction is governed by

$$\begin{aligned} n^S = n^S_{0S} \, \left( {\mathrm{det}}\,{{\mathbf {F}}}_{\! SM} \right) ^{-1} \, {\mathrm{exp}}\, \left( \int _{t_0}^{t} \frac{{{\hat{n}}}^S}{n^S}\, {\mathrm{d}}{{\tilde{t}}} \right) . \end{aligned}$$

(20)

2.4 Constitutive setting

2.4.1 Thermodynamical basis

In order to derive admissible material laws, the second law of thermodynamics in form of the Clausius–Duhem inequality for the overall aggregate is applied. It states the direction, in which a certain process has to occur and reads for the overall aggregate $\varphi$, cf. [24],

$$\begin{aligned} \begin{aligned} \sum _{\alpha } \Big \{&-\rho ^\alpha \big [(\psi ^\alpha )^\prime _\alpha +\theta ^\prime _\alpha \eta ^\alpha \big ] - \hat{{\mathbf {p}}}^\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha - {\hat{\rho }}^\alpha \Big (\psi ^\alpha + \displaystyle \frac{1}{2} \, \overset{\prime }{{\mathbf {x}}}_\alpha \cdot \overset{\prime }{{\mathbf {x}}}_\alpha \Big ) \\&+{{\mathbf {T}}}^\alpha \cdot {{\mathbf {L}}}_\alpha - \frac{1}{\theta } \, {{\mathbf {q}}}^\alpha \cdot \mathrm {grad}\, \theta - {{\mathcal {P}}} (n^\alpha )^\prime _S \Big \} \ge 0, \end{aligned} \end{aligned}$$

(21)

where the saturation condition has already been included as a side condition by means of the time derivative $(\,\cdot \,)^\prime _S$ of (3) multiplied by the Lagrange multiplier ${{\mathcal {P}}}$. In particular, this time derivative yields under consideration of (4) and (14)$_1$ together with the property of material incompressibility of the solid skeleton and the liquid water

$$\begin{aligned} \begin{aligned} \sum _\alpha (n^\alpha )^\prime _S =&\frac{{\hat{\rho }}^S}{\rho ^{SR}} - n^S \, {{\mathbf {I}}}\cdot {{\mathbf {D}}}_S - \frac{n^S}{\rho ^{SR}} \frac{\partial \rho ^{SR}}{\partial \theta } \theta ^\prime _S \\&- \frac{{\hat{\rho }}^S}{\rho ^{LR}} - n^{L} \, {{\mathbf {I}}}\cdot {{\mathbf {D}}}_{L} - \frac{n^{L}}{\rho ^{LR}} \frac{\partial \rho ^{LR}}{\partial \theta } \theta ^\prime _{L} - \mathrm {grad}\, n^{L} \cdot {{\mathbf {w}}}_{L} \\&- n^G \, {{\mathbf {I}}}\cdot {{\mathbf {D}}}_G - \frac{n^G}{\rho ^{GR}}(\rho ^{GR})^\prime _G - \mathrm {grad}\, n^G \cdot {{\mathbf {w}}}_G = 0 \end{aligned} \end{aligned}$$

(22)

with $\mathbf{I}$ as the identity tensor. When (22) is inserted into the entropy inequality (21), the following form is obtained:

$$\begin{aligned} \begin{aligned}&\big ({{\mathbf {T}}}^S + n^S \mathcal {P} \,{{\mathbf {I}}}\big ) \cdot {{\mathbf {D}}}_S - \rho ^S \big (\psi ^S\big )^\prime _S - \rho ^S \left[ \eta ^S - \mathcal {P} \frac{1}{\big (\rho ^{SR}\big )^2} \frac{\partial \rho ^{SR}}{\partial \theta }\right] \, \theta ^\prime _S \\&\quad +\, \big ({{\mathbf {T}}}^L + n^L \mathcal {P} \,{{\mathbf {I}}}\big ) \cdot {{\mathbf {D}}}_L -\rho ^L \big (\psi ^L\big )^\prime _L - \rho ^L \left[ \eta ^L - \mathcal {P} \frac{1}{\big (\rho ^{LR}\big )^2} \frac{\partial \rho ^{LR}}{\partial \theta }\right] \, \theta ^\prime _L \\&\quad +\, \big ({{\mathbf {T}}}^G + n^G \mathcal {P} \,{{\mathbf {I}}}\big ) \cdot {{\mathbf {D}}}_G -\rho ^G \big (\psi ^G\big )^\prime _G - \rho ^G \eta ^G \theta ^\prime _G + \mathcal {P} \frac{n^G}{\rho ^{GR}}\big (\rho ^{GR}\big )^\prime _G \\&\quad -\, \big (\hat{{\mathbf {p}}}^L - \mathcal {P} \, \mathrm {grad}\, n^L\big ) \cdot {{\mathbf {w}}}_L - {\hat{\rho }}^S \overset{\prime }{{\mathbf {x}}}_S \cdot {{\mathbf {w}}}_L - \big (\hat{{\mathbf {p}}}^G - \mathcal {P} \,\mathrm {grad}\, n^G\big ) \cdot {{\mathbf {w}}}_G \\&\quad -\, \frac{1}{\theta } \big ({{\mathbf {q}}}^S + {{\mathbf {q}}}^L + {{\mathbf {q}}}^G\big ) \cdot \mathrm {grad}\, \theta \\&\quad -\, {\hat{\rho }}^S \Big (\psi ^S + \frac{{\mathcal {P}}}{\rho ^{SR}} + \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S \cdot \overset{\prime }{{\mathbf {x}}}_S - \psi ^L - \frac{{\mathcal {P}}}{\rho ^{LR}} - \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_L \cdot \overset{\prime }{{\mathbf {x}}}_L\Big ) \ge 0. \end{aligned} \end{aligned}$$

(23)

Therein, according to [25], extra stresses and the extra momentum productions are identified via

$$\begin{aligned} \begin{aligned}&{{\mathbf {T}}}^\alpha _E = {{\mathbf {T}}}^\alpha + n^\alpha {{\mathcal {P}}} \, {{\mathbf {I}}},\\&\hat{{\mathbf {p}}}^\beta _E = \hat{{\mathbf {p}}}^\beta - {{\mathcal {P}}} \, \mathrm {grad}\, n^\beta . \end{aligned} \end{aligned}$$

(24)

Proceeding from the basic principles of continuum thermodynamics, the set of process variables is chosen to account also for the principle of phase separation [22]. Based on the immiscibility of solid, liquid and gas, the phase-separation principle proceeds from the assumption that each material, as on its microstructure, does only depend on its own constitutive variables. Following this, the free energies $\psi ^\alpha$ and their time derivatives $(\psi ^\alpha )^\prime _\alpha$ read

$$\begin{aligned} \begin{aligned}&\psi ^S = \psi ^S\big (\theta , {{\mathbf {C}}}_S, {\varvec{\mathcal {M}}}^S\big ) \\&\longrightarrow \big (\psi ^S\big )^\prime _S = \frac{\partial \, \psi ^S}{\partial \, \theta } \, \theta ^\prime _S + \frac{\partial \, \psi ^S}{\partial \, {{\mathbf {C}}}_S} \cdot \big ({{\mathbf {C}}}_S\big )^\prime _S, \\&\psi ^L = \psi ^L\big (\theta ,s^L\big ) \\&\longrightarrow \big (\psi ^L\big )^\prime _L = \frac{\partial \, \psi ^L}{\partial \, \theta } \, \theta ^\prime _L + \frac{\partial \, \psi ^L}{\partial \, s^L} \, \big (s^L\big )^\prime _L, \\&\psi ^G = \psi ^G\big (\theta , \rho ^{GR}\big ) \\&\longrightarrow \big (\psi ^G\big )^\prime _G = \frac{\partial \, \psi ^G}{\partial \, \theta } \, \theta ^\prime _G + \frac{\partial \, \psi ^G}{\partial \, \rho ^{GR}} \, \big (\rho ^{GR}\big )^\prime _G. \end{aligned} \end{aligned}$$

(25)

As a result of the above, $\psi ^S$ depends on temperature and deformation, and additionally on the structural tensor ${\varvec{\mathcal {M}}}^S = {{\mathbf {a}}}^S_0 \otimes {{\mathbf {a}}}^S_0$ with one preferred direction, thus indicating a transversely isotropic material. As ${{\mathbf {a}}}^S_0$ is given in the reference configuration, ${{\varvec{\mathcal {M}}}^S}$ is constant yielding $({{\varvec{\mathcal {M}}}^S})^\prime _S\equiv \mathbf{0}$. In addition to $\psi ^S$, $\psi ^L$ depends on temperature and saturation and $\psi ^G$ on temperature and effective density.

As $s^L$ is a function of further kinematic variables, $(s^L)^\prime _L$ has to be found with the aid of (5) and (14)$_1$ together with the property of material incompressibility of the solid skeleton and the liquid water, viz.

$$\begin{aligned} \begin{aligned} \big (s^{L}\big )^\prime _{L}&= \frac{1}{n^F} \left( - \frac{{\hat{\rho }}^S}{\rho ^{LR}} - n^{L} \, {{\mathbf {I}}}\cdot {{\mathbf {D}}}_{L} - \frac{n^{L}}{\rho ^{LR}} \frac{\partial \rho ^{LR}}{\partial \theta } \theta ^\prime _{L} + \frac{s^L {\hat{\rho }}^S}{\rho ^{SR}} - s^L n^S \, {{\mathbf {I}}} \cdot {{\mathbf {D}}}_S \right. \\&\quad \left. - \frac{s^L n^S}{\rho ^{SR}} \frac{\partial \rho ^{SR}}{\partial \theta } \theta ^\prime _S + s^L \mathrm {grad}\, n^S \cdot {{\mathbf {w}}}_{L}\right) . \end{aligned} \end{aligned}$$

(26)

Combining this result with (25)$_2$ and inserting the time derivatives of the free energies into the entropy inequality (23) yield

$$\begin{aligned} \begin{aligned}&\left[ {{\mathbf {T}}}^S_E + n^S \big (s^{L}\big )^2 \rho ^{LR} \, \frac{\partial \psi ^{L}}{\partial s^{L}}{{\mathbf {I}}} - 2 \rho ^S {{\mathbf {F}}}_S \frac{\partial \psi ^S}{\partial {{\mathbf {C}}}_S} {{\mathbf {F}}}_S^T\right] \cdot {{\mathbf {D}}}_S \\&\quad + \left( {{\mathbf {T}}}^{L}_E + n^{L} s^{L} \rho ^{LR} \, \frac{\partial \psi ^{L}}{\partial s^{L}}{{\mathbf {I}}}\right) \cdot {{\mathbf {D}}}_{L} + {{\mathbf {T}}}^G_E \cdot {{\mathbf {D}}}_G \\&\quad - \rho ^S \, \left[ \eta ^S - \frac{{{\mathcal {P}}}}{\big (\rho ^{SR}\big )^2} \frac{\partial \rho ^{SR}}{\partial \theta } - \big (s^{L}\big )^2 \frac{\rho ^{LR}}{\big (\rho ^{SR}\big )^2} \frac{\partial \psi ^{L}}{\partial s^{L}} \frac{\partial \rho ^{SR}}{\partial \theta } + \frac{\partial \psi ^S}{\partial \theta }\right] \, \theta ^\prime _S \\&\quad - \rho ^L \, \left[ \eta ^{L} - {{\mathcal {P}}} \frac{1}{\big (\rho ^{LR}\big )^2} \frac{\partial \rho ^{LR}}{\partial \theta } - s^{L} \frac{1}{\rho ^{LR}} \frac{\partial \psi ^{L}}{\partial s^{L}} \frac{\partial \rho ^{LR}}{\partial \theta } + \frac{\partial \psi ^{L}}{\partial \theta }\right] \, \theta ^\prime _{L} \\&\quad - \rho ^G \, \Big (\eta ^G + \frac{\partial \psi ^G}{\partial \theta }\Big ) \, \theta ^\prime _G + \left( {{\mathcal {P}}} \frac{n^G}{\rho ^{GR}} - \rho ^G \frac{\partial \psi ^G}{\partial \rho ^{GR}}\right) \big (\rho ^{GR}\big )^\prime _G \\&\quad - \frac{1}{\theta } \Big ({{\mathbf {q}}}^S+{{\mathbf {q}}}^L+{{\mathbf {q}}}^G\Big ) \cdot \mathrm {grad}\, \theta \\&\quad -\left[ \hat{{\mathbf {p}}}^{L}_E + \big (s^{L}\big )^2 \rho ^{LR} \, \frac{\partial \psi ^{L}}{\partial s^{L}} \mathrm {grad}\, n^S + {\hat{\rho }}^S \overset{\prime }{{\mathbf {x}}}_S\right] \cdot {{\mathbf {w}}}_{L} - \hat{{\mathbf {p}}}^G_E \cdot {{\mathbf {w}}}_G \\&\quad - {\hat{\rho }}^S \, \left[ \psi ^S + \frac{{{\mathcal {P}}}}{\rho ^{SR}} + \big (s^{L}\big )^2 \frac{\rho ^{LR}}{\rho ^{SR}} \frac{\partial \psi ^{L}}{\partial s^{L}} + \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S\cdot \overset{\prime }{{\mathbf {x}}}_S \right. \\&\quad \left. - \psi ^L - \frac{{{\mathcal {P}}}}{\rho ^{LR}} - s^{L} \frac{\partial \psi ^{L}}{\partial s^{L}} - \frac{1}{2}\overset{\prime }{{\mathbf {x}}}_{L}\cdot \overset{\prime }{{\mathbf {x}}}_{L} \right] \ge 0. \end{aligned} \end{aligned}$$

(27)

According to the Coleman–Noll procedure applied to multiphasic materials, each term of (27) is evaluated separately in order to ensure the independence of the process variables, compare [24]. In particular, the entropy inequality (27) can be split in an equilibrium and a non-equilibrium part, compare [25] for details. The equilibrium part of (27), where the constitutive terms, for which thermodynamical restrictions are used, do not depend on the rates of deformation and temperature, can be found by skipping all terms after line five of (27) and setting the remainder equal to zero. Thus, the terms in front of the rates $\mathbf{D}_S$, $\mathbf{D}_L$ and $\mathbf{D}_G$ as well as in front of $\theta ^\prime _S$, $\theta ^\prime _L$, $\theta ^\prime _G$ and $(\rho ^{GR})^\prime _G$ have to vanish, such that the equation is fulfilled for arbitrary rates in the sense of a sufficient condition. Following this, the Lagrange multiplier is found as

$$\begin{aligned} {{\mathcal {P}}} = \big (\rho ^{GR}\big )^2 \frac{\partial \psi ^G}{\partial \rho ^{GR}} =: p^{GR} \end{aligned}$$

(28)

and can be identified as the excess gas pressure. Note that the fluid stresses ${{\mathbf {T}}}^\beta$ generally consist of equilibrium and non-equilibrium parts; however, only the equilibrium parts are considered here, such that

$$\begin{aligned} \mathbf{T}^G_E=\mathbf{0}\quad \hbox {and}\quad \mathbf{T}^L_E = - n^L s^{L} \rho ^{LR} \frac{\partial \psi ^{L}}{\partial s^{L}}\,{{\mathbf {I}}} \end{aligned}$$

(29)

result by use of (24)$_1$ in the equilibrium fluid stresses

$$\begin{aligned} \begin{array}{c} {{\mathbf {T}}}^G = -n^G p^{GR} \,{{\mathbf {I}}}\,, \\ \displaystyle {{\mathbf {T}}}^{L} = -n^{L} \left( p^{GR} + s^{L} \rho ^{LR} \frac{\partial \psi ^{L}}{\partial s^{L}}\right) \,{{\mathbf {I}}} =: - n^{L} p^{LR}\,{{\mathbf {I}}}\,. \end{array} \end{aligned}$$

(30)

Based on (30)$_2$, a comparison of the gas pressure $p^{GR}$ and the liquid pressure $p^{LR}$ yields

$$\begin{aligned} p^{LR} = p^{GR} + s^{L} \rho ^{LR} \frac{\partial \psi ^{L}}{\partial s^{L}}\,, \end{aligned}$$

(31)

such that the capillary pressure $p^C$ defined as the pressure difference between the non-wetting and the wetting fluid, gas and liquid, is obtained as:

$$\begin{aligned} p^C = p^{GR} - p^{LR} = -\, s^L \rho ^{LR} \frac{\partial \psi ^L}{\partial s^L}. \end{aligned}$$

(32)

The evaluation of the solid skeleton stress $\mathbf{T}^S$ is based on the first line of (27) yielding

$$\begin{aligned} \mathbf{T}^S_E=n^S s^L p^C\,{{\mathbf {I}}} + 2\,\rho ^S\,\mathbf{F}_S\frac{\partial \psi ^S}{\partial \mathbf{C}_S}{} \mathbf{F}^T_S\,, \end{aligned}$$

(33)

such that

$$\begin{aligned} \begin{aligned} \mathbf{T}^S&=-n^S \big (p^{GR}-s^L p^C\big )\,{{\mathbf {I}}} + 2\,\rho ^S\,\mathbf{F}_S\frac{\partial \psi ^S}{\partial \mathbf{C}_S}{} \mathbf{F}^T_S\, \\&=-n^S p^{FR}\,\mathbf{I}\, + \mathbf{T}^S_{E\,\mathrm {mech}}\,, \end{aligned} \end{aligned}$$

(34)

where

$$\begin{aligned} \begin{aligned} p^{FR}:=\big (1-s^L\big )\,p^{GR}+s^Lp^{LR}\,, \\ \mathbf{T}^S_{E\,\mathrm {mech}}=2\,\rho ^S\,\mathbf{F}_S\frac{\partial \psi ^S}{\partial \mathbf{C}_S}{} \mathbf{F}^T_S\,. \end{aligned} \end{aligned}$$

(35)

From (33)–(35), it is seen that the evaluation of the entropy inequality naturally recovers Dalton’s law (35)$_1$, compare [18], with $p^{FR}$ as the effective excess pore pressure, while $\mathbf{T}^S_{E\,\mathrm {mech}}$ from (35)$_2$ yields the basis for the definition of a nonlinear thermo-elastic, anisotropic porous solid material. Note that the case of $p^{FR}=0$ does not necessarily indicate limp tissue cells as only the excess pore pressure vanishes, while the real pressure equals the ambient pressure. Finally, the thermodynamical restrictions for the entropies can be found in case of thermodynamical equilibrium by setting lines three to five of (27) to zero. Thus,

$$\begin{aligned} \begin{aligned}&\eta ^S = \frac{{p^{FR}}}{(\rho ^{SR})^2} \frac{\partial \rho ^{SR}}{\partial \theta } - \frac{\partial \psi ^S}{\partial \theta }, \\&\eta ^{L} = \frac{p^{LR}}{(\rho ^{LR})^2} \frac{\partial \rho ^{LR}}{\partial \theta } - \frac{\partial \psi ^{L}}{\partial \theta }, \\&\eta ^G = - \frac{\partial \psi ^G}{\partial \theta }. \end{aligned} \end{aligned}$$

(36)

In the next step, the non-equilibrium parts of (27) are investigated. Introducing the dissipative expression of the liquid momentum production as

$$\begin{aligned} \begin{aligned} \hat{{\mathbf {p}}}^{L}_{E\,\mathrm {dis}}&\,=\, \hat{{\mathbf {p}}}^{L}_E + (s^{L})^2 \rho ^{LR} \frac{\psi ^{L}}{\partial s^{L}}\, \mathrm {grad}\, n^S + {\hat{\rho }}^S \overset{\prime }{{\mathbf {x}}}_S \\&\,=\, \hat{{\mathbf {p}}}^{L}_E - s^{L}p^C\, \mathrm {grad}\, n^S + {\hat{\rho }}^S \overset{\prime }{{\mathbf {x}}}_S\,, \end{aligned} \end{aligned}$$

(37)

the non-equilibrium terms of (27) can be summarised as the dissipation inequality

$$\begin{aligned} \begin{aligned} \mathcal {D} \,=\,&- \hat{{\mathbf {p}}}^{L}_{E\,\mathrm {dis}} \cdot {{\mathbf {w}}}_{L} - \hat{{\mathbf {p}}}^G_E \cdot {{\mathbf {w}}}_G - \frac{1}{\theta } ({{\mathbf {q}}}^S+{{\mathbf {q}}}^{L}+{{\mathbf {q}}}^G) \cdot \mathrm {grad}\, \theta \\&- {\hat{\rho }}^S \, (\psi ^S + \frac{p^{FR}}{\rho ^{SR}} + \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S\cdot \overset{\prime }{{\mathbf {x}}}_S - \psi ^L - \frac{p^{LR}}{\rho ^{LR}} - \frac{1}{2}\overset{\prime }{{\mathbf {x}}}_{L}\cdot \overset{\prime }{{\mathbf {x}}}_{L}) \, \ge 0. \end{aligned} \end{aligned}$$

(38)

In order to make sure that the inequality (38) is fulfilled for arbitrary processes, the following assumptions are made

$$\begin{aligned} \hat{{\mathbf {p}}}^L_{E\,\mathrm {dis}} \propto - {{\mathbf {w}}}_L, \quad \hat{{\mathbf {p}}}^G_E \propto - {{\mathbf {w}}}_G, \quad {{\mathbf {q}}}^\alpha \,\propto \, -\,\mathrm {grad}\,\theta . \end{aligned}$$

(39)

Furthermore, with regard to the mass interaction, the constraint reads

$$\begin{aligned} \begin{aligned} {\hat{\rho }}^S \propto - \,&\left( \psi ^S + \frac{p^{FR}}{\rho ^{SR}} + \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S\cdot \overset{\prime }{{\mathbf {x}}}_S - \psi ^L - \frac{p^{LR}}{\rho ^{LR}} - \frac{1}{2}\overset{\prime }{{\mathbf {x}}}_{L}\cdot \overset{\prime }{{\mathbf {x}}}_{L}\right) . \end{aligned} \end{aligned}$$

(40)

Based on (39) and (40), the material laws for the respective processes need to be found.

2.4.2 Solid skeleton

The thermo-mechanical extra stress ${{\mathbf {T}}}^S_{E\,\mathrm {mech}}$ has to be found on the basis of (35), such that it describes the required properties of finite deformations, anisotropy, thermoelasticity and mass interactions describing the double porosity. Therefore, the potential energy is additively split into isotropic and anisotropic parts, where the anisotropic part accounts for one preferred direction (the out-of-plane direction for the presented numerical example). Within solid mechanics, the Helmholtz free energy is usually substituted by the strain-energy function $W^S$ per unit reference volume via

$$\begin{aligned} W^S = \rho ^S_{0S} \psi ^S, \quad W^S = W^S_{\mathrm {iso}} + W^S_{\mathrm {aniso}} \end{aligned}$$

(41)

with isotropic and anisotropic parts $W^S_{\mathrm {iso}}$ and $W^S_{\mathrm {aniso}}$ for which an invariant representation has to be found. In the present case, the first invariant and the square root of the third invariant of $\mathbf{C}_S$ are introduced as $I_S = {{\mathbf {C}}}_S \cdot {{\mathbf {I}}}$ and $J_S := \sqrt{III_S}= \sqrt{\mathrm {det}\,{\mathbf {C}}_S}$. Based on various articles, such as [26, 27, 30, 63], the isotropic part of the strain energy is given via

$$\begin{aligned} \begin{aligned} W^S_{\mathrm {iso}} =&\frac{\Lambda ^S}{4} \left( J_S^2-2\,\mathrm {ln}\, J_S -1 \right) -\mu ^S \mathrm {ln}\, J_S + \frac{\mu ^S}{2} \left( I_S-3 \right) \\&- 3 \, \alpha ^S k^S \mathrm {ln}\, J_S \left( \theta -\theta _{0} \right) - \rho ^S_{0S} c_v^S \left[ \theta \, \mathrm {ln}\, \left( \frac{\theta }{\theta _{0}} \right) -\theta +\theta _{0} \right] , \end{aligned} \end{aligned}$$

(42)

where $\Lambda ^S$ and $\mu ^S$ are the Lamé constants, $k^S$ the bulk modulus and $c_v^S$ the specific heat at constant volume. Note that for the description of the thermo-mechanical coupling, the stress-temperature modulus $m^S_\theta = - 3 \,\alpha ^S k^S$ is frequently introduced, cf. [27].

Note in passing that with a decreasing solid volume fraction $n^S$, as in the present case of cell dehydration, the consideration of the so-called compaction point of porous solids [26], namely the local total closure of the pore space, is not necessary. The remaining anisotropic part $W^S_{\mathrm {aniso}}$ of the solid strain energy, cf. [51, 53], has been developed for transversely isotropic materials, such that

$$\begin{aligned} W^S_{\mathrm {aniso}}\,=\, \frac{1}{2} \alpha _{S1} \left( J_{S4} - 1 \right) ^{\alpha _{S2}}. \end{aligned}$$

(43)

Therein, the stretch in the preferred direction is denoted by $J_{S4} = {{\mathbf {a}}}^S \cdot {{\mathbf {a}}}^S$, $\alpha _{S1}$ and $\alpha _{S2}$ are material parameters. Note that ${{\mathbf {a}}}^S=\mathbf{F}_S\,\mathbf{a}^S_0$ is pointing into the preferred direction of the current configuration. The anisotropic extension may account for hard and soft tissues, as it accounts for extension and shrinking as well. In case that a deformation-independent stiffness in the direction of anisotropy has to be considered, $\alpha _{S2} = 2$ needs to be chosen, since the material tangent is constant in that particular case. Based on the considerations so far, the mechanical extra-stress is derived via

$$\begin{aligned} {{\mathbf {T}}}^S_{E\,\mathrm {mech}} = {{\mathbf {T}}}^S_{\mathrm {iso}} + {{\mathbf {T}}}^S_{\mathrm {aniso}}, \end{aligned}$$

(44)

where the respective parts are given by

$$\begin{aligned} \begin{aligned}&{{\mathbf {T}}}^S_{\mathrm {iso}} = \frac{\rho ^S}{\rho ^S_{0S}} \left[ \frac{\Lambda ^S}{2}\,\left( J_S^2-1\right) \, {{\mathbf {I}}} + 2\,\mu ^S \, {{\mathbf {K}}}_S - 3\, \alpha ^S k^S \left( \theta -\theta _{0}\right) \, {{\mathbf {I}}} \,\right] , \\&{{\mathbf {T}}}^S_{\mathrm {aniso}} = \frac{\rho ^S}{\rho _{0}} \, \alpha _{S1} \, \alpha _{S2} \, \left( J_{S4} - 1\right) ^{\alpha _{S2}-1} \, \left( {\mathbf {a}}^S \otimes {\mathbf {a}}^S\right) . \end{aligned} \end{aligned}$$

(45)

Thus, it is easily seen from (16) that the dehydration of the tissue cells is naturally included in (45) by $\hat{\rho }^S$, also compare the similar ansatz by [52] dealing with porous materials with a growing (or shrinking) solid skeleton.

2.4.3 Fluid constituents

Within the macro-pore space, two immiscible pore fluids are present, namely, materially incompressible water (liquid) and materially compressible air (gas). For the liquid component, the thermal dependency needs to be accounted for, where the density anomaly at $3.98^{\circ } \mathrm {C}$ has also to be taken into consideration. A discussion of this effect can be found in [7, 67] with experimentally based approaches on how to handle this temperature dependency. Here, the approach from [67] is used, such that the effective liquid density is given via

$$\begin{aligned} \rho ^{LR}({\tilde{\theta }})\,=\, a_5 \left[ 1-\frac{({\tilde{\theta }}+a_1)^2({\tilde{\theta }}+a_2)}{a_3({\tilde{\theta }}+a_4)} \right] , \end{aligned}$$

(46)

where $a_1-a_5$ are fitting parameters. Note that ${\tilde{\theta }}$ in (46) is given in $\mathrm{^\circ C}$, such that the conversion ${\tilde{\theta }} [\mathrm{^\circ C}] = \theta [\mathrm{K}] -273.15$ has to be used.

For the gaseous component, the equation of state of the ideal gas law in the form of Boyle–Mariotte is applied

$$\begin{aligned} \rho ^{GR}(\theta ,p^{GR})\,=\,\frac{p^{GR}+p_0}{\bar{R}^G \theta }, \end{aligned}$$

(47)

where $\bar{R}^G$ is the specific gas constant and $p_0$ the ambient pressure. Therefore, $p^{GR}$ is the excess pressure compared to the ambient pressure.

Based on the evaluation of the entropy inequality given in (39)$_1$ and (39)$_2$, the following approaches are admissible for the direct momentum productions

$$\begin{aligned} \begin{aligned} \hat{{\mathbf {p}}}^L_{E\,\mathrm {dis}} = -\left( n^L\right) ^2 \rho ^{LR}g\left( {{\mathbf {K}}}^L_r\right) ^{-1} {{\mathbf {w}}}_L, \quad \hat{{\mathbf {p}}}^G_E = -\left( n^G\right) ^2 \rho ^{GR}g\left( {{\mathbf {K}}}^G_r\right) ^{-1} {{\mathbf {w}}}_G \end{aligned} \end{aligned}$$

(48)

with $g=|\mathbf{g}|$. Furthermore, the relative permeabilities ${{\mathbf {K}}}^{L}_r$ and ${{\mathbf {K}}}^{G}_r$ are given as functions of the hydraulic conductivities $\mathbf{K}^L$ and $\mathbf{K}^G$ after [13] as

$$\begin{aligned} {{\mathbf {K}}}^{L}_r = \kappa _r^L {{\mathbf {K}}}^{L}, \quad {{\mathbf {K}}}^{G}_r = \kappa _r^G {{\mathbf {K}}}^G, \end{aligned}$$

(49)

where $\kappa ^L_r$ and $\kappa ^G_r$ are the so-called relative permeability factors given by [13]

$$\begin{aligned} \kappa _r^L=\left( s^L_{\mathrm {eff}}\right) ^{\frac{2+3\lambda _c}{\lambda _c}}, \quad \kappa _r^G=\left( 1-s^L_{\mathrm {eff}}\right) ^2 \left[ 1-\left( s^L_{\text {eff}}\right) ^{\frac{2+\lambda _c}{\lambda _c}} \right] . \end{aligned}$$

(50)

Therein, $\lambda _c$ is the pore-size distribution index and $s^L_{\mathrm {eff}}$ the effective liquid saturation defined after [69] as

$$\begin{aligned} s^L_{\mathrm {eff}} = \frac{s^L-s^L_{\mathrm {res}}}{1-s^L_{\mathrm {res}}-s^G_{\mathrm {res}}} \end{aligned}$$

(51)

with the residual saturations $s^L_{\mathrm {res}}$ and $s^G_{\mathrm {res}}$. As plant tissue is an unsaturated porous material, the effect of capillarity has to be taken into account for the macro-pore space. Choosing the relation between the (effective) liquid saturation $s^L_{\mathrm {eff}}$ and the capillary pressure $p^C$, compare (32), according to [13] to ensure thermodynamic consistency, cf. [27], one makes use of

$$\begin{aligned} s^L_{\mathrm {eff}} = \left( \frac{p_d}{p^C}\right) ^{\lambda _c} \quad \Longleftrightarrow \quad p^C = p_d \left( s^L_{\mathrm {eff}}\right) ^{-\tfrac{1}{\lambda _c}} \end{aligned}$$

(52)

for $p^C \ge p_d$, where $p_d$ is the bubbling pressure. As there are no experimental data available concerning the action of capillarity, $\lambda _c$ and $p_d$ are chosen in order to account for the initial fluid volume fractions and an assumption regarding the pore size distribution, as will be discussed in Sect. 4. In addition, the conductivity tensors $\mathbf{K}^L$ and $\mathbf{K}^G$, measured in $\mathrm{m^3/(m^2 \,s)}$, can be related to the intrinsic permeability $\mathbf{K}^S$, measured in $\mathrm{m^2}$, through

$$\begin{aligned} \mathbf{K}^L=\frac{\rho ^{LR}g}{\mu ^{LR}}{} \mathbf{K}^S, \quad \mathbf{K}^G=\frac{\rho ^{GR}g}{\mu ^{GR}}{} \mathbf{K}^S \end{aligned}$$

(53)

with $\mu ^{\beta R}$ as the shear viscosity of the respective fluid. The intrinsic permeability ${{\mathbf {K}}}^{S}$ in the current configuration can capture an anisotropic (macro-)pore space, where the anisotropy is assumed to be inherent and not caused by the deformation. Following [42],

$$\begin{aligned} {{\mathbf {K}}}^S = \left( \frac{1-n^S}{1-n^S_{0S}} \, \frac{n^S_{0S}}{n^S} \right) ^\pi {{\mathbf {K}}}^{S}_{0S}\,, \end{aligned}$$

(54)

where $\mathbf{K}^S$ is related to the intrinsic permeability $\mathbf{K}^S_{0S}$ of the solid reference configuration with $\pi$ as a nonlinearity parameter covering the change in available pore space, while the term in parentheses accounts for an opening or shrinking of the pore space.

Based on thermodynamic considerations and the derivation of [27], the Helmholtz liquid free energy is given as

$$\begin{aligned} \psi ^L = \frac{\lambda _c \, p^C}{\rho ^{LR}} - c_v^{LR} \left[ \theta \, \mathrm {ln}\,\left( \frac{\theta }{\theta _{0}}\right) - \theta + \theta _{0}\right] + \psi ^L_0, \end{aligned}$$

(55)

where $c_v^{LR}$ is the specific heat at constant volume. Furthermore, the reference energy potential is chosen as $\psi ^L_0 = p^{FR}_0/\rho ^{SR}_0 - \lambda _c p^C_0 / \rho ^{LR}_0 - p^{LR}_0/\rho ^{LR}_0$, such that there is initially no solid mass production . Note that with $(\,\cdot \,)_0$ or $(\,\cdot \,)_{0\alpha }$ indicated quantities refer to the reference state.

The filter velocity $n^L {{\mathbf {w}}}_L$ is determined by solving the quasi-static momentum balance (14)$_2$ by use of the restrictions (24)$_2$, (37), (48)$_1$–(53)$_1$ with respect to $n^L {{\mathbf {w}}}_L$. Thus,

$$\begin{aligned} \begin{aligned} n^L {{\mathbf {w}}}_{L} = - \frac{\kappa _r^L}{\mu ^{LR}} {{\mathbf {K}}} \left( \, \mathrm {grad}\, p^{LR} - \rho ^{LR} {{\mathbf {g}}} - \frac{p^C}{s^L} \mathrm {grad}\, s^L + \frac{{{\hat{\rho }}}}{n^L} \overset{\prime }{{\mathbf {x}}}\right) \end{aligned} \end{aligned}$$

(56)

Therein, $\mathrm {grad}\, p^{LR}$ is the contribution arising from water pressure, $\rho ^{LR} {{\mathbf {g}}}$ from the gravitational potential and $(p^C/s^L)\, \mathrm {grad}\, s^L$ from the interaction between the pore fluids in the macro pores, which is mainly determined by the capillary pressure. Finally, $({{\hat{\rho }}}^S/n^L)\, \overset{\prime }{{\mathbf {x}}}_S$ accounts for the momentum production that is due to mass interaction in terms of cell dehydration, which is orders of magnitude smaller than the contributions from the pressures. Concerning the numerical example of this article that is applied to Equisetum hyemale, a standard cross section of the plant is considered in horizontal direction, such that gravitational forces do not matter, such that this contribution as well as $({{\hat{\rho }}}^S/n^L)\, \overset{\prime }{{\mathbf {x}}}_S$ can be neglected.

According to the introduction of the water potential to address the flow of water in biology-based articles, such as [39, 43] or [45], one might introduce a potential function $h_W$ based on (56), where the gradient of the water potential consists of four parts:

$$\begin{aligned} \mathrm {grad}\,h_W = \mathrm {grad}\, p^{LR} - \rho ^{LR} {{\mathbf {g}}} - \frac{p^C}{s^L} \mathrm {grad}\, s^L + \frac{{\hat{\rho }}^S}{n^L} \overset{\prime }{{\mathbf {x}}}_S. \end{aligned}$$

(57)

For the water potential itself, this implies

$$\begin{aligned} h_W = p^{LR} + U_g + \lambda _c \, p^C + U_{{\hat{\rho }}^S} \end{aligned}$$

(58)

with the water pressure potential $p^{LR}$, the gravitational potential $U_g$, the capillarity potential $\lambda _c \, p^C$ and a potential $U_{{\hat{\rho }}^S}$ accounting for the local momentum production that is due to solid mass production.

For the gas phase, the Helmholtz free energy $\psi ^G$ is based on thermodynamic considerations and is given according to [25, 27] via

$$\begin{aligned} \begin{aligned} \psi ^G =&\bar{R^G} \theta \left[ \mathrm {ln} \left( \frac{\rho ^{GR}}{\rho ^{GR}_{0G}}\right) + \frac{\rho ^{GR}_{0G}}{\rho ^{GR}} -1\right] - c_v^{GR} \left[ \theta \, \mathrm {ln} \left( \frac{\theta }{\theta _{0}}\right) - \theta + \theta _{0}\right] . \end{aligned} \end{aligned}$$

(59)

Therein, $c_v^{GR}$ is the specific heat at constant volume.

As the filter velocity of the pore liquid, the filter velocity $n^G {{\mathbf {w}}}_G$ and accordingly the seepage velocity ${{\mathbf {w}}}_G$ are also determined from the quasi-static momentum balance (14)$_2$ and the restrictions (24)$_2$ and (48)$_2$–(53)$_2$ resulting in

$$\begin{aligned} n^G {{\mathbf {w}}}_G = - \frac{\kappa _r^G}{\mu ^{GR}} {{\mathbf {K}}}^{S} \left( \mathrm {grad} \, p^{GR} - \rho ^{GR} {{\mathbf {g}}} \right) . \end{aligned}$$

(60)

2.4.4 Cell dehydration

In addition to the fluids within the macro-pore space, the dehydration of the tissue cells occurs via micro pores in the cell walls. The impact of this structural property in relation to the water status on the tissue cells has already been tackled by [16, 39, 49, 65, 66], however, rather in the context of cell growth.

The double-porosity feature is also an interesting property from an engineering point of view, inter alia discussed by [21, 33], within the group of Borja, cf. [10, 11, 14, 15, 70], or within the framework of a THM model by [37]. As was outlined in the introduction, these articles proceed from individual fluid constituents at two porosity scales with individual states of motion and mass exchange of a respective fluid across the scales. However, within this contribution, the water is either confined to the solid skeleton or mobile within the macro-pore space. Hence, the introduction of an independent state of motion of water in the micro-pore space is not necessary, but mass exchange is still admissible. Thus, the micro-pore flow is solely considered in terms of a mass-exchange or a mass-production term, respectively, describing the flow of water through the cell wall from the intracellular space to the surrounding intercellular (macro-pore) space. Note that the articles mentioned above regarding double porosity treat the flow already at the continuum scale as the result of a virtual homogenisation. For a mathematically rigorous treatment of the homogenisation, the interested reader is referred to [3] or [55].

Similar to the flow of water at the macroscale represented by (56), the flow of water at the microscale is also considered as a non-equilibrium process. Thus, the evaluation of the entropy inequality (38) suggests that the flow at the microscale through the cell wall, here represented by the mass production ${\hat{\rho }}^S$, is governed by the difference in chemical potentials between the solid skeleton and the water in the macro-pore space. The chemical potential is constituted by three parts, which are the Helmholtz free energy $\psi ^S$, a pressure-driven part ${p^{FR}}/\rho ^{SR}$ and the mass-specific kinetic energy $\displaystyle \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S\cdot \overset{\prime }{{\mathbf {x}}}_S$ of the solid skeleton. In this context, the chemical potential may be also denoted as water potential at the microscale. Thus, the mass production is introduced by the admissible ansatz via

$$\begin{aligned} \begin{aligned} {\hat{\rho }}^S = - \kappa ^{SL} \omega \, L_p \, \left( \rho ^{SR}\right) ^2 \,\left( \psi ^S + \frac{p^{FR}}{\rho ^{SR}} + \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S\cdot \overset{\prime }{{\mathbf {x}}}_S - \psi ^L - \frac{p^{LR}}{\rho ^{LR}} - \frac{1}{2}\overset{\prime }{{\mathbf {x}}}_{L}\cdot \overset{\prime }{{\mathbf {x}}}_{L}\right) . \end{aligned} \end{aligned}$$

(61)

Therein, the hydraulic conductivity $L_p$ of the cell wall accounts for the micro-pore morphology and the physical properties of the pore water. Furthermore, $\omega$ is the specific surface of the cells in natural conditions and is the result of a virtual homogenisation accounting for the cell wall, through which the cell dehydrates in relation to its representative volume. Finally, $\kappa ^{SL} \ge 0$ is the effective permeability coefficient that accounts for the change in effective surface area. This term is introduced as

$$\begin{aligned} \kappa ^{SL} = {\left\{ \begin{array}{ll} \displaystyle 39.0625 \, (n^S-0.1)^2 \, (0.9-n^S)^2 \qquad \text {for}:\, 0.1 \le n^S \le 0.9, \qquad \text {else}.\\ 0&{}\\ \end{array}\right. } \end{aligned}$$

(62)

This ansatz is admissible and accounts for the constraints of dehydration and rehydration, compare Fig. 2.

Fig. 2

Effective permeability coefficient $\kappa ^{SL}$ for a varying state of hydration of the tissue cells. In particular, $\kappa ^{SL}$ has a maximum value of 1.0 and it accounts for the residual solid volume fraction $n^S_{\mathrm {res}}$ and the residual pore space $n^F_{\mathrm {res}}$

The constraints are given by a residual solid tissue material of approximately $n^S_{\mathrm {res}} \approx 0.1$ and a residual pore space of about $n^F_{\mathrm {res}} \approx 0.1$. Thus, the water status of the cells in a freezing environment is described by a local change in solid volume fraction. Note that the kinetic energy terms in (61) will be ignored for the numerical example as they are of negligible size.

2.4.5 Heat flux

Based on (39)$_3$, the dissipative character of the respective partial heat fluxes is accounted for by

$$\begin{aligned} {{\mathbf {q}}}^\alpha \,=\, -n^\alpha H^{\alpha R}\, \mathrm {grad}\, \theta , \end{aligned}$$

(63)

resulting in Fourier’s law for each constituent $\varphi ^\alpha$. Therein, $H^{\alpha R} \ge 0$ is the effective constituent-specific heat conductivity.

3 Numerical solution via the FEM

3.1 Strong formulation of the governing equations

For the thermo-mechanical ternary model under study, numerical computations are based on the mass, momentum and energy balances of the constituents, which basically result in six scalar and three vectorial equations, such that a system of fifteen coupled scalar equations have to be solved. These equations are reduced to only one momentum and one energy balance of the overall aggregate accompanied by mass balances of the individual constituents, such that the general system of fifteen equations reduces to a system of only seven scalar equations.

In a quasi-static setting, the acceleration terms on the left-hand side of (14)$_2$ are dropped. Following this, an addition of the momentum balances of all constituents yields

$$\begin{aligned} {{\mathbf {0}}} \,=\, \mathrm {div}\, ({{\mathbf {T}}}^S_{E\,\mathrm {mech}} - p^{FR}\, {\mathbf {I}}) + \rho \,{{\mathbf {g}}} + {\hat{\rho }}^{S}\,{{\mathbf {w}}}_{L}, \end{aligned}$$

(64)

where (15)$_2$ has been used. This equation represents the momentum balance of the overall medium and acts, at the same time, as the representative of the solid skeleton. Summing up the momentum balances is necessary, as the external forces, as a result of the internal coupling between the constituents, cannot be split into portions acting on the individual components of the model. Note that in the quasi-static setting, the individual momentum balances of the pore water and the air drop off, as they have already been used to determine the seepage or filter velocities of the fluids, compare (56) and (60), that can be inserted into the mass balances of the constituents. Based on (14)$_1$ and (15)$_1$, these balances read

$$\begin{aligned} \begin{aligned}&(\rho ^S)^\prime _S + \rho ^S \mathrm {div} \, ({{\mathbf {u}}}_S)^\prime _S - {\hat{\rho }}^S\,=\, 0, \\&(\rho ^L)^\prime _S + \mathrm {div} \,(\rho ^L{{\mathbf {w}}}_L) + \rho ^L \mathrm {div} \,({{\mathbf {u}}}_S)^\prime _S + {\hat{\rho }}^S \,=\, 0, \\&(\rho ^G)^\prime _S + \mathrm {div} \, (\rho ^G {{\mathbf {w}}}_G) + \rho ^G \mathrm {div} \, ({{\mathbf {u}}}_S)^\prime _S \,=\, 0. \end{aligned} \end{aligned}$$

(65)

As the overall system is governed by only one temperature, the energy balances (14)$_3$ of all constituents can be summed up to yield

$$\begin{aligned} \begin{aligned}&\rho ^S (\varepsilon ^S)^\prime _S + \rho ^L (\varepsilon ^L)^\prime _L + \rho ^G (\varepsilon ^G)^\prime _G - {{\mathbf {T}}}^S_{E \,\mathrm {mech}} \cdot {{\mathbf {L}}}_S + p^{FR} \, \mathrm {div} \,({{\mathbf {u}}}_S)^\prime _S \\&- n^L \mathrm {grad}\, p^{LR} \cdot {{\mathbf {w}}}_L - n^G \mathrm {grad}\, p^{GR} \cdot {{\mathbf {w}}}_G \\&+ \mathrm {div} \, ({{\mathbf {q}}}^S + {{\mathbf {q}}}^L + {{\mathbf {q}}}^G + n^L p^{LR} \, {{\mathbf {w}}}_L + n^G p^{GR} \, {{\mathbf {w}}}_G) \\&+ \hat{{\mathbf {p}}}^L_{E \,\mathrm {dis.}} \cdot {{\mathbf {w}}}_L + n^F p^C \mathrm {grad}\, s^L \cdot {{\mathbf {w}}}_L + \hat{{\mathbf {p}}}^G_E \cdot {{\mathbf {w}}}_G \\&+ {\hat{\rho }}^S (\varepsilon ^S + \displaystyle \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S \cdot \overset{\prime }{{\mathbf {x}}}_S - \varepsilon ^L - \displaystyle \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_L \cdot \overset{\prime }{{\mathbf {x}}}_L) \,=\, 0, \end{aligned} \end{aligned}$$

(66)

where (15)$_3$ has been taken into consideration.

Together with the initial and boundary conditions, the above system of coupled equations defines the strong formulation of the IBVP under discussion. This system of equations has to be solved for the primary variables that can be summarised in the solution vector ${\varvec{y}}\!=\!(\mathbf{u}_S, n^S, p^{LR}, p^{GR}, \theta )^T$ corresponding individually to the following equations:

• Momentum balance of $\varphi ,~{\mathrm{cf.}}$ (64)    $\longrightarrow \quad {\mathbf {u}}_S$

• Mass balance of $\varphi ^S,~{\mathrm{cf.}}$ (65)$_1$    $\longrightarrow \quad n^S$

• Mass balance of $\varphi ^L,~{\mathrm{cf.}}$ (65)$_2$ $\quad \longrightarrow \quad p^{LR}$

• Mass balance of $\varphi ^G,~{\mathrm{cf.}}$ (65)$_3$    $\longrightarrow \quad p^{GR}$

• Energy balance of $\varphi ,~{\mathrm{cf.}}$ (66)    $\longrightarrow \quad \theta$

As the problem is formulated in a quasi-static setting, thus neglecting all inertia terms, the fluid momentum balances have not been used in the numerical setting so far. Thus, the fluid velocities and therewith the filter velocities $n^\beta {\mathbf {w}}_\beta$ reduce to secondary variables that can be found from the fluid momentum balances, compare (56) and (60).

Note that all other secondary variables occurring in the above equations can be determined with the aid of the primary variable and the corresponding constitutive equations.

3.2 Weak formulation of the governing equations

Based on the strong formulation of the problem, it is usually impossible to find a solution for ${\varvec{y}}$. Thus, the strong formulation is transformed to its weak counterpart by use of the Bubnov–Galerkin method. Following this, the governing equations (64)–(66) are multiplied by test functions $\delta {{\mathbf {u}}}_S$, $\delta n^S$, $\delta p^{LR}$, $\delta p^{GR}$ and $\delta \theta$ with a subsequent integration over the solid domain $\Omega$. Specifically, multiplying (64) with $\delta {{\mathbf {u}}}_S$, the weak formulation of the momentum balance of the overall aggregate $\varphi$ reads after some manipulations

$$\begin{aligned} \begin{aligned} \mathcal {G}_{{{\mathbf {u}}}_S}\, \equiv&\int \limits _\Omega ({{\mathbf {T}}}^S_{E\,\mathrm {mech}} - p^{FR} {{\mathbf {I}}}) \cdot \mathrm {grad}\, \delta {{\mathbf {u}}}_S \, \mathrm {d}v \\&- \int \limits _\Omega (\rho \, {\mathbf {g}} + {\hat{\rho }}^S \mathbf{w}_L) \cdot \delta {{\mathbf {u}}}_S \, \mathrm {d}v - \int \limits _{\Gamma ^{{{\tilde{\mathbf {t}}}}}_N} {{\tilde{{\mathbf {t}}}}} \cdot \delta {{\mathbf {u}}}_S \, \mathrm {d}a = 0, \end{aligned} \end{aligned}$$

(67)

where ${{{{\tilde{\mathbf {t}}}}}} = ({{\mathbf {T}}}^S + {{\mathbf {T}}}^L + {{\mathbf {T}}}^G)\, {{\mathbf {n}}}$ is the surface traction acting at the Neumann boundary ${\Gamma ^{{{{{\tilde{\mathbf {t}}}}}}}_N}$, while $\mathbf{n}$ is the outward-oriented unit surface normal. Regarding the primary variables $n^S$, $p^{LR}$ and $p^{GR}$, the corresponding weak formulations of the mass balances read

$$\begin{aligned} \begin{aligned} \displaystyle \mathcal {G}_{n^S} \equiv&\int \limits _\Omega \left[ \big (n^S\big )^\prime _S + \frac{n^S}{\rho ^{SR}} \big (\rho ^{SR}\big )^\prime _S - \frac{{\hat{\rho }}^S}{\rho ^{SR}} - \big ({{\mathbf {u}}}_S\big )^\prime _S \cdot \mathrm {grad}\, n^S \right] \, \delta n^S \, \mathrm {d}v \\&-\int \limits _\Omega n^S \big ({{\mathbf {u}}}_S\big )^\prime _S \cdot \mathrm {grad}\, \delta n^S \, \mathrm {d}v + \int \limits _{\Gamma ^{{\tilde{v}}^S}_{N}} {\tilde{v}}^S \delta n^S \, \mathrm {d}a = 0, \\ \displaystyle \mathcal {G}_{p^{LR}} \equiv&\int \limits _\Omega \left[ \big (n^L\big )^\prime _S \, \rho ^{LR} + n^L \big (\rho ^{LR}\big )^\prime _S + \rho ^L \,\mathrm {div}\, \big ({{\mathbf {u}}}_S\big )^\prime _S + {\hat{\rho }}^S \right] \, \delta p^{LR} \, \mathrm {d}v \\&-\int \limits _\Omega \rho ^L {{\mathbf {w}}}_L \cdot \mathrm {grad}\, \delta p^{LR} \, \mathrm {d}v + \int \limits _{\Gamma ^{{\tilde{m}}^L}_{N}} {\tilde{m}}^L \delta p^{LR} \, \mathrm {d}a = 0, \\ \displaystyle \mathcal {G}_{p^{GR}} \equiv&\int \limits _\Omega \left[ \big (n^G\big )^\prime _S \, \rho ^{GR} + n^G \big (\rho ^{GR}\big )^\prime _S + \rho ^G \,\mathrm {div}\, \big ({{\mathbf {u}}}_S\big )^\prime _S \right] \, \delta p^{GR} \, \mathrm {d}v \\&-\int \limits _\Omega \rho ^G {{\mathbf {w}}}_G \cdot \mathrm {grad}\, \delta p^{GR} \, \mathrm {d}v + \int \limits _{\Gamma ^{{\tilde{m}}^G}_{N}} {\tilde{m}}^G \delta p^{GR} \, \mathrm {d}a = 0. \end{aligned} \end{aligned}$$

(68)

Therein, ${\tilde{m}}^\beta = \rho ^\beta \, {{\mathbf {w}}}_\beta \cdot {{\mathbf {n}}}$ is the mass efflux of the respective fluid through the Neumann boundary ${\Gamma ^{{\tilde{m}}^\beta }_{N}}$. However, for the mass balance of the solid skeleton, there is no volume efflux through the Neumann boundary ${\Gamma ^{{\tilde{v}}^S}_{N}}$ that has to be considered, meaning that ${\tilde{v}}^S = n^S ({{\mathbf {u}}}_S)^\prime _S \cdot {{\mathbf {n}}} = 0$. This is due to the Lagrangian formulation of the solid skeleton. However, the fluid effluxes ${\tilde{m}}^\beta$ are generally non-vanishing based on the fact that the fluids are described by a modified Eulerian description with an Eulerian domain that deforms with the solid skeleton. Finally, the temperature $\theta$ is computed by use of the weak form of the energy balance of the overall aggregate $\varphi$. Thus,

$$\begin{aligned} \begin{aligned} \displaystyle \mathcal {G}_{\theta } \equiv&\int \limits _\Omega \left[ \rho ^S \big (\varepsilon ^S\big )^\prime _S + \rho ^L \big (\varepsilon ^L\big )^\prime _L + \rho ^G \big (\varepsilon ^G\big )^\prime _G - {{\mathbf {T}}}^S_{E \, \mathrm {mech}} \cdot {{\mathbf {L}}}_S + p^{FR} \,\mathrm {div}\, \big ({{\mathbf {u}}}_S\big )^\prime _S \right. \\&\quad \left. - n^L \mathrm {grad}\, p^{LR} \cdot {{\mathbf {w}}}_L - n^G \mathrm {grad} p^{GR} \cdot {{\mathbf {w}}}_G \right. \\&\quad \left. + \hat{{\mathbf {p}}}^L_{E \, \mathrm {dis}} \cdot {{\mathbf {w}}}_L + n^F p^C \mathrm {grad} s^L \cdot {{\mathbf {w}}}_L + \hat{{\mathbf {p}}}^G_E \cdot {{\mathbf {w}}}_G\right. \\&\quad \left. + {\hat{\rho }}^S \left( \varepsilon ^S + \displaystyle \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_S \cdot \overset{\prime }{{\mathbf {x}}}_S - \, \varepsilon ^L - \displaystyle \frac{1}{2} \overset{\prime }{{\mathbf {x}}}_L \cdot \overset{\prime }{{\mathbf {x}}}_L\right) \right] \, \delta \theta \, \mathrm {d}v \\&\quad - \int \limits _\Omega \left( {{\mathbf {q}}}^S + {{\mathbf {q}}}^L + {{\mathbf {q}}}^G + n^L p^{LR} \, {{\mathbf {w}}}_L + n^G p^{GR} \, {{\mathbf {w}}}_G\right) \cdot \mathrm {grad}\, \delta \theta \, \mathrm {d}v \\&\quad + \int \limits _{\Gamma ^{{\tilde{q}}}_{N}} {\tilde{q}} \, \delta \theta \, \mathrm {d}a = 0. \end{aligned} \end{aligned}$$

(69)

Therein, ${\tilde{q}} = ({{\mathbf {q}}}^S + {{\mathbf {q}}}^L + {{\mathbf {q}}}^G + n^L p^{LR} \, {{\mathbf {w}}}_L + n^G p^{GR} \, {{\mathbf {w}}}_G) \cdot {{\mathbf {n}}}$ is the heat influx through the Neumann boundary $\Gamma ^{{\tilde{q}}}_{N}$.

Following the above explanations, the overall problem is based on the abstract formulation

$$\begin{aligned} \boxed { \begin{array}{c} \hbox {``Find a weak solution } {\varvec{y}} \hbox { that solves the function} \\ {{\mathcal {G}}} ({\varvec{y}})=({{\mathcal {G}}}_{\mathbf{u}_S},{{\mathcal {G}}}_{n^S},{{\mathcal {G}}}_{p^{LR}},{{\mathcal {G}}}_{p^{GR}}, {{\mathcal {G}}}_\theta )^T=0 \\ \hbox {and satisfies the boundary conditions''.} \end{array}} \end{aligned}$$

(70)

As the bio-physical processes included in ${{\mathcal {G}}}({\varvec{y}})$ are nonlinear and strongly coupled, a monolithic solution scheme based on the FEM is applied, where the spatial discretisation proceeds from Taylor–Hood elements with a quadratic approximation function for the solid displacement ${\mathbf {u}}_S$ and linear approximation functions for $n^S$, $p^{LR}$, $p^{GR}$ and $\theta$, thus fulfilling the Ladyzhenskaya–Babuška–Brezzi (LBB) condition to avoid stability issues. As we use the finite-element solver PANDAS¹ for our numerical computations, the spatial discretisation is based on an implicit solution scheme proceeding from the Newton–Raphson method, while the temporal domain proceeds from a backward Euler scheme.

4 Application to Equisetum hyemale

The geometry of the underlying plant tissue is depicted in Fig. 3. Therein, Fig. 3a displays a quarter of the cross section in its natural, unfrozen state. The same quarter of the idealised cross section is shown in Fig. 3b, where the boundaries are labelled with Roman numerals, for which the Dirichlet or Neumann boundary conditions are defined according to Table 1. The boundary I is directed to the pith cavity, while the boundary $\mathrm {V}$ is directed to the vallecular canals. Furthermore, three more or less arbitrary points within the cross section of the stem are indicated with their coordinates given in $[\mathrm{mm}]$ as $\mathrm {A}\,(1.3435, 1.3435, 0.025)$, $\mathrm {B}\,(1.452, 1.447, 0.025)$ and $\mathrm {C}\,(1.583, 1.586, 0.025)$.

Fig. 3

Quarter of Equisetum hyemale, shown a in natural condition and b as an idealised but representative cross section with dimensions of the hollow cylinder and vallecular canals with $0.22\,\mathrm {mm}$ diameters. For the discussion of results, points A, B and C are indicated, the boundaries are denoted by Roman numerals ($\mathrm {I}$–$\mathrm {VII}$), cf. Table 1 for the corresponding boundary condition

As was reported in [61], the main locations of ice formation in Equisetum hyemale were found at the surface to the pith cavity and the vallecular canals. The diameter of the pith cavity is $3.8\,\mathrm {mm}$, the vallecular canals with diameter $0.22\,\mathrm {mm}$ are distributed according to Fig. 3b within a cross section with a width of $0.5\,\mathrm {mm}$. The following computations are based on the TPM model described in the preceding sections and make use of a symmetry assumption, such that only one quarter of the cross section has to be considered. The symmetry conditions, meaning no fluxes and no circumferential displacements, along with the other boundary conditions are listed in Fig. 3b and Table 1 for an idealised cross section.

Table 1 Boundary conditions for the boundaries $\mathrm {I}$–$\mathrm {VII}$ of Fig. 3b applied to the initial-boundary-value-problems for Equisetum hyemale

Furthermore, initial conditions have to be defined for the primary variables $\mathbf{u}_S$, $n^S$, $p^{LR}$, $p^{GR}$ and $\theta$ reading

$$\begin{aligned} \begin{aligned} {\mathbf {u}}_{S0} = {\mathbf {0}},~ n^S_{0S} = 0.6,~ p_0^{LR} = -\,0.11 \, [\mathrm {MPa}],~ p_0^{GR} = 0,~ \theta _{0}&= 278.15 \, [\mathrm {K}]. \end{aligned} \end{aligned}$$

(71)

Finally, the Dirichlet boundary condition in terms of temperature ($\theta _D$) at boundary II, cf. Table 1, is prescribed by

$$\begin{aligned} \theta _D = {\left\{ \begin{array}{ll} \displaystyle \theta _{0} - \frac{\theta _{\mathrm {diff}}}{t_{\mathrm {diff}}} \, t \qquad \text {for}:\, t \le t_{\mathrm {diff}}, \\ \displaystyle \theta _{0} - \theta _{\mathrm {diff}} \qquad \quad \text {else} \\ \end{array}\right. } \end{aligned}$$

(72)

with $\theta _{\mathrm {diff}} = 15 \,\mathrm {K}$ and $t_{\mathrm {diff}} = 27,000\,\mathrm{s}$, compare [61]. On the other hand, the onset of ice formation introduces a drop in the water potential, cf. [46]. Hence, this is applied at the boundaries of the extracellular space. Here, the interpolated values from the experimental investigation of [46] have been used for the Dirichlet boundary condition ${\Gamma ^{p^{LR}}_{D}}$ at the vallecular canals (boundary V) as well as at the pith cavity (boundary I):

$$\begin{aligned} p^{LR}_D = {\left\{ \begin{array}{ll} \displaystyle p_0^{LR} + p_{\mathrm {diff}} \, (\theta -\theta _{\mathrm {fr}}) \quad \text {for}:\, \theta < \theta _{\mathrm {fr}}, \\ \displaystyle p_0^{LR} \qquad \text {else} \\ \end{array}\right. } \end{aligned}$$

(73)

with $\theta _{\mathrm {fr}} = 273.15\,\mathrm {K}$ and $p_{\mathrm {diff}} = 0.007\, \mathrm {MPa}$. This shows exemplarily the strong coupling of the underlying THM model, as the pressure drop is a function of the local temperature at the given boundaries of ice formation.

Finally, the model requires the knowledge of all material parameters, as given in (Table 2). The approximated initial effective liquid saturation of $s^L_{\mathrm {eff},\,0} = 0.1$ is implemented by considering a pore size distribution index of $\lambda _c = 0.6$, as for rather poorly sorted porous materials, cf. [35], and an emerging bubbling pressure of $p_d = 0.00237 \, \mathrm {MPa}$. In addition to the mechanical anisotropy, also the hydraulic anisotropy at the macro scale is taken into account by the intrinsic permeability tensor ${\mathbf {K}}^S_{0S}$:

$$\begin{aligned} {\mathbf {K}}^S_{0S} \,=\, \begin{pmatrix} 10^{-8} &{} 0 &{} 0 \\ 0 &{} 10^{-8} &{} 0 \\ 0 &{} 0 &{} 10^{-6} \end{pmatrix} \, {\mathbf {e}}_i \otimes {\mathbf {e}}_j \, [\mathrm {mm}^2]. \end{aligned}$$

(74)

Given (74), the assumption has been made according to the mechanical behaviour of the solid skeleton, that there is one preferred permeability direction, namely perpendicular to the investigated cross section. Moreover, note that due to the high water content of the solid skeleton, the material parameters of the solid skeleton and the pore water are mostly similar.

5 Results and discussion

The following numerical results have been computed with the finite-element solver PANDAS. Of special interest are the temperature-induced effects on the hydraulics in the micro- and macro-pore space in Equisetum hyemale, thus exhibiting the double-porosity feature.

5.1 Deformation

The deformation of an underlying cross section of Equisetum hyemale is shown in Fig. 4 exhibiting a representative in-plane deformation pattern at time t = 50,000 s. This type of behaviour has also been observed by [61] as a result of their experiments indicating an oval shape of the vallecular canals in the deformed configuration, compare also Figs. 4 and 7. This behaviour is the result of the coupled thermo-mechanical deformation and the outflow of water into the vallecular canals and the pith cavity, thus leading to shrinkage of the plant.

Fig. 4

The undeformed cross section of Equisetum hyemale is depicted in light blue; the representative in-plane deformation taken at time t = 50,000 s is shown with the adapted mesh in black. (Color figure online)

5.2 Thermally induced hydraulics in the micro- and macro-pore space

Based on the initial condition of the effective water pressure $p^{LR}$, cf. (71)$_3$, Fig. 5 exhibits the distribution of $p^{LR}$ in the displayed cross section for meaningful time steps. These are 6000 s; 16,000 s; 28,000 s; 36,000 s; 56,000 s; 72,000 s; 116,000 s and 180,000 s increasing from top left to bottom right. When the temperature $\theta$ within the vallecular canals and the pith cavity reaches a value below $273.15\,\mathrm {K}$, the described change in water pressure according to (73) is due to the onset of ice formation. From there on, the pressure drop propagates from boundaries $\mathrm {I}$ and $\mathrm {V}$ through the cross section. The values of the permeability ${\mathbf {K}}^S_{0S}$ determine how fast the pressure drop propagates. The minimal value of $p^{LR}$ is $-\,0.18\,\mathrm {MPa}$, which has been chosen according to the experiments of [46]

Fig. 5

Water pressure $p^{LR}$ distribution within the cross section of Equisetum hyemale for time steps 6000 s; 16,000 s; 28,000 s; 36,000 s; 56,000 s; 72,000 s; 116,000 s and 180,000 s increasing from top left to bottom right regarding time

The change in water pressure has an effect on the macro-pore flow, as depicted in Fig. 6. In particular and according to (56), the pressure distribution in Fig. 5 causes a Darcy flow of water into the vallecular canals and the pith cavity, which is depicted in Fig. 6 in a detailed view. Note that this is a representative pressure and flow pattern and, therefore, displayed in a qualitative manner. It shows the distribution of the filter velocity $n^L {\mathbf {w}}_L$ in the region of interest. Note in passing that the gas exchange at the pith cavity and the vallecular canals is neglected, as stated in Table 1.

Table 2 Material parameters of the ternary TPM plant-tissue model

Fig. 6

Characteristic water pressure $p^{LR}$ distribution and resulting filter velocity $n^L{{\mathbf {w}}}_L$ (blue arrows) in an enlarged view. Smaller arrows indicate a lower filter velocity; longer arrows indicate a higher velocity. (Color figure online)

In Fig. 7, a detailed view of the cross section in frozen condition is given. Therein, the vallecular canal is depicted in light green and the surrounding porous tissue material in dark green. The pith cavity is shown at the lower end, where parts are filled with frozen water, which is also depicted in light green. In frozen condition, it is shown in Fig. 7 that water left the porous tissue into the vallecular canals and the pith cavity, which is qualitatively reflected by the numerical results in Fig. 6. Furthermore, as shown in Fig. 7, the shape of the vallecular canals is oval in frozen condition. Note that this shape change is very pronounced for the displayed cross section; however, this strong expression is not always the case. Nevertheless, it fits qualitatively to the displayed deformation pattern in Fig. 4.

Fig. 7

By courtesy of [59]. Detailed view of the cross section in frozen condition. The indicated vallecular canal is depicted in light green and contains frozen water, whereas the surrounding porous tissue exhibits dark green colour. The pith cavity at the lower end contains partly frozen water, which is also depicted in light green. (Color figure online)

The change in water pressure $p^{LR}$ has also consequences for the micro-pore flow. The flow in the micro-pore space has been considered as cell dehydration in the modelling approach via the mass production term ${\hat{\rho }}^S$ according to (61), which is displayed in Fig. 8. It shows the onset of the micro-pore flow, when the pressure decreases at the surface to the pith cavity and the vallecular canals. As the pressure gradient close to these regions is higher, also shown by the higher filter velocity in Fig. 6, the cell dehydration is more pronounced in this region, cf. Fig. 8. Obviously, the change in dehydration follows the pressure gradient.

Fig. 8

Cell dehydration in terms of ${\hat{\rho }}^S$ for time steps 6000 s; 16,000 s; 28,000 s; 36,000 s; 56,000 s; 72,000 s; 116,000 s and 180,000 s increasing from top left to bottom right regarding time. Cell dehydration is governed by the arising pressure gradient and the involved micro porosity, in particular, the effective hydraulic conductivity of the cell wall

The initial solid volume fraction $n^S_{0S}$ is 0.6. It reduces due to the water loss to a minimal value of 0.1, which corresponds to the (solid) biological tissue material. This shows essentially the dehydration of the solid skeleton as the consequence of the involved thermo-hydro-mechanical processes, which has been recognised as one of the key features of frost hardiness of plant tissues (Fig. 9).

Fig. 9

Evolution of the solid volume fraction $n^S$ for time steps 6000 s; 16,000 s; 28,000 s; 36,000 s; 56,000 s; 72,000 s; 116,000 s and 180,000 s increasing from top left to bottom right regarding time. It approaches a minimum value of 0.1, which corresponds to the amount of biological tissue material in the solid skeleton

The liquid volume fraction $n^L$ is even more insightful, as it arises due to the micro- and macro-pore flow, compare Fig. 10. The initial slight drop in liquid volume fraction $n^L$ at time frame 2 close to the pith cavity as well as the vallecular canals, cf. also Fig. 15, is due to the macro-pore flow, as the drop in water pressure instantaneously leads to an efflux of water into the pith cavity and the vallecular canals, the cell dehydration cannot compensate for this. However, this depends crucially on the relation between the involved hydraulic conductivities at the micro- and the macroscale. The increase in liquid volume fraction for later time steps is due to the increasing cell dehydration.

Fig. 10

Evolution of the liquid volume fraction $n^L$ for time steps 6000 s; 16,000 s; 28,000 s; 36,000 s; 56,000 s; 72,000 s; 116,000 s and 180,000 s increasing from top left to bottom right regarding time. The slight drop in liquid volume fraction at time frame 2, cf. also Fig. 15, is due to the macro-pore flow, the increase in the following time frames is mainly due to the micro-pore flow, i. e. cell dehydration

The described qualitative results are supported by the evaluation of the nodal values at three points A, B and C within the cross section. The nodal values of the temperature $\theta$, the water pressure $p^{LR}$, the solid mass production ${\hat{\rho }}^S$, the solid volume fraction $n^S$ and the liquid volume fraction $n^L$ are given for these points, as introduced in Fig. 3.

The temperature evolution is very similar for all points, indicating that there is just a small temperature gradient over the cross section for a certain time step, compare Fig. 11.

Fig. 11

Temperature $\theta$ close to the pith cavity (Point A), in the middle of the cross section (Point B) and close to the outer surface (Point C). Although the temperature is prescribed at the outer surface, the respective temperature evolution over time is very similar for all locations within the cross section (the curves overlap)

As Point A is at the boundary to the pith cavity, where the water pressure drop is prescribed, the pressure drops linearly, when the temperature is below the freezing point $\theta _{\mathrm {fr}}$ according to (73). When the location of interest is a bit further away from that boundary, as in case of the points B and C, the pressure arises with a time shift as a consequence of the THM processes, but also mainly due to the chosen macro permeability, compare Fig. 12.

Fig. 12

Water pressure $p^{LR}$ for the selected points A, B and C. The pressure at Point A is defined by the prescribed pressure $p^{LR}_D$, when the temperature is below the freezing point $\theta _{\mathrm {fr}}$. The pressure at Point B and Point C is the result of the internally coupled processes, in particular, the chosen permeability at the macroscale

The respective time evolution of the local pressure is one factor that influences the dehydration of the tissue cells, and the other is given by the effective permeability at the microscale. In particular, at Point A, the pressure drop leads to the dehydration of the tissue cells, displayed in Figs. 13 and 14, respectively, and, initially to a drop in liquid volume fraction within the macro-pore space, cf. Fig. 15. But as the dehydration is getting more pronounced, there is an increase in liquid volume fraction. In fact, this illustrates the double-porosity feature of the material, where the decrease is a macro effect, i. e. the cells do not dehydrate fast enough to compensate, and the increase of liquid volume fraction is rather a micro-effect, as the cells dehydrate then faster and overcompensate the efflux of water into the pith cavity in the specific case of the location at Point A. The mentioned time shift of the pressure at locations B and C is also reflected in a similar time shift in the mass production and the corresponding volume fractions. The reduced intensity of the mass production in Fig. 13 indicates a lower gradient of the local pressures.

Fig. 13

The solid mass production ${\hat{\rho }}^S$ at points A, B and C evolves as a function of the effective permeability at the microscale and the pressure difference, which propagates through the cross section, but loses its intensity. Therefore, the maximum value for the solid mass production goes from A over B to C, from the inside, from where the pressure is prescribed, to the outside, with a lower maximum value but longer duration

Fig. 14

Solid volume fraction $n^S$ at points A, B and C. Since the solid mass production ${\hat{\rho }}^S$ has a time shift, that has been shown in Fig. 13, the cells close to the pith cavity dehydrate faster in relation to locations B and C

Fig. 15

Liquid volume fraction $n^L$ at points A, B and C. The double-porosity effect is evident for the time evolution of $n^L$ at Point A, where, after the pressure drop initiates, the water flows immediately into the pith cavity, where the cell dehydration in terms of ${\hat{\rho }}^S$ cannot keep up. For later time steps, the dehydration leads to an increase of water in the macro-pore space. The initial effect is hardly visible at Point B, and not at all at Point C, as the flow in the macro-pore space is small at locations B and C

5.3 Comparison of the bulk flows into the vallecular canals and the pith cavity

Finally, the water efflux into the vallecular canals and the pith cavity are compared. This comparison indicates, where the larger ice body can be found. Thus, the fraction of the total amount of water leaving the porous tissue material into the vallecular canals is 0.73, the efflux into the pith cavity is, consequently, 0.27. This supports the experimental findings by [61] showing that the larger ice body is found in the vallecular canals, as it attracts more water in freezing conditions. The causal volumetric efflux $\mathrm {[mm^3/s]}$ of water into the vallecular canals and the pith cavity is compared over time, which is shown in Fig. 16. It shows that the volumetric efflux into the vallecular canals is higher compared to the volumetric efflux into the pith cavity at any instant of time.

Fig. 16

Comparison of the volumetric flow rate [mm$^3$/s] into the vallecular canals, marked in blue, and the pith cavity, marked in red. (Color figure online)

One reason for the higher total amount of water in the vallecular canals is the higher surface area. Hence, the volumetric efflux of water in relation to its surface in terms of $\mathrm {[mm^3/mm^2 \,s]}$ to the vallecular canals and the pith cavity, respectively, is compared, cf. Fig. 17. This shows that the higher total amount of water in the vallecular canals is not just an effect that can be traced back to its larger surface area, as the volumetric flow per surface area is also larger in this case for the vallecular canals.

Fig. 17

Comparison of the volumetric flow rate per surface area [mm$^3$/mm$^2$ s] into the vallecular canals, marked in blue, and the pith cavity, marked in red. (Color figure online)

6 Conclusion

Extracellular ice formation in plant tissues occurs dispersed in intercellular spaces as well as in larger cavities. Both cases lead to the dehydration of living cells as an essential component of frost hardiness. In order to assess the coupled phenomena that lead to the dehydration of the tissue cells as a consequence of dispersed ice formation, the inclusion of the phase transformation of the water to ice in the macro-pore space is necessary. However, in case that ice formation on internal surfaces and cavities is considered as in this contribution, the impact of the formation of ice on the plant tissue material can be addressed by appropriate boundary conditions.

Therefore, the particular strength of the presented model is the establishment of a modelling framework that allows for addressing the coupled thermo-hydro-mechanical phenomena upon freezing conditions including a proof of thermodynamic consistency, here by application to Equisetum hyemale. The model exhibits a quasi-double-porosity feature in terms of its hydraulics in the micro- and macro-pore space. This type of water management is crucial for the survivability of the plant. This is reflected by the numerical results. First of all, the deformation shows the anticipated pattern. Furthermore, the dehydration of the tissue cells is determined as a function of the ice-induced drop in water pressure at the surface to the pith cavity and the vallecular canals. But this pressure drop causes also the efflux of macro-pore water towards the locations of ice formation. It has furthermore been shown by comparison of these flows that the main location of extracellular ice is found in the vallecular canals.