Hot Cabbage: Thermoregulation in Eastern Skunk Cabbage
BEE 4530: Computer-Aided Engineering: Applications to Biological Processes
Ashley Altera, Anastacia Dressel, Johannes van Osselaer, and Nathan Souchet
May 10th, 2023
©2023 Ashley Altera, Anastacia Dressel, Johannes van Osselaer, and Nathan Souchet
Table of Contents
1.0 Executive Summary.................................................................................................................2
2.0 Introduction..............................................................................................................................2
2.1 Background and Significance.............................................................................................. 2
2.2 Problem Statement...............................................................................................................3
2.3 Model Objectives................................................................................................................. 4
3.0 Methods.....................................................................................................................................5
3.1 Assumptions.........................................................................................................................5
3.1.1 Biological assumptions............................................................................................... 5
3.1.2 Physical assumptions.................................................................................................. 5
3.2 Geometry and Schematic..................................................................................................... 6
3.3 Parameters............................................................................................................................6
3.4 Governing Equations........................................................................................................... 7
3.5 Boundary and Initial Conditions..........................................................................................8
3.6 Heat Generation Term..........................................................................................................9
3.7 Mesh Convergence.............................................................................................................11
4.0 Results..................................................................................................................................... 13
5.0 Validation................................................................................................................................15
5.1 Validation of Generation Term at a Given Temperature.................................................... 15
5.2 Validation of Model Against Changing Data.....................................................................16
6.0 Sensitivity Analysis................................................................................................................ 18
6.1 Varying outer convective heat transfer coefficient............................................................ 18
6.2 Varying Properties of the Spadix and Spathe.....................................................................21
7.0 Discussion............................................................................................................................... 22
7.1 Discussion of Results.........................................................................................................22
7.2 Potential Design Applications............................................................................................23
8.0 Conclusion.............................................................................................................................. 24
9.0 Appendix.................................................................................................................................25
9.1 Derivation of the Heat Generation Term........................................................................... 25
9.1.1 Initial model..............................................................................................................25
9.1.2 Revised Model.......................................................................................................... 26
9.2 Model Parameters.............................................................................................................. 28
9.2.1 Outer convective heat transfer coefficient................................................................ 28
9.2.2 Parameters list...........................................................................................................29
9.2 Computational Costs..........................................................................................................30
10.0 References.............................................................................................................................31
1 of 32
1.0 Executive Summary
Maintaining a safe temperature is essential for nearly all forms of life. Animals have developed
complex systems to maintain an ideal temperature such as shivering and sweating, but such
thermoregulatory systems are not unique animals. The eastern skunk cabbage (Symplocarpus
foetidus) utilizes controlled metabolic heating to maintain a specific elevated temperature in its
flower. Other plant species have evolved the ability to generate heat, but the skunk cabbage is
unique in its ability to maintain its temperature within approximately 1℃ from its ideal 23.6℃,
even in variable and freezing conditions. We have mechanistically modeled skunk cabbage
thermogenesis by developing a finite element model of heat transfer and fluid flow in the
blooming skunk cabbage. The model includes a variable heat generation function which mimics
the natural heating behavior under changing outside conditions, allowing for close analysis of the
thermoregulatory behavior and the energetic costs of heating. In simulated outside temperatures
from 10℃ to -10℃, the generation ranged from 0.07W to 0.17W to maintain ideal temperatures
in the center of the flower. This model displays the impact that changes in temperature can have
on the skunk cabbage due to the dramatically increased energetic cost of heating at colder
temperatures. The complex behavior of the skunk cabbage has evolved to maintain temperatures
within an acceptable range at a minimized energetic cost. Understanding this thermoregulatory
system provides insights that may be used in future development of surfaces used to prevent ice
accumulation and other systems that require low cost thermal regulation.
Keywords: Eastern skunk cabbage, Symplocarpus foetidus, thermogenesis, thermoregulation,
COMSOL, finite element method
2.0 Introduction
2.1 Background and Significance
In early spring, the eastern skunk cabbage (Symplocarpus foetidus) blooms from the mud, often
before the snow has started to melt. The flower has a large hood, the spathe, that encases a thick
stalk that holds and receives pollen, the spadix. As spring progresses, the spathe opens up for
pollination. During blooming, skunk cabbage maintains the temperature of its spadix around
23℃ regardless of outside fluctuating temperatures and even in freezing temperatures (Ito et al.,
2003; Seymour, 2004). This heating is the most precise thermoregulation that has been observed
in plants (Ito-Inaba et al., 2008). Thermoregulation in plants serves multiple essential purposes,
including pollinator attraction, frost protection, and pollen function. Multiple plant species heat
their flowers to increase scent volatility and pollinator attraction. Thermogenesis also provides
direct protection from freezing damage, allowing skunk cabbage to bloom earlier in the season
than most plants. In addition to these heating benefits, maintaining an ideal temperature
dramatically improves skunk cabbage pollen function. Germination rates in skunk cabbage
2 of 32
improve by close to 50% when the spadix is maintained at 23℃ compared to 5℃ warmer or
cooler (Seymour et al., 2009). Although attracting pollinators and melting snow does not require
the exact thermoregulation observed in the eastern skunk cabbage, the specificity of the ideal
temperature for pollen function in addition to the other benefits of heating has driven the skunk
cabbage to evolve extremely accurate thermoregulatory behavior.
Previous research has focused on the thermoregulatory mechanism and biological reasoning for
thermoregulation, but no studies have closely evaluated the heat transfer out from the spadix and
the heat production required to maintain the constant spadix temperature. It has been estimated
that heat generation varies in response to air temperatures, but these estimations did not
accurately calculate heat loss in different conditions. We created a model using the finite element
method to evaluate the heat loss to the outside air and the variable heat generation inside the
spadix. The model takes into account internal radiative heat transfer, outside air convection,
natural convection in the air space between the spathe and spadix, and a dynamic generation term
to accurately describe the heating process. This accuracy allows for close estimation of energy
costs for heating in differing temperatures, and could be used to help analyze the impacts of
widely variable temperatures on the eastern skunk cabbage.
2.2 Problem Statement
The process of heating and maintaining a specific elevated temperature in early spring is unique
to the eastern skunk cabbage and a few closely related species (Figure 1). This thermoregulatory
process is conducted to allow for early blooming and ensure that the spadix is kept at an optimal
temperature for pollen tube growth and germination (Takahashi et al., 2007). Modeling this
process will allow for a greater depth of understanding of eastern skunk cabbage
thermoregulation. We intend to create an accurate 2D axisymmetric model of the eastern skunk
cabbage that actively demonstrates natural spadix behavior. This model can serve as a basis for
analysis of heat generation in the eastern skunk cabbage and elucidate the impact of external
conditions on the energetic cost of thermoregulation.
3 of 32
Figure 1: The top panel shows a three dimensional view of skunk cabbage. Panel a shows an
interior slice of the plant heating in a warmer environment of 10℃. Panel b shows increased heat
generation in a colder environment of -10℃ while maintaining constant spadix temperature.
2.3 Model Objectives
The objective of this project is to model the thermoregulation and heat transfer of the eastern
skunk cabbage exposed to cold outside temperatures to provide insight into the plant behavior
and the biological costs of heat generation. To achieve this overarching goal, we have four
sub-goals.
1. Derive a novel function to represent generation in a changing environment.
2. Create an accurate model of the eastern skunk cabbage thermal regulation.
3. Test the model under variable environmental conditions.
4. Analyze the energetic costs of heat generation from the spadix of the eastern skunk
cabbage using the model.
4 of 32
3.0 Methods
3.1 Assumptions
3.1.1 Biological assumptions
We assume that the plant geometry is axially symmetric and that the spathe is fully closed around
the spadix. This assumption is close to reality in the earlier stages of flower growth when the
spadix is starting to heat. We assume that properties of the spathe and spadix are uniform across
their respective domains. Heat generation in the spadix is simplified to be uniform throughout
the domain for ease of calculation and analysis.
3.1.2 Physical assumptions
We assume all physics acts in 2D due to the axially symmetric geometry. Radiative exchange is
assumed to be uniform across the spadix and spathe surfaces to allow for the use of the
blackbody exchange equations. The convective heat transfer coefficient is constant along the
outer boundary of the spathe which is likely close to the physical reality as airflow is similar
along the surface. Fluid flow inside of the plant is assumed to be laminar since fluid velocities
are very low. The ground below the plant is assumed to be perfectly insulating.
These assumptions allowed us to develop the following schematic to model heat transfer in
eastern skunk cabbage (Figure 2)
5 of 32
3.2 Geometry and Schematic
Figure 2: Model schematic showing various domains. Measurements are for a “young stage” of
the skunk cabbage (Seymour & Blaylock, 1999).
3.3 Parameters
The convective heat transfer coefficient (h) at the spathe boundary was calculated for windless
conditions using the equation for natural convection over a vertical slab. The calculated value
was h = 4.2 [W/m2*K] (see appendix 9.2.1 for details).
Measured thermal properties of the spadix are shown in Table 1 (Takahashi et al., 2007). Spathe
thermal properties, unavailable from literature, were assumed to be identical to those of the
spadix, as tissue from the same species provides a better approximation than data from other
plant tissues (Jayalakshmy & Philip, 2010).
6 of 32
Table 1: Spathe and spadix thermal properties
Thermal conductivity 0.3
[W/(m*K)]
Specific heat [J/(kg*K)] 1660
Density [kg/m3] 1690
Thermal properties of the air domain came directly from the COMSOL materials database and
were checked to be reasonable.
3.4 Governing Equations
We used a transient version of the governing heat transfer equation in cylindrical coordinates for
the spathe domain. This equation excludes both convection and heat generation.
1 ∂ ∂𝑇 ∂ ∂𝑇 ∂𝑇
𝑟 ∂𝑟 (𝑘𝑟 ∂𝑟 ) + ∂𝑧 (𝑘 ∂𝑧 ) =  ρ𝑐 ∂𝑡 (4)𝑝
For the air domain, we used a transient version of the governing heat transfer equation in
cylindrical coordinates that includes convective heat transfer.
1 ∂ (𝑘𝑟 ∂𝑇𝑟 ∂𝑟 ∂𝑟 ) +
∂ ∂𝑇
∂𝑧 (𝑘 ∂𝑧 ) =  ρ𝑐 (
∂𝑇 ∂𝑇 ∂𝑇
𝑝 ∂𝑡
+ 𝑣
𝑟 ∂𝑟
+ 𝑣 ) (5)
𝑧 ∂𝑧
The spadix domain, also in cylindrical coordinates, includes heat generation. This generation
term is added to the governing equation as a variable Q’ in W/m3.
1 ∂
𝑟 ∂𝑟 (𝑘𝑟
∂𝑇 ) + ∂∂𝑟 ∂𝑧 (𝑘
∂𝑇
∂𝑧 ) + 𝑄' =  ρ𝑐
∂𝑇 (6)
𝑝 ∂𝑡
Due to the large difference in spadix and spathe temperatures, we included natural air convection
in the air domain. For this fluid domain the governing equations are the Navier- Stokes equations
in cylindrical coordinates where 𝑢 and 𝑢 represent velocities in the r and z direction.
𝑟 𝑧
Additionally, ρ is considered a function of temperature.
∂𝑢 ∂𝑢 ∂𝑢 ∂𝑝 1 ∂ ∂𝑢 ∂
2𝑢
ρ(𝑇)( 𝑟 𝑟 𝑟∂𝑡 + 𝑢 ∂𝑟 + 𝑢 ∂𝑧 ) =− ∂𝑟 + µ( 𝑟 ∂𝑟 (𝑟
𝑟
∂𝑟 )) +
𝑟 ) (7)
𝑟 𝑧 ∂𝑧2
∂𝑢 ∂𝑢 ∂𝑢 ∂𝑢 ∂2𝑢
ρ(𝑇)( 𝑧∂𝑡 + 𝑢
𝑧
∂𝑟 + 𝑢
𝑧 ) =− ∂𝑝 1 ∂ 𝑧 𝑧
𝑟 𝑧 ∂𝑧 ∂𝑧
+ µ( 𝑟 ∂𝑟 (𝑟 ∂𝑟 ) + 2 ) + ρ(𝑇)𝑔 (8)∂𝑧 𝑧
1 ∂
𝑟 ∂𝑟 (𝑟𝑢 ) +
∂
∂𝑧 (𝑢 ) = 0   (9)𝑟 𝑧
For the fluid dynamics, we assume an incompressible fluid with constant viscosity. For boundary
conditions, we will refer to 𝑢 as the 2D velocity vector composed of 𝑢 and 𝑢 .
𝑟 𝑧
7 of 32
Due to the large difference in temperature between the spadix and the spathe, we included
radiative exchange between these surfaces as shown in Figure 3. Since we do not have access to
COMSOL’s radiative heat transfer between two surfaces, we approximated radiation between the
spadix and the spathe as a blackbody exchange. This radiative exchange was implemented as a
uniform heat flux across the outer surface of the spadix and inner surface of the spathe.
𝑞  =  σ𝐴 𝐹 (𝑇4 − 𝑇4) (10)
1−2 1 1−2 1 2
Since the spadix is nearly fully enclosed by the spathe, 𝐹  ≈ 1. Equation 10 calculates
𝑆𝑝𝑎𝑑𝑖𝑥−𝑆𝑝𝑎𝑡ℎ𝑒
net radiative exchange from surface 1 to surface 2. Implementing these equations as internal
boundary conditions for spadix and spathe surfaces, respectively, yields the following equations:
𝑞  =  − σ𝐴 (𝑇4 − 𝑇4 ) (11)
𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑢𝑟𝑓𝑎𝑐𝑒 𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑝𝑎𝑡ℎ𝑒
𝑞  =  σ𝐴 (𝑇4 − 𝑇4 ) (12)
𝑆𝑝𝑎𝑡ℎ𝑒 𝑆𝑢𝑟𝑓𝑎𝑐𝑒 𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑝𝑎𝑡ℎ𝑒
Figure 3: Radiative heat transfer from the spadix to the spathe as described by equations 11 and
12.
3.5 Boundary and Initial Conditions
The ground was considered an insulating boundary distinguished by a zero flux condition..
8 of 32
∂𝑇
∂𝑧 ｜ =  0 (13)𝑧 = 0
At the center of the plant, our model is assumed to have axisymmetric geometry and therefore
has a zero flux boundary:
∂𝑇
∂𝑟 ｜ =  0 (14)𝑟 = 0, 0 ≤ 𝑧 ≤ 0.0581𝑚
where z is the height of the spadix in meters.
For the surface of the spathe domain that contacts the outside air, we used a convective boundary
with the convective heat transfer coefficient (ℎ) and outside temperature ( 𝑇 ) as a function of
∞
time:
− 𝑘 · ∂𝑇∂𝑟 ｜ =  ℎ (𝑇 − 𝑇 (𝑡)) (15)𝑟=𝑒𝑑𝑔𝑒 𝑜𝑓 𝑠𝑝𝑎𝑡ℎ𝑒 𝑆𝑝𝑎𝑡ℎ𝑒 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 ∞
At the axis of symmetry, the velocity flux was zero due to symmetry.
∂𝑢
∂𝑟｜ =  0 (16)𝑟 = 0, 0 ≤ 𝑧 ≤ 0.0581𝑚
where z is the height of the spadix in meters.
All other fluid dynamics boundaries of the air domain were set as no slip walls.
𝑢｜ =  0 (17)
𝑏𝑜𝑢𝑛𝑑𝑎𝑟𝑦
Temperatures in all domains were initialized to be equal to the initial outside temperature of
10℃.
𝑇 =  10℃ (18)
𝑖
The fluid dynamics were initialized with zero velocity at 1 atm of pressure.
𝑢 = 0 (19)
𝑖
𝑝 = 1 [𝑎𝑡𝑚] (20)
𝑖
3.6 Heat Generation Term
In order for the model to accurately represent the eastern skunk cabbage’s behavior in a changing
environment, a dynamic generation term was needed. Skunk cabbage is able to increase its heat
generation in cold environments and decrease the term as temperatures warm with the purpose of
maintaining an optimal inner spadix temperature. Engineering models for generation have
previously been published hypothesizing that generation scales with the rate of temperature
change of the spadix. In other words, generation (Q) defined in units of watts, is presented as a
∂𝑇
function of 𝑆𝑝𝑎𝑑𝑖𝑥∂𝑡 (Takahashi et al., 2007). For our model, the implementation of this approach
9 of 32
presented a high level of complexity. In order to avoid this difficulty, a novel function describing
generation in Symplocarpus foetidus was designed.
To design an appropriate Q term, data collected by Ito and colleagues in 2004 was analyzed. Two
key characteristics of how Symplocarpus foetidus responds to rapid jumps in temperature formed
the basis for our mathematical approximation of generation. The first is overshooting behavior,
where the plant “overheats itself” when recovering after a temperature drop. The second is a
“range of temperature stability”, where the plant is satisfied with the spadix being at a range of
temperatures around 23.63℃ as opposed to one exact temperature.
Using these observations heat generation was modeled using a piecewise differential equation
(see appendix for details). This approach displayed both overshooting behavior and a range of
temperature stability when used in the model.
To visualize our piecewise differential equation for Q, ∂𝑄∂𝑡 was graphed in Figure 4 against
potential spadix temperatures.
Figure 4: Potential spadix temperatures and the resulting change in generation at that spadix
temperature, as formulated by our piecewise function.
This visualization showcases the spadix sensitivity range where our function displays a nearly
zero slope to span the range of temperatures between 22.73℃ and 24.53℃. With this generation
10 of 32
model we hope to accurately represent how Symplocarpus foetidus is able to adapt to fluctuating
outside temperatures.
3.7 Mesh Convergence
The model calculations are run using the COMSOL physics generated mesh with 2126 elements
shown below (Figure 5). This setting has a maximum difference in element size (maximum
element size - minimum element size) of 5.29mm for the spathe and spadix domains and
1.35mm in the air domain. Although this mesh may not be fully converged in all cases, the
differences in the results between the finest and coarsest mesh are minute (±0.3K), as shown in
Figure 6, and computation time dramatically increases for a finer mesh. Due to these
considerations, we believe that this mesh is the best choice for this model and provides a result
that is accurate within the needs of the model.
Figure 5: Final 2126 element mesh used for all validation and final model runs. Element size was
reduced in areas where velocity and temperature gradients are largest and where geometry is
most complex to ensure accuracy.
11 of 32
Figure 6: Mesh convergence analysis performed by comparing spadix and spathe temperatures
for different mesh sizes. Spadix and spathe temperature data was collected after 1000 seconds.
Mesh appears to converge at a maximum element size difference (maximum element size -
minimum element size) of ~4mm.
12 of 32
4.0 Results
The model was run under windless conditions for 100 minutes with changing air temperatures
starting at 10℃ before dropping to -10℃ and increasing back to 5℃. In order to initialize
temperature, the model was run at a constant 10℃ temperature for 50 minutes to ensure all
domains were at a stable temperature at the start of the simulation. These conditions were
selected to showcase the model’s function in adjusting generation to changing temperature and
provide insight into the energy requirements at varying temperatures.
Figure 7: A temperature profile of the skunk cabbage after 100 minutes. The spadix maintains a
temperature between 26℃ and 24℃. The spathe has a temperature of around 11℃, much closer
to the ambient temperature (10℃).
13 of 32
Figure 8: Air speed in between the spathe and the spadix after 100 minutes. The largest amount
of movement is directly above the spadix, which is where heat is generated, as hot air flows up
from the spadix before cooling and flowing down along the inside of the spathe.
As displayed in Figure 7 and Figure 8, the heat of the spadix causes natural convection in the air
domain, leading to higher air velocities and increased temperatures above the spadix as air near
the spadix is heated and rises. Figure 9 displays the average spadix temperature and spadix heat
generation over 100 minutes with changing outside air temperatures. As outside air temperature
changes, the spadix temperature is initially affected and begins to drop. When the spadix
temperature gets too low, generation begins to ramp up and stabilizes the spadix temperature
back into the optimal range. Spadix temperature stabilizes approximately 20 minutes after
outside temperatures change, as generation does not start to change significantly until spadix
temperature has exited the optimal range. This behavior is clearly evident in the outside
temperature change from -10℃ to 5℃ where generation does not start to decrease for nearly ten
minutes.
14 of 32
Figure 9: Values predicted by the model for spadix temperature and spadix heat generation over
time as a result of changes in the air temperature. The skunk cabbage aims to be in the optimal
range of spadix temperature, and changes the spadix heat generation accordingly.
Beyond the dynamics of generation change with time, the magnitude of generation provides
insight into the heat loss in the different outside temperature conditions. Notably, even with a
240% increase in generation, spadix temperature stabilizes at a lower temperature when the
outside temperature is -10℃ than in the 10℃ condition. This comparison shows the dramatic
increase in heat loss at lower outside temperatures, as the added energy from the increased
generation is lost to the outside air.
5.0 Validation
5.1 Validation of Generation Term at a Given Temperature
In order to validate the magnitude of the spadix generation term, generation was compared to
approximations of generation based on experimental measurements of CO2 production.
15 of 32
Constant generation inside the spadix was derived from CO2 production data for an outside
temperature of -10℃ from Seymour, 2004 together with the relation between CO2 production
and heat generation stated in Takahashi et al., 2007.
Experimental generation estimation:
For a plant in a windless environment at -10℃
Energy production for average 4.36g spadix at -10℃ = 0.6W (Seymour, 2004)
Heat production of spadix = 0.6W * 0.2 (20%) = 0.12W (Takahashi et al., 2007)
In the above calculation, it should be noted that the assumption that 20% of total energy
production is used for heat production is not fully validated and was taken from Takahashi et al.
2007. Comparing the calculated values to our model, in identical conditions, our model stabilizes
at a generation of 0.16W. This generation is slightly higher than the estimated real world values;
however, this calculation allows us to confirm our model's heat generation is reasonable and
within the logical range. Further validation would provide greater insight, but experimental data
of spadix generation is exceptionally rare and difficult to measure.
5.2 Validation of Model Against Changing Data
Figure 10: A comparison of experimental and model predicted spadix temperatures with
changing air temperature graphed against time for approximately 770 minutes (Takahashi et al.,
2007).
16 of 32
In Figure 10, our model’s predicted spadix temperature is plotted against collected data
(Takahashi et al., 2007). The predicted spadix temperature follows the real data values very
closely up until about 650 minutes when we notice a deviation. At this time point, the model
predicted the spadix temperature rises about 1 degree above recorded temperature. This
discrepancy can be explained from the inaccuracies in our heat generation term. Generation in
the model is programmed to settle at any value within a range for the spadix optimal
temperature. In reality, the plant does not settle at an arbitrary value within this range, but settles
based on the energy required to stay at that temperature. This difference caused the model to
plateau at a higher temperature and as a result get hotter than the experimental data between 650
and 700 minutes. Additionally, the model’s overshooting/undershooting behavior seems to be
slightly more extreme than the experimental data. This indicates that the constants in the
generation function may need to be fine-tuned to match the plant’s behavior more exactly.
Figure 11: A comparison of the experimental and model predicted spadix temperatures with
changing air temperature graphed against time for approximately 770 minutes (Seymour, 2004).
No parameters were changed from validation with the previous experimental data (Takahashi et
al., 2007)). In Figure 11, there are major discrepancies between the experimental data and the
model, in particular the difference at 380 minutes. These differences are likely due to inaccuracy
17 of 32
in the generation function when spadix temperatures move far outside of optimal. Due to the
large decrease in outside temperature at 380 minutes from 2℃ to -14°C spadix temperature
decreases far below the optimal range, and modeled generation increases more quickly than is
biologically possible in the skunk cabbage.
In each of the experiments conducted in literature, the outside temperature was changed
dramatically (25 to -14℃), and scientists recorded the change in the plant’s spadix temperature
shown in gray (Figures 10 and 11). The orange represents our COMSOL model’s results under
the same temperature changes. As can be seen, the resulting change in temperature of the spadix
from the model is very similar to the recorded reality. Although not perfect (discrepancies
discussed above), we believe this indicates that our implementation of the generation term in our
heat transfer model accurately mimics the plant’s response to the changes in spadix temperature
based on the outside temperature. However, the plant’s adjustment of the spadix temperature to
within the optimal temperature range occurs more quickly than the experimental data, especially
after dramatic temperature drops. This inaccuracy is likely due to the rate at which the model
generation changes when spadix temperatures are far removed from the ideal. Based on this data,
it appears that change in modeled generation is faster at these extremes than the real plant’s
response indicating that the constants of the generation function may need to be adjusted.
6.0 Sensitivity Analysis
6.1 Varying outer convective heat transfer coefficient
We wanted to see how changes in wind speed from the outside air affected the amount of
generation required for the plant to maintain its spadix temperature. This change in wind was
modeled as a change in the convective heat transfer coefficient (h). Based on the model, it is
observed that, for h values above 15 [W/(m²*K)], a threshold is reached where the amount of
generation does not change drastically as the h term continues to increase (Figure 12). This
threshold changes for different ambient temperatures, and varies more when the temperature
changes more rapidly.
18 of 32
Figure 12: Change in generation of the spadix over time with different h values as the ambient
temperature changes. As the h value increases, Q reaches a stable value at different ambient
temperatures.
By integrating the generation over time, we can determine the total amount of energy in joules
needed for the skunk cabbage to maintain its optimal temperature range for the given timestep
and compare these values (Figure 13).
19 of 32
Figure 13: Total energy required for the plant to maintain optimal spadix temperature range at
different convective heat transfer (h) values. The h values were varied from the lowest predicted
value for windless conditions (h = 3) and increased by 6 W/(m2*K) for every subsequent test to
model the effect of increasing wind speeds. As h increases, the difference in energy between h
values becomes smaller.
As displayed in Figure 13, total energy required to maintain temperature increased as h
increased. At a high h, the surface temperature of the spathe is approximately equal to the outside
air temperature. As h increases further it will only approach outside air temperature meaning that
increases in h above 30 will have a minimal effect on total heat generation required.
20 of 32
6.2 Varying Properties of the Spadix and Spathe
Figure 14: Percent change in total energy required to heat Symplocarpus foetidus from a 10%
increase or decrease in each spadix parameter.
Figure 15: Percent change in total energy required to heat Symplocarpus foetidus from a 10%
increase or decrease in each spathe parameter.
21 of 32
Our sensitivity analysis on the spathe and spadix provided several insights about the relationship
that these properties have to the plants overall ability to maintain its temperature. For the spadix,
it is clear that thermal conductivity has the largest effect on required energy. This effect is likely
due to the impact of thermal conductivity on the time taken by generated heat to diffuse from the
center of the spadix to the surface where heat is quickly dissipated through convection.
Our spathe analysis showed that all three parameters had very little effect on total required
energy compared to spadix parameters. Of the three, it seems that the specific heat and density
has the largest effect on total energy required, indicating that the main function of the spathe is
not as an insulator but as a heat sink. This data also helped to validate that our choice to assume
spathe parameters were equal to spadix parameters was not a bad assumption since spathe
parameter values have such a low impact on total required energy.
7.0 Discussion
7.1 Discussion of Results
The results and validation demonstrate the effectiveness of this model in portraying the reality of
how the eastern skunk cabbage responds to changes in the ambient temperature. The generation
function computed the change in heat generation as a function of average spadix temperature at
any given time. This effectively stabilized the average spadix temperature within the specified
optimal temperature range and closely follows natural behavior. Although average spadix
temperature is accurate to experimental measurements, the center of the spadix tends to be
warmer (approximately 26℃) than the optimal temperature range. This may not be accurate as
experimental findings did not specify whether temperature readings were taken on the spadix
surface or inside the spadix.
Outside of the spadix, the fluid simulation accurately modeled natural convection. In Figure 7,
the air just around and directly above the spadix is 4-8℃ warmer than the air close to the spathe.
This temperature difference causes the fluid flow displayed in Figure 8. The air directly above
the spadix has a much higher velocity than other parts of the air domain, as here the air is heated
the most and becomes less dense. The rest of Figure 8 shows the circular airflow pattern that
developed from these differences in air density. This is what we were expecting to see, as hot air
rises and then falls as it is cooled. The magnitude of air velocity within the plant has never been
recorded, but the velocities reached in the model are within reasonable expectation for the
relatively small temperature differences between the spadix and the spathe. This airflow
increases the heat loss at the surface of the spadix as the convective flow quickly pushes hot air
away to the colder spathe surface.
22 of 32
Figure 9 displays the behavior of our model during two drastic temperature changes (a 20℃
decrease and a 15℃ increase) over 100 minutes. This exhibited the mechanics of the generation
term and how it interacted with other elements of our model. When the air temperature changes,
the spadix temperature exits the range of ideal temperatures, found in the literature and dictated
in our generation equation, before the heat generation starts to change. This is precisely what is
observed in reality, and the spadix temperature trends closely match experimental data (Figures
10 and 11). Although there is little experimental data outlining the amount of power generated by
the spadix with changing temperatures, the magnitude of the values closely follows the
approximate generation measurements that were used for comparison. Takahashi et al. 2007,
hypothesized that skunk cabbage has a sensory delay where temperature changes are not sensed
for approximately 30 minutes (Takahashi et al. 2007). We were unable to implement this aspect
of dynamic heat generation; despite this, our dynamic generation term is able to accurately
model the experimental data, which implies that the sensory delay may not be the real biological
mechanism. There were no models created with the sensory delay, but it would be interesting to
compare a model with this delay to our current model to elucidate the true biological mechanism.
The 2D axisymmetric geometry was a simplification necessary for us to be able to implement the
physics that needed to be included for a functioning model. However, this is likely not as
accurate as a 3D geometry. The young skunk cabbage has a closed spathe, but a majority of the
time that the skunk cabbage is thermoregulating its spadix, the spathe is open, and there’s a gap
that exposes the inside of the plant to the outside air (Seymour & Blaylock, 1999). This would
require a 3D geometry to account for the opening, and would change the movement of the air
inside the plant. The skunk cabbage would lose more heat due to the opening of the plant, and it
would be interesting and valuable to see the energy costs associated with this adjustment.
However, with the time and resources we had available to us, our working model and results are
reasonable and present a lot of insight into the thermogenesis of the eastern skunk cabbage.
7.2 Potential Design Applications
Ice formation, accumulation, and adhesion to surfaces can cause both economic and safety
issues, such as electrical failures and dangerous roads. Keeping critical surfaces ice free can be
costly and challenging, and freezing damages have cost the United States 32.8 billion USD in the
last ~40 years (Smith, 2022). Heating an icy surface forms a thin layer of water between the ice
and the surface that can allow the ice to slide away. Combining heating with hydrophobicity,
forces this water layer to disconnect from the surface and can effectively prevent ice buildup.
Thermoresponsive coatings are among the most common coatings to prevent ice buildup; these
materials will respond to stimuli to generate heat. To be successful, coatings need to be able to
maintain an optimal temperature and harness hydrophobic properties (Shamshiri et al., 2021).
Thermogenesis in plants, namely the eastern skunk cabbage, has inspired a lot of the research
23 of 32
into thermoresponsive coatings. Our model is a clear example of a thermoresponsive system that
can maintain temperatures while minimizing energetic costs. Although these materials have a lot
of potential, the technology itself is still in the early stages of development. Biomimetics offers
an avenue to emulate the techniques nature uses to repel ice in order to improve human-made
systems (Shamshiri et al., 2021). Future research may allow for the development of
thermoregulated coatings that dynamically generate the heat required to prevent ice
accumulation mimicking the processes of Symplocarpus foetidus.
8.0 Conclusion
This model of the eastern skunk cabbage provided major novel insight into the complex behavior
of thermoregulation. The dynamic heat generation term allowed the model to closely reflect the
plant’s response to dramatic temperature changes while providing data of the generation
required. This data allows for accurate analysis of the energetic impacts of thermoregulation
which is essential in understanding the evolution and fitness benefits of this behavior. The
biological dynamics of variable heat generation and temperature sensing in the skunk cabbage
remain unknown, but this model provides insight into these functions. Our results imply that
generation change likely functions similarly to the model implementation.
For future directions, further analysis and improved models that account for the 3D geometry
and the open spathe of the mature flower can provide additional insights and improved accuracy.
The assumptions made to simplify this model allowed it to be completed within one semester,
but future projects with more available time and computing power can expand on the foundation
we have laid. Pairing modeling with experimental data can further improve the simulation by
allowing accurate measurements of real plants. These improved models and the results discussed
in this report may help researchers reach a deeper understanding of skunk cabbage
thermoregulation and inspire innovation in human thermoregulatory technologies.
24 of 32
9.0 Appendix
9.1 Derivation of the Heat Generation Term
9.1.1 Initial model
As the spadix temperature drops below the spadix optimal temperature (𝑇 ), the skunk
𝑜𝑝𝑡𝑖𝑚𝑎𝑙
cabbage cannot immediately change its generation. As the plant cools due to lower outside
temperature, generation begins to “ramp up”, heating the plant to get back to optimal
temperature. This increase in generation continues until optimal temperature is reached, at which
time, the plant begins to dial back. In biology, instantaneous change is impossible, so once the
plant realizes optimal temperature has been reached there is a delay in the time required to ramp
down generation causing spadix temperature to often rise above 𝑇 . This phenomenon can
𝑜𝑝𝑡𝑖𝑚𝑎𝑙
be described as overshooting behavior and can be modeled using a differential equation for the
rate of change of generation (Q).
When skunk cabbage is far from its optimal temperature, generation changes more rapidly than
when temperatures are close to optimal (Takahashi et al., 2007). Symplocarpus foetidus does not
have a sensor for outside temperature and therefore its response to a changing environment has
to be based on changes in spadix temperature. This implies that the rate of change in generation
is proportional to the difference between current spadix temperature and optimal temperature.
∂𝑄
∂𝑡  ∝ |𝑇 − 𝑇 | (A1)𝑜𝑝𝑡𝑖𝑚𝑎𝑙 𝑠𝑝𝑎𝑑𝑖𝑥
These observations lead to the following differential equation, where 𝑘 is the generation scaling
1
constant and 𝑇 is the average temperature of the spadix at any given time t.
𝑆𝑝𝑎𝑑𝑖𝑥
∂𝑄
∂𝑡  =  𝑘 * (𝑇 − 𝑇 (𝑡)) (A2)1 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑝𝑎𝑑𝑖𝑥
In order to calculate values for our differential equation, we began by determining the spadix
optimal temperature. Using plot digitizer, stabilized spadix temperature at four different outside
air temperatures (10.84, 12.86, 15.80, and 15.84℃) were recorded from data collected by Ito et
al. 2004. These values were averaged to generate one spadix average optimal temperature
determined to be 23.63℃.
𝑇 =  23. 63℃ (A3)
𝑜𝑝𝑡𝑖𝑚𝑎𝑙
25 of 32
Next, we needed to find a constant (𝑘 ) that would reflect real behavior observed in experimental
1
data. Experimental results from Ito et al. 2004 helped to guide our understanding of how quickly
the skunk cabbage would change its generation yielding a generation scaling constant, 𝑘 .
1
𝑘 =  0. 00005 (A4)
1
Plugging these values into our differential equation and converting to Kelvin yields our first
model of generation:
∂𝑄
∂𝑡  =  0. 00005 * (296. 78− 𝑇 (𝑡)) (A5)𝑆𝑝𝑎𝑑𝑖𝑥
Additionally, in the biological context generation (Q) can never be a negative value and thus a
hard limit was set to specify that regardless of ∂𝑄∂𝑡 , 𝑄 ≥  0.
This preliminary equation for generation displayed clear overshooting behavior but failed to
display accurate spadix temperatures in more complex situations. This failure was due to a lack
of accommodation for the spadix’s allowed range of temperature stability.
9.1.2 Revised Model
Depending on outside air temperature, spadix temperature will stabilize at different values within
an acceptable range. This pattern follows a trend where the lower the outside air temperature, the
lower the stabilization temperature value of the spadix. We hypothesize this behavior is due to
the fact that Symplocarpus foetidus performs a minimization calculation weighing the cost of
energy generation (Q) against the benefits of having the spadix at a higher, more optimal
temperature. Unfortunately, a calculation to replicate the cost-benefit analysis done by the plant
has never been attempted in literature and the required values are unavailable.
As an alternative, we designed our revised function for Q to stabilize in a range of acceptable
spadix temperatures defined by a spadix maximum optimal temperature (𝑇 ) and a
𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑎𝑥−𝑂𝑝𝑡
spadix minimum optimal temperature (𝑇 ). This range was calculated based on our
𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑖𝑛−𝑂𝑝𝑡
already established value for optimal spadix temperature (𝑇 ) and a spadix temperature
𝑜𝑝𝑡𝑖𝑚𝑎𝑙
sensitivity range (𝑇 ). Spadix sensitivity range was presented in the literature
 𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑒𝑛𝑠𝑖𝑡𝑖𝑣𝑖𝑡𝑦 𝑅𝑎𝑛𝑔𝑒
to be ± 0. 9℃ (Ito et al., 2004).
𝑇 = 𝑇 + 𝑇 (A6)
𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑎𝑥−𝑂𝑝𝑡 𝑜𝑝𝑡𝑖𝑚𝑎𝑙  𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑒𝑛𝑠𝑖𝑡𝑖𝑣𝑖𝑡𝑦 𝑅𝑎𝑛𝑔𝑒
𝑇 = 𝑇 − 𝑇 (A7)
𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑖𝑛−𝑂𝑝𝑡 𝑜𝑝𝑡𝑖𝑚𝑎𝑙  𝑆𝑝𝑎𝑑𝑖𝑥 𝑆𝑒𝑛𝑠𝑖𝑡𝑖𝑣𝑖𝑡𝑦 𝑅𝑎𝑛𝑔𝑒
𝑇  =  296. 78 +  0. 9 =  297. 68 (A8)
𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑎𝑥−𝑂𝑝𝑡
26 of 32
𝑇  =  296. 78 −  0. 9 =  295. 88 (A9)
𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑖𝑛−𝑂𝑝𝑡
One error in this approximation compared to real data is that stabilization temperatures will not
be correlated with outside air temperature.
We needed a function for ∂𝑄∂𝑡 that within the spadix sensitivity range would allow the change in
generation to stay very close to zero. Then, if temperature exits the minimum or maximum
temperature, the generation rate needs to increase or decrease in an approximately linear fashion
as it did in the previous model. To match these parameters a piecewise function was designed
with a cubic function defined around the spadix sensitivity range and two linear functions at
higher or lower temperatures. Our cubic function was defined by spadix optimal temperature (
𝑇 ) and a sensitivity scaling constant, 𝑘 .
𝑜𝑝𝑡𝑖𝑚𝑎𝑙 2
∂𝑄
∂𝑡 =  𝑘 * (𝑇 − 𝑇 (𝑡))
3 (A10)
2 𝑜𝑝𝑡𝑖𝑚𝑎𝑙 𝑆𝑝𝑎𝑑𝑖𝑥
For temperatures around spadix sensitivity range
We wanted a 𝑘 that allowed the part of the cubic function that has a nearly zero slope to span
2
the range of temperatures that skunk cabbage allows its spadix to be at. To accomplish this, 𝑘
2
was set to be 0.000002.
The two linear functions for ∂𝑄∂𝑡 at high and low temperatures were based on the initial model we
developed but centered around 𝑇 and 𝑇 respectively instead of
𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑎𝑥−𝑂𝑝𝑡 𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑖𝑛−𝑂𝑝𝑡
𝑇 .
𝑜𝑝𝑡𝑖𝑚𝑎𝑙
∂𝑄
∂𝑡  =  𝑘 * (𝑇 − 𝑇 (𝑡)): For High Temperatures (A11)1 𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑎𝑥−𝑂𝑝𝑡 𝑆𝑝𝑎𝑑𝑖𝑥
∂𝑄
∂𝑡  =  𝑘 * (𝑇 − 𝑇 (𝑡)): For low Temperatures (A12)1 𝑆𝑝𝑎𝑑𝑖𝑥−𝑀𝑖𝑛−𝑂𝑝𝑡 𝑆𝑝𝑎𝑑𝑖𝑥
Using an equation solver we are able to determine the two intersection points of the three
functions listed above defined as 𝑇 and 𝑇 respectively. These
𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡 𝐸𝑄10𝑎−11𝑎 𝑖𝑛𝑡𝑒𝑟𝑠𝑒𝑐𝑡 𝐸𝑄10𝑎−12𝑎
values were used as the bounds of the piecewise function in order to ensure a continuous
function.
27 of 32
Applying known values, the final function for ∂𝑄∂𝑡 can be expressed as:
9.2 Model Parameters
9.2.1 Outer convective heat transfer coefficient
The value for convective heat transfer coefficient (h) on the outside boundary of the spathe was
approximated using the following equation for natural convection over a vertical plane. The
shape of the spathe is approximately cylindrical, and since equation 1 is satisfied, the cylinder
can be further approximated as a vertical plane where
𝐷/𝐿 >  35/𝐺𝑟1/4 (A15)
𝐿
𝑁𝑢 = (0. 825 + ((0. 387𝑅𝑎1/6)/[1 + (0. 492/𝑃𝑟)9/16]8/27)2 (A16)
𝐿 𝐿
𝑁𝑢 =  ℎ𝐿/𝑘 (A18)
𝐿 𝑎𝑖𝑟
Where D is the diameter of the spathe, L is the height of the spathe, GrL is the Grashof number,
RaL is the Rayleigh number for height L, Pr is the Prandtl number, h is the convective heat
transfer coefficient, and kair is the thermal conductivity of the outside air. The values of these
variables change based on air temperature (Tair) and the temperature difference between the
spathe surface and outside air (ΔT). Table A1 displays h values calculated at the extremes of Tair
and ΔT commonly experienced by the plant as concluded from experimental data (Seymour &
Blaylock, 1999).
Table A1: Convective heat transfer coefficient [W/m2*K] for the outside of the spathe in
windless conditions
Tair=15℃ Tair=-10℃
ΔT=2℃ 3.33 3.40
ΔT=10℃ 4.91 5.05
Based on the data in Table A1, we concluded that, for a wind free environment, a constant h =
4.2 is accurate for our model. This h varies based on wind conditions, but for the purposes of
comparison to the wind-free experimental data available, calculating this h value was essential.
28 of 32
9.2.2 Parameters list
Table A2: List of parameter values with sources.
Parameter name Value Source
Stefan-Bolzmann constant 5.6704*10-8 [W/(m2*K)] Bionumbers: ID 101926
View factor spadix-spathe 1 The spathe nearly
encapsulates the spadix
meaning view factor can be
simplified
Spadix height 13.9 [mm] Seymour & Blaylock, 1999
Spadix diameter 10.5 [mm] Seymour & Blaylock, 1999
Spadix Volume 1.2036 [m3] Calculated from geometry
Spadix exposed surface area 545.105 [mm] Calculator from spadix
geometry
Spathe height 58.1 [mm] Seymour & Blaylock, 1999
Spathe diameter 26.05 [mm] Seymour & Blaylock, 1999
Spadix thermal conductivity 0.3 [W/(m*K)] Takahashi et al., 2007
Spadix density 1690 [kg/m3] Takahashi et al., 2007
Spadix heat capacity 1660 [J/kg*K] Takahashi et al., 2007
Spathe thermal conductivity 0.3 [W/(m*K)] Assumed to be equal to
spadix
Spathe density 1690 [kg/m3] Assumed to be equal to
spadix
Spathe heat capacity 1660 [J/kg*K] Assumed to be equal to
spadix
Convective heat transfer 4.2 [W/(m2*K)] Calculated as described in
coefficient 9.2.1
29 of 32
9.2 Computational Costs
The model run performed for the results section of this report took 220 seconds to run and used
1.32 GB of physical memory and 1.62 GB of virtual memory. This short computation time
allowed for many repeated runs during development and sensitivity analysis. The largest factors
in achieving such an efficient model were the mesh and the time step. The mesh was set to be as
large as possible without causing error which allowed for a computation with only 15831 degrees
of freedom.
30 of 32
10.0 References
Ito, K., Ito, T., Onda, Y., & Uemura, M. (2004). Temperature-triggered periodical thermogenic
oscillations in skunk cabbage (Symplocarpus foetidus). Plant and Cell Physiology, 45(3),
257–264. https://doi.org/10.1093/pcp/pch038
Ito, K., Onda, Y., Sato, T., Abe, Y., & Uemura, M. (2003). Structural requirements for the
perception of ambient temperature signals in homeothermic heat production of skunk
cabbage (Symlocarpus foetidus). Plant, Cell & Environment, 26(6), 783–788.
https://doi.org/10.1046/j.1365-3040.2003.00989.x
Ito-Inaba, Y., Hida, Y., Ichikawa, M., Kato, Y., & Yamashita, T. (2008). Characterization of the
plant uncoupling protein, SrUCPA, expressed in spadix mitochondria of the thermogenic
skunk cabbage. Journal of Experimental Botany, 59(4), 995–1005.
https://doi.org/10.1093/jxb/ern024
Jayalakshmy, M. S., & Philip, J. (2010). Thermophysical properties of plant leaves and their
influence on the environment temperature. International Journal of Thermophysics, 31,
2295–2304.
Seymour, R. S. (2004). Dynamics and precision of thermoregulatory responses of eastern skunk
cabbage Symplocarpus foetidus. Plant, Cell & Environment, 27(8), 1014–1022.
https://doi.org/10.1111/j.1365-3040.2004.01206.x
Seymour, R. S., & Blaylock, A. J. (1999). Switching off the heater: Influence of ambient
temperature on thermoregulation by eastern skunk cabbage Symplocarpus foetidus.
Journal of Experimental Botany, 50(338), 1525–1532.
https://doi.org/10.1093/jexbot/50.338.1525
Seymour, R. S., Ito, Y., Onda, Y., & Ito, K. (2009). Effects of floral thermogenesis on pollen
31 of 32
function in Asian skunk cabbage Symplocarpus renifolius. Biology Letters, 5(4),
568–570. https://doi.org/10.1098/rsbl.2009.0064
Shamshiri, M., Jafari, R., & Momen, G. (2021). Potential use of smart coatings for icephobic
applications: A review. Surface and Coatings Technology, 424, 127656.
https://doi.org/10.1016/j.surfcoat.2021.127656
Smith, A. (2022, January 24). 2021 U.S. billion-dollar weather and climate disasters in
historical context. Climate.gov.
http://www.climate.gov/news-features/blogs/beyond-data/2021-us-billion-dollar-weather-
and-climate-disasters-historical
Takahashi, K., Ito, T., Onda, Y., Endo, T., Chiba, S., Ito, K., & Osada, H. (2007). Modeling of the
thermoregulation system in the skunk cabbage: Symplocarpus foetidus. Physical Review
E, 76(3), 031918. https://doi.org/10.1103/PhysRevE.76.031918
32 of 32