This tutorial covers how to construct energy balances for applications such as estimating the body temperature of organisms. The tutorial first reviews general functions to estimate forms of heat exchange and balance the exchanges. Finally, we present functions to implement specialized energy balances for particular organisms.
Solar radiation, Qabs
The solar radiation absorbed by animals, \(Q_{abs} (W)\), is the sum of direct \(S_{dir}\), diffuse \(S_{dif}\), and reflected \(S_{ref}\) solar radiation (\(W/m^2\)). The sum is weighted by the organism’s surface area \(A (m)\) exposed to the radiation sources. Additionally, all forms of incoming radiation are multiplied by the solar absorptivity of the animal surface (\(a\) proportion) to estimate absorbed radiation. The summation of incoming solar radiation is thus as follows:
\[Qabs= a \times A_{dir} \times S_{dir} + a \times A_{dif} \times S_{dif} + a \times A_{ref} \times S_{ref},\]
where \(A_{dir}\), \(A_{dif}\), and \(A_{ref}\) are the surface areas exposed to direct, diffuse, and reflected solar radiation, respectively.
The summation is available in R as follows:
plot(x = seq(100, 1000, 100),
y = Qradiation_absorbed(a = 0.9,
A = 1,
psa_dir = 0.4,
psa_dif = 0.4,
psa_ref = 0.4,
S_dir = seq(100, 1000, 100),
S_dif = 200,
rho = 0.5),
type = "l",
xlab = expression("direct solar radiation" ~ (W/m^{2})),
ylab = "solar radiation absorbed (W)")
points(x = seq(100, 1000, 100),
y = Qradiation_absorbed(a = 0.9,
A = 1,
psa_dir = 0.2,
psa_dif = 0.4,
psa_ref = 0.4,
S_dir = seq(100, 1000, 100),
S_dif = 200,
rho = 0.5),
type = "l",
lty = "dashed")
points(x = seq(100, 1000, 100),
y = Qradiation_absorbed(a = 0.45,
A = 1,
psa_dir = 0.4,
psa_dif = 0.4,
psa_ref = 0.4,
S_dir = seq(100, 1000, 100),
S_dif = 200,
rho = 0.5),
type = "l",
lty = "dotted")
points(x = seq(100, 1000, 100),
y = Qradiation_absorbed(a = 0.9,
A = 1,
psa_dir = 0.4,
psa_dif = 0.4,
psa_ref = 0.4,
S_dir = seq(100, 1000, 100),
S_dif = 200,
rho = 0.2),
type = "l",
lty = "dotdash")
legend(x = "topleft",
title = "parameters",
legend = c("a = 0.9, psa_dir = 0.4, rho = 0.5", "a = 0.9, sa_dir = 0.2, rho = 0.5", "a = 0.45, psa_dir = 0.4, rho = 0.5", "a = 0.9, psa_dir = 0.4, rho = 0.2"),
lty = c("solid", "dashed", "dotted", "dotdash"))
In the functions, psa_dir, psa_dif, and psa_ref are the proportions of surface area exposed to direct, diffuse, and reflected solar radiation, respectively. The functions proportion_silhouette_area() and proportion_silhouette_area_shapes() provide assistance calculating the proportional of surface area exposed to direct radiation via estimating silhouette area. Diffuse (or scattered) solar radiation is received by the upward facing surface and solar radiation reflected off the ground is received by the downward facing surface, minus any area in contact with the ground. Additionally, see the microclimate tutorial for tools for estimating incoming solar radiation.
Convection, Qconv
Animals exchange heat with the air or water they are immersed in at a rate determined by the difference between the temperature of the animal, \(T_b (K)\), and that of the air or water, \(T_a (K)\). The rate is also determined by the animal’s surface area exposed to the air or water, which we estimate from the input parameters of the animal’s surface area, \(A (m^2)\), and the proportion of the surface area in contact with the surrounding fluid. A heat transfer coefficient, \(H_L (W m^{-2} K^{-1})\), which can be estimated using the relations below, is used to quantify the rate of heat exchange. An enhancement factor multiplier, \(ef\), can also be incorporated to account for increases in heat exchange resulting from air turbulence in field conditions. Conduction can be estimated as follows:
\[Q_{conv}= ef\cdot H_L(A\cdot \mbox{proportion})(T_a-T_b).\]
The function is available in R:
plot(x = 293:313,
y = Qconvection(T_a = 303,
T_b = 293:313,
H = 5,
A = 0.0025,
proportion = 0.6),
type = "l",
xlab = "ground temperature (K)",
ylab = "Q convection (W)")
points(x = 293:313,
y = Qconvection(T_a = 303,
T_b = 293:313,
H = 15,
A = 0.0025,
proportion = 0.6),
type = "l",
lty = "dashed")
points(x = 293:313,
y = Qconvection(T_a = 303,
T_b = 293:313,
H = 5,
A = 0.0025,
proportion = 0.3),
type = "l",
lty = "dotted")
points(x = 293:313,
y = Qconvection(T_a = 303,
T_b = 293:313,
H = 15,
A = 0.0025,
proportion = 0.3),
type = "l",
lty = "dotdash")
legend(x = "topleft",
title = "parameters",
legend = c("H = 5, proportion = 0.6", "H = 15, proportion = 0.6", "H = 5, proportion = 0.3", "H = 15, proportion = 0.3"),
lty = c("solid", "dashed", "dotted", "dotdash"))
The convective heat transfer coefficient depends on whether the heat exchange is free or forced (due to fluid flow). The variables involved in heat transfer can be combined into dimensionless groups that are used to compare the scales of heat exchange processes and thus whether modeling free or forced convection is more appropriate. We provide the primary dimensionless numbers.
The Nusselt number provides a dimensionless estimate of conductance and is estimated as \(Nu = H_L \cdot D / K\), where \(H_L\) is the convective heat transfer coefficient (\(W m^{-2} K^{-1}\)), \(D\) is the characteristic dimension (\(m\)), and \(K\) is thermal conductivity (\(W K^{-1}m^{-1}\)). The characteristic dimension quantifies the size of the organism relevant to convective heat exchange and it depends on orientation relevant to the fluid flow, usually wind. The most common approximation of the characteristic dimension is the cube root of the volume of the animal (Mitchell 1976).
The Prandtl number describes the ratio of kinematic viscosity to thermal diffusivity as \(Pr= c_p \mu /K\), where \(c_p\) is the specific heat at constant pressure (\(J mol^{-1} K^{-1}\)) and \(\mu\) is dynamic viscosity (\(mol \:s^{-1}m^{-1}\)).
The Reynolds Number indicates the amount of turbulence in air or water flows and is estimated as the ratio of internal viscous forces: \(Re = uD / \nu\), where \(u\) is wind speed (\(m/s\)), \(D\) is the characteristic dimension (\(m\)), and \(\nu\) is the kinematic viscosity (\(m^2 s^{-1}\)). Kinematic viscosity is the ratio of dynamic viscosity to the density of the fluid.
The Grashof number describes the ability of a parcel of fluid warmer or colder than the surrounding fluid to rise against or fall with the attractive force of gravity. It is estimated as the ratio of a buoyant force times an inertial force to the square of a viscous force: \(Gr = g D^3 \frac{\mid Tg-Ta \mid}{Ta {\nu}^2}\), where \(g\) is gravitational acceleration (\(g= 9.8 m/s\)).
The dimensionless groups are available in R as follows:
plot(x = 1:30,
y = Nusselt_number(H_L = 1:30,
D = 0.01,
K = 0.5),
type = "l",
xlab = expression("heat transfer coefficient" ~ (W ~ m^{-2} ~ K^{-1})),
ylab = "Nusselt number")
points(x = 1:30,
y = Nusselt_number(H_L = 1:30,
D = 0.01,
K = 1.5),
type = "l",
lty = "dashed")
points(x = 1:30,
y = Nusselt_number(H_L = 1:30,
D = 0.05,
K = 0.5),
type = "l",
lty = "dotted")
points(x = 1:30,
y = Nusselt_number(H_L = 1:30,
D = 0.05,
K = 1.5),
type = "l",
lty = "dotdash")
legend(x = "topleft",
title = "parameters",
legend = c("D = 0.01, K = 0.5", "D = 0.01, K = 1.5", "D = 0.05, K = 0.5", "D = 0.05, K = 1.5"),
lty = c("solid", "dashed", "dotted", "dotdash"))
plot(x = 1:30,
y = Prandtl_number(c_p = 1:30,
mu = 0.00001,
K = 0.5),
type = "l",
xlab = expression("specific heat at constant pressure" ~ (J ~ mol^{-1} ~ K^{-1})),
ylab = "Prandtl number")
points(x = 1:30,
y = Prandtl_number(c_p = 1:30,
mu = 0.00002,
K = 0.5),
type = "l",
lty = "dashed")
points(x = 1:30,
y = Prandtl_number(c_p = 1:30,
mu = 0.00001,
K = 1.5),
type = "l",
lty = "dotted")
points(x = 1:30,
y = Prandtl_number(c_p = 1:30,
mu = 0.00002,
K = 1.5),
type = "l",
lty = "dotdash")
legend(x = "topleft",
title = "parameters",
legend = c("mu = 0.00001, K = 0.5", "mu = 0.00002, K = 0.5", "mu = 0.00001, K = 1.5", "mu = 0.00002, K = 1.5"),
lty = c("solid", "dashed", "dotted", "dotdash"))
plot(x = seq(0, 5, 0.2),
y = Reynolds_number(u = seq(0, 5, 0.2),
D = 0.001,
nu = 1),
type = "l",
xlab = "wind speed (m/s)",
ylab = "Reynolds number")
points(x = seq(0, 5, 0.2),
y = Reynolds_number(u = seq(0, 5, 0.2),
D = 0.001,
nu = 1.5),
type = "l",
lty = "dashed")
points(x = seq(0, 5, 0.2),
y = Reynolds_number(u = seq(0, 5, 0.2),
D = 0.002,
nu = 1),
type = "l",
lty = "dotted")
points(x = seq(0, 5, 0.2),
y = Reynolds_number(u = seq(0, 5, 0.2),
D = 0.002,
nu = 1.5),
type = "l",
lty = "dotdash")
legend(x = "topleft",
title = "parameters",
legend = c("D = 0.001, nu = 1", "D = 0.001, nu = 1.5", "D = 0.002, nu = 1", "D = 0.002, nu = 1.5"),
lty = c("solid", "dashed", "dotted", "dotdash"))
plot(x = seq(0, 0.01, 0.001),
y = Grashof_number(T_a = 30,
T_g = 35,
D = seq(0, 0.01, 0.001),
nu = 1.0),
type = "l",
xlab = "characteristic dimension (m)",
ylab = "Grashof number")
points(x = seq(0, 0.01, 0.001),
y = Grashof_number_Gates(T_a = 30,
T_g = 35,
beta = 0.00367,
D = seq(0, 0.01, 0.001),
nu = 1.0),
type = "l",
col = "red")
points(x = seq(0, 0.01, 0.001),
y = Grashof_number(T_a = 30,
T_g = 35,
D = seq(0, 0.01, 0.001),
nu = 1.4),
type = "l",
lty = "dashed")
points(x = seq(0, 0.01, 0.001),
y = Grashof_number_Gates(T_a = 30,
T_g = 35,
beta = 0.00367,
D = seq(0, 0.01, 0.001),
nu = 1.4),
type = "l",
col = "red",
lty = "dashed")
legend(x = "topleft",
title = "parameters",
legend = c("Grashof_number, nu = 1", "Grashof_number_Gates, nu = 1", "Grashof_number, nu = 1.4", "Grashof_number_Gates, nu = 1.4"),
lty = c("solid", "dashed", "dotted", "dotdash"))
Relations among the dimensionless groups can also be used for estimation. Empirically-derived relationships can be used to estimate the Nusselt number \(Nu\) from the Reynolds number \(Re\) or the Grashof number \(Gr\). Comparisons of \(Gr\) and \(Re\) allow determining whether modeling free or forced convection is appropriate. Forced convection is appropriate when \(Gr \leq 16 Re^2\). For convenience, we provide the relationships as R functions:
plot(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "sphere"),
type = "l",
xlab = "Re",
ylab = "Nu",
ylim = c(0,4))
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "cylinder"),
type = "l",
lty = "dashed")
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "frog"),
type = "l",
lty = "dotted")
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "lizard_traverse_to_air_flow"),
type = "l",
lty = "dotdash")
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "lizard_parallel_to_air_flow"),
type = "l",
lty = "dotdash",
col = "red")
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "lizard_surface"),
type = "l",
lty = "dotdash",
col = "orange")
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "lizard_elevated"),
type = "l",
lty = "dotdash",
col = "green")
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "flyinginsect"),
type = "l",
lty = "longdash")
points(x = 1:10,
y = Nusselt_from_Reynolds(Re = 1:10,
taxon = "spider"),
type = "l",
lty = "twodash")
legend(x = "topright",
title = "taxa",
legend = c("sphere", "cylinder", "frog", "lizard_traverse_to_air_flow", "lizard_parallel_to_air_flow", "lizard_surface", "lizard_elevated", "flyinginsect", "spider"),
lty = c("solid", "dashed", "dotted", "dotdash", "dotdash", "dotdash", "dotdash", "longdash", "twodash"),
col = c("black", "black", "black", "black", "red", "orange", "green", "black", "black" ))
Nusselt_from_Grashof(Gr = 5)
## [1] 0.7177674
free_or_forced_convection(Gr = 100,
Re = 5)
## [1] "intermediate condition, mixed convection based on Nusselt numbers is appropriate"
We offer methods to estimate the convective heat transfer coefficient based on either empirical measurements (heat_transfer_coefficient()) or approximating the animal shape as a sphere (heat_transfer_coefficient_approximation()). Approximating the animal shape as a sphere enables simplification while also producing a reasonable approximation (Mitchell 1976). The functions approximate forced convective heat transfer as a function of windspeed \(u (m/s)\), the characteristic dimension \(D (m)\), the thermal conductivity of the air \(K (W m^{-1} K^{-1})\), the kinematic viscosity of the air \(\nu (m^2 s^{-1})\), and the taxa or a generic shape. An additional, simplified function (heat_transfer_coefficient_simple()) provides a reasonable approximation based on \(u\) and \(D\) for most environmental conditions.
oldpar <- par()
par(mar = c(5, 5, 3, 2))
plot(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "sphere"),
type = "l",
xlab = "Air velocity (m/s)",
ylab = expression("heat transfer coefficient," ~ H[L] ~ (W ~ m^{-2} ~ K^{-1})))
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "cylinder"),
type = "l",
lty = "dashed")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "frog"),
type = "l",
lty = "dotted")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "lizard_surface"),
type = "l",
lty = "dotdash")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "lizard_elevated"),
type = "l",
lty = "dotdash")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "flyinginsect"),
type = "l",
lty = "longdash")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "spider"),
type = "l",
lty = "twodash")
legend(x = "bottomright",
title = "taxa",
legend = c("sphere", "cylinder", "frog", "lizard_surface", "lizard_elevated", "flyinginsect", "spider"),
lty = c("solid", "dashed", "dotted", "dotdash", "dotdash", "dotdash", "dotdash", "longdash", "twodash"))
par(mar = c(5, 5, 3, 2))
plot(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_approximation(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "sphere"),
type = "l",
xlab = "Air velocity (m/s)",
ylab = expression("heat transfer coefficient," ~ H[L] ~ (W ~ m^{-2} ~ K^{-1})))
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_approximation(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "frog"),
type = "l",
lty = "dashed")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_approximation(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "lizard"),
type = "l",
lty = "dotted")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_approximation(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "flyinginsect"),
type = "l",
lty = "dotdash")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_approximation(u = seq(0, 3, 0.25),
D = 0.05,
K = 25.7 * 10^(-3),
nu = 15.3 * 10^(-6),
taxon = "spider"),
type = "l",
lty = "longdash")
legend(x = "bottomright",
title = "taxa",
legend = c("sphere", "frog", "lizard", "flyinginsect", "spider"),
lty = c("solid", "dashed", "dotted", "dotdash", "longdash"))
plot(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_simple(u = seq(0, 3, 0.25),
D = 0.05,
type = "Spotila"),
type = "l",
xlab = "Air velocity (m/s)",
ylab = expression("heat transfer coefficient," ~ H[L] ~ (W ~ m^{-2} ~ K^{-1})))
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_simple(u = seq(0, 3, 0.25),
D = 0.03,
type = "Spotila"),
type = "l",
lty = "dashed")
points(x = seq(0, 3, 0.25),
y = heat_transfer_coefficient_simple(u = seq(0, 3, 0.25),
D = 0.07,
type = "Spotila"),
type = "l",
lty = "dotted")
legend(x = "bottomright",
title = "characteristic dimension (m)",
legend = c(0.03, 0.05, 0.07),
lty = c("dashed", "solid", "dotted"))
Conduction, Qcond
Animals exchange heat with the substrate and other surfaces via contact proportional to the difference between the surface temperature of the animal, \(T_b (K)\), and the substrate temperature, \(T_s (K)\). This conduction occurs via diffusion of heat. The extent of energy exchange is additionally determined by the area of contact and the thickness of the skin or other animal surface (assumes a well-mixed interior, see NicheMapR for approaches without this assumption). Finally, the rate of diffusion depends on the thermal conductivity of the materials in contact. Our two basic functions for estimating conduction differ in whether the rate-limiting step is the thermal conductivity of the animal [conductance_animal()] or the substrate [conductance_substrate()]. The thermal conductivity of the ground is generally greater than that of animal tissues, so animal thermal conductivity is generally the rate-limiting step and most applications should use conductance_animal().
Both functions take as parameters the animal’s surface area, \(A (m^2)\), and the proportion of the surface area in contact with the substrate, \(\mbox{proportion}\). When animal conductance is the rate limiting step, \(Q_{cond}\) can be estimated as follows:
\[Q_{cond}= A\cdot \mbox{proportion} \cdot K(T_b-T_g)/d, \]
where \(K\) is thermal conductivity (\(W K^{-1} m^{-1}\)), \(T_g\) is ground (surface) temperature \((K)\), \(T_b\) is body temperature \((K)\), and \(d\) is the mean thickness of the animal skin (surface, \(m\)). However, the extremely thin skins of most animals will result in high levels of conduction. To estimate steady state conditions, one viable approach is to assume that heat is exchanged down to a given soil levels such that \(d=0.25m\). This formulation assumes the organism has a well-mixed interior rather than an interior temperature gradient.
When substrate thermal conductivity is the rate-limiting step, \(Q_{cond}\) can be estimated as follows:
\[Q_{cond} = A\cdot \mbox{proportion} \cdot (2K_g/D)(T_b-T_g),\]
where \(K_g\) is the thermal conductivity of substrate (\(W K^{-1} m^{-1}\)) and \(D\) is the characteristic dimension of the animal (\(m\)).The most common approximation of the characteristic dimension is the cube root of the volume of the animal (Mitchell 1976).
The functions are available in R:
plot(x = 293:313,
y = Qconduction_animal(T_g = 293:313,
T_b = 303,
d = 0.025,
K = 0.3,
A = 10^-3,
proportion = 0.2),
type = "l",
xlab = "ground temperature (Tb, K)", ylab = "Q conductance (W)")
points(x = 293:313,
y = Qconduction_animal(T_g = 293:313,
T_b = 303,
d = 0.025,
K = 0.5,
A = 10^-3,
proportion = 0.2),
type = "l",
lty = "dashed")
points(x = 293:313,
y = Qconduction_animal(T_g = 293:313,
T_b = 303,
d = 0.025,
K = 0.3,
A = 10^-3,
proportion = 0.4),
type = "l",
lty = "dotted")
points(x = 293:313,
y = Qconduction_animal(T_g = 293:313,
T_b = 303,
d = 0.025,
K = 0.5,
A = 10^-3,
proportion = 0.4),
type = "l",
lty = "dotdash")
legend(x = "topright",
title = "parameters",
legend = c("K = 0.3, proportion = 0.2", "K = 0.5, proportion = 0.2", "K = 0.3, proportion = 0.4", "K = 0.5, proportion = 0.4"),
lty = c("solid", "dashed", "dotted", "dotdash"))
plot(x = 293:313,
y = Qconduction_substrate(T_g = 293:313,
T_b = 303,
D = 0.01,
K_g = 0.3,
A = 10^-2,
proportion = 0.2),
type = "l",
xlab = "ground temperature (Tb, K)",
ylab = "Q conductance (W)")
points(x = 293:313,
y = Qconduction_substrate(T_g = 293:313,
T_b = 303,
D = 0.01,
K_g = 0.5,
A = 10^-2,
proportion = 0.2),
type = "l",
lty = "dashed")
points(x = 293:313,
y = Qconduction_substrate(T_g = 293:313,
T_b = 303,
D = 0.01,
K_g = 0.3,
A = 10^-2,
proportion = 0.4),
type = "l",
lty = "dotted")
points(x = 293:313,
y = Qconduction_substrate(T_g = 293:313,
T_b = 303,
D = 0.01,
K_g = 0.5,
A = 10^-2,
proportion = 0.4),
type = "l",
lty = "dotdash")
legend(x = "topright",
title = "parameters",
legend = c("K = 0.3, proportion = 0.2", "K = 0.5, proportion = 0.2", "K = 0.3, proportion = 0.4", "K = 0.5, proportion = 0.4"),
lty = c("solid", "dashed", "dotted", "dotdash"))
Evaporation, Qevap
We provide functions to estimate heat loss associated with evaporation. We offer an empirically-derived relationship for a lizard (Porter et al. 1973; in Gates 2012) and a function based on resistances for amphibians (Spotila, O’Connor, and Bakken 1992). The function for amphibians is based on the latent heat of vaporization of water: evaporation of water results in heat loss at a rate of \(2.44 \times 10^6 J kg^{-1}\) at biological temperatures. The rate of water loss for an amphibian with fully wetted skin can be expressed as
\[ E_c= A (\rho_s-h \rho_a)/r_e,\]
where \(A\) is surface area (\(m^2\)), \(r_e\) is the external (convective) resistance to water vapor transport (\(s/m\)), \(\rho_s\) is the saturation water vapor density at skin surface (\(kg/m^3\)), \(h\) is relative humidity (between 0 and 1), and \(\rho_a\) is the saturation water vapor density in ambient air (\(kg/m^3\)).
The rate of water loss for an amphibian without fully wetted skin can be expressed as \[E_c= A (\rho_s-h\rho_a)/(r_i+r_e),\] where \(r_i\) is the internal (cutaneous) resistance to vapor transport (\(s/m\)). We provide R functions for heat loss through evaporation as the product of \(E_c\) and the latent heat of vaporization of water:
# kPa
vp <- saturation_vapor_pressure(293:313)
# convert to kg/m^3
e_s <- vp * 0.032
# kPa
vp <- saturation_vapor_pressure(293:313-10)
# convert to kg/m^3
e_a <- vp * 0.032
temps <- 293:313
Qevaps <- rep(NA,21)
Qevaps_wet <- rep(NA,21)
for (ind in 1:21) {
Qevaps[ind] <- Qevaporation(A = 0.1,
T_b = temps[ind],
taxon = "amphibian",
e_s = e_s[ind],
e_a = e_a[ind],
hp = 0.5,
H = 20,
r_i = 50)
Qevaps_wet[ind] <- Qevaporation(A = 0.1,
T_b = temps[ind],
taxon = "amphibian_wetskin",
e_s = e_s[ind],
e_a = e_a[ind],
hp = 0.5,
H = 20,
r_i = 50)
}
plot(x = temps,
y = Qevaps,
type = "l",
xlab = "body temperature (K)",
ylab = "evaporative heat loss (W)",
lty = "dashed",
ylim = c(0, 400))
points(x = temps,
y = Qevaps_wet,
type = "l",
lty = "dotted")
Qevap <- unlist(lapply(293:313,
FUN = Qevaporation,
A = 0.1,
taxon = "lizard"))
points(x = 293:313,
y = Qevap,
type = "l" )
legend(x = "right",
title = "taxa",
legend = c("amphibian", "amphibian_wetskin", "lizard"),
lty = c("dashed", "dotted", "solid"))
We estimate external resistance to water vapor transport, \(r_e\), using the Lewis rule describing the relationships among the coefficients for heat and mass transport: \(r_e= 0.93 \rho c_p/H_L\), where \(\rho\) is the density of air (\(kg m^{-3}\)), \(c_p\) is the specific heat of air at constant pressure (\(J kg^{-1}K^{-1}\)), and \(H_L\) is the convective heat transfer coefficient (\(W m^{-2} K^{-1}\)). We approximate \(\rho c_p\)as 1,200 \(J m^{-3}KC^{-1}\) as commonly done in biological applications (Spotila, O’Connor, and Bakken 1992). We provide the relationship for \(r_e\) as a function:
par(mar = c(5, 5, 3, 2))
plot(x = 10:30,
y = external_resistance_to_water_vapor_transfer(H = 10:30),
type = "l",
xlab = expression("heat transfer coefficient" ~ (W ~ m^{-2} ~ K^{-1})),
ylab = expression("external resistance, " ~ r[e] ~ (sm^{-1})))
At constant temperatures, the relationship between vapor pressure and vapor density is linear and vapor pressure can be approximated at temperatures between \(0°C\) and \(40°C\) as \(e_s= 10^{0.02604 T_a+2.82488}\), where \(T_a\) is air temperature (\(°C\)). We provide the approximation as a function:
plot(x = 10:30,
y = saturation_water_vapor_pressure(T_a = 10:30),
type = "l",
xlab = "temperature (°C)",
ylab = "saturation water vapor pressure (Pa)")
Estimating body temperatures, Tb
TrenchR primarily includes models for estimating the body temperatures (generally the operative temperatures, \(T_e\)) of organisms that have reached an equilibrium with their environment (“steady-state,” no heat retention). The functions for estimating body temperatures all require input temperatures in degrees Celsius and they output temperatures in degrees Celsius for accessibility, but our explanations below are based on Kelvin since the unit is primarily used for heat balance calculations. We use K subscripts to specify the use of Kelvin.
The following function uses the relationships above to solve the energy balance for \(T_b (°C)\).
t.seq <- lapply(20:40,
FUN = Tb_Gates,
A = 0.1,
D = 0.025,
psa_dir = 0.6,
psa_ref = 0.4,
psa_air = 0.6,
psa_g = 0.0,
T_g = 30,
Qabs = 10,
epsilon = 0.95,
H_L = 10,
ef = 1.3,
K = 0.5)
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature, Ta (°C)",
ylab = "body temperature, Tb (°C)",
xlim = c(22, 42),
ylim = c(22, 42))
abline(a = 0,
b = 1,
col = "gray")
t.seq <- lapply(20:40,
FUN = Tb_Gates,
A = 0.1,
D = 0.025,
psa_dir = 0.6,
psa_ref = 0.4,
psa_air = 0.6,
psa_g = 0.0,
T_g = 30,
Qabs = 0,
epsilon = 0.95,
H_L = 10,
ef = 1.3,
K = 0.5)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tb_Gates,
A = 0.1,
D = 0.025,
psa_dir = 0.6,
psa_ref = 0.4,
psa_air = 0.6,
psa_g = 0.0,
T_g = 30,
Qabs = 20,
epsilon = 0.95,
H_L = 10,
ef = 1.3,
K = 0.5)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
legend(x = "bottomright",
title = "Radiation, Qabs",
legend = c(0, 10, 20),
lty = c("dashed", "solid", "dotted"))
Campbell and Norman (2012) use a somewhat simplified energy balance to express \(T_{bK}\) as a function of \(T_{aK}\) plus or minus a temperature increment determined by absorbed radiation, wind speed, and animal morphology:
\[ T_{bK}=T_{aK}+(S_{abs}-Q_{emit})/(c_p(g_r+g_{Ha})),\]
where \(T_{aK}\) is air temperature \((K)\), \(S_{abs}\) is the solar and thermal radiation absorbed (\(W m^{-2}\)), \(Qemit\) is emitted thermal radiation (\(W m^{-2}\)), \(c_p\) is the specific heat of air (\(J mol^{-1} K^{-1}\)), \(g_{r}\) is radiative conductance \((mol\:m^{-2} s^{-1})\), and \(g_{Ha}\) is the boundary conductance \((mol\:m^{-2} s^{-1})\). The model is based on estimating \(T_{bK}\) for a blackbody (perfectly absorbing) cavity with air temperatures equal to surface temperatures (\(T_{aK} = T_{sK}\)). The model assumes that the cavity provides the same heat load as the natural environment and thus equates metabolic heat production and evaporative cooling in the two environments. The formulation assumes a well-mixed interior of the animal and also omits conduction with the substrate. In this scenario, organisms emit thermal radiation from their surface proportional to the fourth power of \(T_{aK}\):
\[Q_{emit}= \epsilon \sigma T_{aK}^4 ,\]
where \(\epsilon\) is the proportional longwave infrared emissivity of skin [0.95-1 for most animals, (Gates 2012)] and \(\sigma\) is the Stefan-Boltzmann constant (\(5.673*10^{-8} W m^{-2} K^{-4}\)).
The thermal radiation exchanged between the animals and the walls of the cavity is proportional to the temperature differences and the radiative conductance. Campbell and Norman (2012) thus use a denominator term that combines thermal radiative exchange with convective heat exchange. The radiative conductance describing the heat exchange between the core of the animal and the environment is estimated as \(g_r= 4 \sigma T_{aK}^3/c_p(mol\:m^{-2} s^{-1})\). The boundary conductance (heat exchange with the air via convection) is estimated assuming forced convection driven by naturally turbulent wind and an empirical approximation of animal shapes (Mitchell 1976):
\[g_{Ha}(mol \:m^{-2} s^{-1})= 1.4 \times 0.135 \sqrt{(u/D)},\]
where 1.4 is a factor to account for increased convection in natural environments, \(u\) is wind speed (\(m/s\)), and \(D\) is the characteristic dimension of the animal (\(m\)).
The function to estimate \(T_b\) (often referred to as \(T_e\)) is available in R as follows (including the above calculations of \(Q_{emit}\), \(g_r\), and \(g_{Ha}\)):
plot(x = 20:40,
y = Tb_CampbellNorman(T_a = 20:40,
T_g = 30,
S = 600,
a_s = 0.7,
a_l = 0.96,
epsilon = 0.96,
c_p = 29.3,
D = 0.17,
u = 1),
type = "l",
xlab = "air temperature (C)",
ylab = "body temperature (C)",
xlim = c(23,37),
ylim = c(17,67))
abline(a = 0,
b = 1,
col = "gray")
points(x = 20:40,
y = Tb_CampbellNorman(T_a = 20:40,
T_g = 30,
S = 200,
a_s = 0.7,
a_l = 0.96,
epsilon = 0.96,
c_p = 29.3,
D = 0.17,
u = 1),
type = "l",
lty = "dashed")
points(x = 20:40,
y = Tb_CampbellNorman(T_a = 20:40,
T_g = 30,
S = 400,
a_s = 0.7,
a_l = 0.96,
epsilon = 0.96,
c_p = 29.3,
D = 0.17,
u = 1),
type = "l",
lty = "dotted")
legend(x = "bottomright",
title = expression("Radiation, Sabs" ~ (W ~ m^{-2})),
legend = c(200,400,600),
lty = c("dashed", "dotted", "solid"))
Specialized energy balance models and functions
In addition to providing general functions that can be used to build heat budget models predicting the body temperatures for a variety of organisms, we implement specialized operative temperature models that have been developed for a variety of organisms.
Lizards
We adapted a biophysical model (Tb_lizard()) to estimate body temperatures, \(T_b (°C)\), for lizards from Buckley (2007). The function follows the modeling approach contained in Campbell and Norman (2012).
Day Length and Sun Angles
The function constrains the lizard’s foraging to hours of daylight and uses the relations below to calculate day length and the time of solar noon. The function calculates the solar declination, \(\delta\) (rad), the angular distance of the sun north or south of the earth’s equator:
\[\delta = \arcsin[0.39795 \cos(0.21631+2 \arctan (0.967 \tan [0.0086(-186+J)]))] \]
where \(J\) is the calendar day, with \(J= 1\) on January 1. One can then calculate day length, \(h_{day} (h)\), using the CBM model (Forsythe et al. 1995): where \(\phi\) is latitude (rad).
The time of solar noon (\(h\)) is calculated as \(t_0= 12-LC-ET\) where \(LC\) is the longitude correction (h) and \(ET\) is the equation of time (h). The \(LC\) is a +1/15-h correction for each degree that a location is west of a standard meridian, which occurs at 0°, 15°, …, 345°. The ET corrects for the difference between sun time and clock time based on calendar day as follows:
\[ET=\frac{-104.7\sin f+596.2 \sin{2f} +4.3 \sin{3f} -12.7 \sin{4f} -429.3 \cos{f} -2.0 \cos{2f}+19.3 \cos{3f}}{3600},\\f=\frac{\pi}{180}(279.575+0.9856J).\]
In addition, the offset needs to be taken into consideration. Although meridians are set every 15° and the time zones are supposed to equal to the time of the meridian that is closest to them, (thus anywhere between 7.5°E to 22.5°E should equal the time at the meridian at 15°E) the actual time zones do not necessarily follow the pattern. For, example, Madrid is located at 3.7°W but the time zone in Spain is set to that of 15°E (GMT+1). This produces an error in the calculations because noon at 3.7°W would be considered equal to noon at 0° when in reality, it should be 11 am at 0°. To eliminate this error, an offset value can be input.
The zenith angle, \(\psi\) (rad), is the sun angle measured from vertical, \(\cos{\psi}=\sin\delta \sin\phi + \cos\delta \cos\phi\cos\frac{\pi}{12(h-t_0)}\), where \(h\) is hour.
Azimuth angle, \(AZ\) (rad), is the angle between the south vector and the sun’s vector on the horizontal plane. The azimuth angle of the sun can be calculated from
\[cos(AZ) = \frac{-(sin\delta-cos\psi sin\phi)}{cos\phi sin\psi}.\]
Within the two values that are produced from this formula, the correct angle is determined based on the longitude of the measured place.
Radiation and Environmental Temperature
Central to the thermal influence of radiation is the Stefan-Boltzmann law, which expresses the total radiant energy over all wavelengths admitted per unit surface area of a blackbody radiator. The law yields the emitted flux density, \(B\) (\(W m^2\)), \(B= \sigma(T_{aK})^4\), where \(T_{aK}\) is air temperature (K) and \(\sigma\) is the Stefan-Boltzmann constant (\(5.67×10^{−8} W m^{−2} K^{−4}\)). Emissivity, \(\epsilon(\lambda)\), where \(\lambda\) is the wavelength, is the fraction of blackbody emittance at a given wavelength emitted by the surface of a material. Gray bodies are those with no wavelength dependence of the emissivity. Thus, emitted energy of a gray body is \(\phi (W m^{−2})= \epsilon B\). We assume lizards are gray bodies, which is reasonable for most natural surfaces. The emissivity of most natural surfaces ranges from 0.95 to 1.0. The function uses the emissivity value of 0.965 measured by Bartlett and Gates (1967) for a Sceloporus lizard. However, the emissivity of a clear atmosphere is substantially lower. We use the approximation by Swinbank (1963) to estimate clear sky emissivity, \(\epsilon_{ac}= 9.2 \times 10^{-6}(T_{aK})^2.\)
We first consider convective heat transport between the lizard’s body and the environment.
The boundary conductance of air (\(mol \: m^{−2} s^{−1}\)) is expressed as
\[g_{Ha}= 1.4 \times 0.135 \sqrt{u/d} ,\]
where \(u\) is wind velocity (\(m s^{−1}\)) and \(d\) is the characteristic dimension (\(m\)). The function uses the cube root of volume to approximate the characteristic dimension (Mitchell 1976). A factor of 1.4 is introduced to account for the increased convection due to environmental turbulence (Mitchell 1976).
We additionally consider radiative conductance, the exchange of thermal radiation between the lizard and the environment proportional to temperature differences. The radiative conductance (\(mol \:m^{−2} s^{−1}\)) is expressed as
\[g_r=\frac{4 \epsilon \sigma (T_{aK})^3}{c_p},\]
where \(c_p\) is the specific heat of air (\(29.3 J \: mol^{−1} K^{−1}\)).
We then turn to calculate the components of the radiant energy budget of a lizard. We estimate the solar (shortwave) component of this quantity by aggregating flux densities for four radiation streams: the direct irradiance on a surface perpendicular to the beam, \(S_p (W m^{−2})\); the diffuse sky irradiance on a horizontal plane, \(S_d (W m^{−2})\); the total irradiance of a horizontal surface, \(S_t (W m^{−2})\); and the reflected radiation from the ground, \(S_r (W m^{−2})\). The calculation of these flux densities requires the introduction of several additional quantities. The atmospheric transmissivity, \(\tau\), ranges between 0.6 and 0.7 for typical clear sky conditions (Gates 1980). We thus assume \(\tau\) = 0.65.
The solar constant, \(S_{p0}\), indicates extraterrestrial flux density to be \(1,360 W m^{−2}\). Optical air mass number, \(m_a\), is the ratio of slant path length through the atmosphere to zenith path length and is a function of atmospheric pressure, \(p_a (kPa)\):
\[p_a = 101.3\exp(-E/8200),\]
where \(E\) is the elevation in meters above sea level. We can then calculate \(m_a\) as \(m_a = \frac{p_a}{101.3 \cos\psi}\).
Direct irradiance is a function of the distance a solar beam travels through the atmosphere; the transmittance of the atmosphere, \(\tau\); and the incident flux density, \(S_{p0}\): \(S_p=S_{p0}\tau^{m_a}\). Sky diffuse radiation can be approximated using an empirical relation (Liu and Jordan 1960), \(S_d= 0.3(1-\tau^{m_a})S_{p0}\cos\psi\). Solar irradiance is the sum of diffuse sky irradiance and the beam irradiance on a surface: \(S_t=S_d+S_p \cos\psi\). Finally, reflected radiation is the product of albedo, \(\rho_S\), and the total shortwave irradiance: \(S_r=\rho_S S_t\). Albedo was empirically derived from satellite images and is the surface reflectance for the solar waveband.
The longwave component of a lizard’s radiant energy budget can be estimated using the Stefan-Boltzmann law. The longwave flux density from atmosphere, \(L_a (W m^{−2})\), is computed as \(L_a = \epsilon_{ac}\sigma(T_{aK})^4\), where \(\epsilon_{ac}\) is clear sky emissivity. The longwave flux density from the ground is \(L_g=\epsilon_s \sigma T_{sK}^4\), where \(T_{sK}\) is surface temperature (K) and \(\epsilon_s\) is surface emissivity.
One can then aggregate short‐ and longwave radiation to compute absorbed radiation:
\[R_{abs}+\alpha_s(F_PS_p+F_dS_d+F_rS_r)+\alpha_L(F_aL_a+F_gL_g),\]
where \(\alpha_S\) and \(\alpha_L\) are the absorptivities in the solar and thermal wavebands, respectively, and \(F_p\), \(F_d\), \(F_r\), \(F_a\), and \(F_g\) are view factors between the surface of the lizard and sources of radiation. Solar absorptivity, \(\alpha_S\), is approximately 0.9 for lizards (Gates 2012). Because absorptivity in a given waveband is equal to emissivity in that waveband, we use the thermal absorptivity, \(\alpha_L\), of 0.965, as measured by Bartlett and Gates (1967).
View factors, also known as configuration factors, refer to the fraction of radiation that is intercepted by the lizard. The beam view factor, \(F_p\), for a lizard is the ratio of the projected area perpendicular to the solar beam, \(A_p\), to the total animal area, \(A\), \(F_p=A_p/A\). We use an empirically derived relationship that was developed for a standing Callisaurus (Muth 1977). We assume a relative azimuth angle of 90°, which indicates that the long axis of the lizard is perpendicular to incoming solar radiation:
\[ A_p= \frac{[(-1.1756\times10^{-4})\psi^2+(-9.2594\times10^{-2})\psi+26.2409]A}{100}\]
Roughgarden et al. (1981) provide an expression for total area based on empirical lizard data from Norris (1967) and Porter and James (1979): \(A = 0.121m^{0.688}\) \((r^2= 0.996),\) where \(m\) is mass (\(g\)), and the relation holds for a variety of lizards from 2 to 50 g.
For diffuse shortwave and longwave radiation, the sky can be approximated as a hemisphere. The diffuse radiation view factor, \(F_d\), for a standing lizard was found to be 0.8 by Bartlett and Gates (1967). We assume that one-half of thermal radiation is received from both the ground and the sky.
The atmospheric thermal radiation factor, \(F_a\), is thus 0.5, and the ground thermal radiation factor, \(F_g\), is 0.5. We likewise assume that the reflected solar radiation view factor, \(F_r\), is 0.5. The operative environmental temperature is calculated within each grid cell as
\[T_{eK}=T_{aK}+\frac{Rabs-\epsilon_s \sigma T_{aK}^4}{c_p(g_r+g_{Ha})}\].
We implement the model in TrenchR as follows:
# sun, surface
t.seq <- lapply(20:40,
FUN = Tb_lizard,
T_g = 30,
u = 0.1,
svl = 60,
m = 10,
psi = 34,
rho_s = 0.7,
elev = 500,
doy = 200,
sun = TRUE,
surface = TRUE,
a_s = 0.9,
a_l = 0.965,
epsilon_s = 0.965,
F_d = 0.8,
F_r = 0.5,
F_a = 0.5,
F_g = 0.5)
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)",
xlim = c(20,40),
ylim = c(20,51))
abline(a = 0,
b = 1,
col = "gray")
# shade, surface
t.seq <- lapply(20:40,
FUN = Tb_lizard,
T_g = 30,
u = 0.1,
svl = 60,
m = 10,
psi = 34,
rho_s = 0.7,
elev = 500,
doy = 200,
sun = FALSE,
surface = TRUE,
a_s = 0.9,
a_l = 0.965,
epsilon_s = 0.965,
F_d = 0.8,
F_r = 0.5,
F_a = 0.5,
F_g = 0.5)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
# sun, above surface
t.seq <- lapply(20:40,
FUN = Tb_lizard,
T_g = 30,
u = 0.1,
svl = 60,
m = 10,
psi = 34,
rho_s = 0.7,
elev = 500,
doy = 200,
sun = TRUE,
surface = FALSE,
a_s = 0.9,
a_l = 0.965,
epsilon_s = 0.965,
F_d = 0.8,
F_r = 0.5,
F_a = 0.5,
F_g = 0.5)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
# shade, above surface
t.seq <- lapply(20:40,
FUN = Tb_lizard,
T_g = 30,
u = 0.1,
svl = 60,
m = 10,
psi = 34,
rho_s = 0.7,
elev = 500,
doy = 200,
sun = FALSE,
surface = FALSE,
a_s = 0.9,
a_l = 0.965,
epsilon_s = 0.965,
F_d = 0.8,
F_r = 0.5,
F_a = 0.5,
F_g = 0.5)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotdash")
legend(x = "topright",
title = "parameters",
legend = c("sun, surface", "shade, surface", "sun, above surface", "shade, above surface"),
lty = c("solid", "dashed", "dotted", "dotdash"))
There are several transient energy budget models available to predict lizard body temperatures. We include a transient model from Fei at al. (2012) that can be implemented as demonstrated below. Please refer to the initial publication for the model details. R code for another transient model of lizard body temperature is available in the following repository: https://github.com/JRubalcaba/BodyTemp.
t.seq <- lapply(20:40,
FUN = Tb_lizard_Fei,
T_g = 27,
S = 800,
lw = 30,
shade = 0.5,
m = 10.5,
Acondfact = 0.05,
Agradfact = 0.4)
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)")
abline(a = 0,
b = 1,
col = "gray")
t.seq <- lapply(20:40,
FUN = Tb_lizard_Fei,
T_g = 27,
S = 800,
lw = 30,
shade = 0.0,
m = 10.5,
Acondfact = 0.05,
Agradfact = 0.4)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tb_lizard_Fei,
T_g = 27,
S = 800,
lw = 30,
shade = 1.0,
m = 10.5,
Acondfact = 0.05,
Agradfact = 0.4)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
legend(x = "bottomright",
title = "proportion shade",
legend = c(0.0,0.5,1.0),
lty = c("dashed", "solid", "dotted"))
Butterflies
We adapted a biophysical model for Colias that was developed and field validated by Kingsolver (1983) to predict thoracic body temperature (\(T_b, °C\)) based on thermoregulatory traits (body size, basal ventral hind wing solar absorptivity, and thoracic fur thickness), behavioral posture (basking and heat-avoidance), and environmental conditions. We briefly describe the model, which was updated and detailed in Buckley and Kingsolver (2012). Adults behaviorally thermoregulate to achieve the body temperatures needed for flight and do not use endogenous heat production to elevate body temperatures. We assume that butterflies select the body temperature closest to their thermal optima (35°C) with available body temperatures bracketed by those in full sun (lateral basking posture with wings closed and the ventral hindwing surfaces oriented perpendicular to the sun) and full shade (no direct radiation). Assuming that butterflies select body temperatures bracketed by full sun and full shade eliminates the need to attempt to model microclimate availability in detail.
We describe the steady-state energy flux balance of a butterfly at rest on vegetation as \[Q_s= Q_t + Q_c,\] where \(Q_s\) is the total solar radiative heat flux (\(W\)), \(Q_t\) is the thermal radiative heat flux (\(W\)), and \(Q_c\) is the convective heat flux (\(W\)). Conduction of heat between the body and vegetation and evaporative heat loss are considered to be negligible. The solar radiative heat flux is
\[ Q_s= Q_{s,dir}+Q_{s,dif}+Q_{s,ref},\\ Q_s= \alpha A_{s,dir}H_{s,dir}/cos(z) +\alpha A_{s,ref}H_{s,dif}+\alpha r_g A_{s,ref}H_{s,ttl}, \]
where \(Q_{s,dir}\),\(Q_{s,dif}\),\(Q_{s,ref}\) are the direct, diffuse, and reflected solar radiative fluxes, respectively; \(H_{s,dir}\),\(H_{s,dif}\),\(H_{s,ttl}\) are the direct, diffuse, and total solar radiative horizontal flux densities (\(W/m^2\)), respectively; \(A_{s,dir}\),\(A_{s,ref}\),\(A_{s,ttl}\) are the direct, reflected, and total solar radiative heat transfer surface areas (\(m^2\)), respectively; \(\alpha\) is wing solar absorptivity; \(r_g\) is substrate solar reflectivity; and \(z\) is the zenith angle (rad). We assume \(A_{s,dir}=A_{s,ref}= A_{s,ttl}\).
Thermal radiative flux including both downward radiation and reflected solar radiation is estimated as follows:
\[Q_t= 0.5A_t \epsilon \sigma (T_{bK}^4-T_{skyK}^4)+0.5A_t \epsilon \sigma (T_{bK}^4-T_{gK}^4),\]
where \(A_t\) is the thermal radiative heat transfer surface area (\(m^2\)), \(T_{bK}\) is the body temperature (\(K\)), \(T_{gK}\) is the ground surface temperature (\(K\)), \(T_{skyK}\) is the equivalent black body sky temperature ($K), \(\epsilon\) is butterfly thermal emissivity, and \(\sigma\) is the Stefan-Boltzmann constant (\(5.67\times10^{-8} W m^{-2} K^{-4}\)).
The convective heat flux is given by:
\[Q_c= h_T A_c (T_{bK}-T_{aK}), \]
where \(A_c\) is the convective heat transfer surface area and \(T_{aK}\) is the air temperature (\(K\)). We assume \(A_c=A_t= A_{s,ttl}\). The total convective heat transfer coefficient, \(h_T Wm^{-2}K^{-1}\), is calculated as the boundary layer conductance \(h_c Wm^{-2}K^{-1}\) and the fur layer conduction in series:
\[1/h_T= 1/h_c+\frac{(r_i+\delta)ln[(r_i+\delta)/r_i]}{k_e} ,\]
where \(\delta\) is the thoracic fur thickness (\(m\)) and \(k_e\) is the thermal conductivity of the fur (\(Wm^{-1}K^{-1}\)).
The boundary layer conductance, \(h_c\), can be estimated using the relationship between two non-dimensional numbers. The Nusselt number, \(nu =h_c D/ka\), is the ratio of convective to conductive heat transfer, where \(k_a\) is the thermal conductivity of air (\(Wm^{-1}K^{-1}\)). We used the maximum width of the mesothorax as the characteristic dimension of the butterfly, \(D (m)\). The Reynolds number, \(Re =uD/v\), is the ratio of inertial forces to viscous forces, where \(u\) is wind speed (\(m/s\)) and \(v\) is kinematic viscosity (\(m^2/s\)). We used the Nu-Re relation for a cylinder, \(nu = 0.6Re^{0.5}\), which is a reasonable approximation for Colias. We estimate the energy flux balance to estimate \(T_b\).
We provide the operative temperature model for butterflies as follows:
t.seq <- lapply(20:40,
FUN = Tb_butterfly,
T_g = 25,
T_sh = 20,
u = 0.4,
S_sdir = 300,
S_sdif = 100,
z = 30,
D = 0.36,
delta = 1.46,
alpha = 0.6,
r_g = 0.3)
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)")
abline(a = 0,
b = 1,
col = "gray")
t.seq <- lapply(20:40,
FUN = Tb_butterfly,
T_g = 25,
T_sh = 20,
u = 0.4,
S_sdir = 100,
S_sdif = 100,
z = 30,
D = 0.36,
delta = 1.46,
alpha = 0.6,
r_g = 0.3)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tb_butterfly,
T_g = 25,
T_sh = 20,
u = 0.4,
S_sdir = 500,
S_sdif = 100,
z = 30,
D = 0.36,
delta = 1.46,
alpha = 0.6,
r_g = 0.3)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
legend(x = "bottomright",
title= expression("direct radiation," ~ S[sdif]),
legend = c(100,300,500),
lty = c("dashed", "solid", "dotted"))
Grasshoppers
We provide a biophysical model (Tb_grasshopper()) to predict grasshopper body temperatures (Buckley, Nufio, and Kingsolver 2014), partially adapted from (Samietz, Salser, and Dingle 2005). We use an energy budget to describe the flow of energy between the grasshopper and the environment: \(Q_s=Q_t+Q_c+Q_{cond}\). Here \(Q_s\) is the total input of heat due to solar radiation (\(W\)). \(Q_t\) describes the flux of thermal radiative heat due to both incoming thermal radiation (ground and sky) and that emitted by the grasshopper (\(W\)). \(Q_c\) is the flux of heat between the grasshopper and the surrounding fluid (air) via convection (\(W\)). \(Q_{cond}\) is the flux of heat between the grasshopper’s body and the solid surfaces with which the grasshopper’s body is in contact via conduction (\(W\)). We omit evaporative heat loss as it should be negligible for the grasshopper (Anderson, Tracy, and Abramsky 1979).
The solar radiative heat flux is estimated as the sum of direct (\(Q_{s,dir}\)), diffuse (\(Q_{s,dif}\)), and reflected (\(Q_{s,ref}\)) components (Kingsolver 1983):
\[ Q_s= Q_{s,dir}+Q_{s,dif}+Q_{s,ref},\\ Q_s= \alpha A_{s,dir}H_{s,dir}/cos(z) +\alpha A_{s,ref}H_{s,dif}+\alpha r_g A_{s,ref}H_{s,ttl}. \]
Each component is calculated as the product of the solar absorptivity of the grasshopper [we assume \(\alpha = 0.7\), (Anderson, Tracy, and Abramsky 1979)], the horizontal flux of solar radiation (\(H_{s,dir}\), \(H_{s,dif}\), and \(H_{s,ttl}\) for the direct, diffuse, and total fluxes, respectively), and the silhouette area of the grasshopper exposed to solar radiation (\(A_{s,dir}\), \(A_{s,dif}\), and \(A_{s,ttl}\) for the direct, diffuse, and total surface areas, respectively). The direct radiation is adjusted for the zenith angle (\(z\), degrees), which is the angle of the sun away from the vertical.
We calculate the surface area by approximating the body of a female grasshopper as a rotational ellipsoid (Samietz, Salser, and Dingle 2005). The major axis is equal to the grasshopper’s length. We calculate the semi-minor axis (half of the grasshopper’s width) as \(a = (0.365+0.241\times1000 L)/1000 m\) using a regression from Lactin and Johnson (1998). If \(e=\sqrt{1-a^2/c^2}\), the surface area can be calculated as follows:
\[A = 2 \pi a^2+\frac{2 \pi a c}{e \arcsin{e}}.\] The ratio of silhouette area to surface area of a grasshopper is a linear function of zenith angle: \(A_s/A = 0.19-0.00173z\). Thus, \(A_{s,dir}=A_{s,ref}=(0.19-0.00173z)A\). We partitioned the observed total radiation (\(H_{s,ttl}\)) into diffuse (\(H_{s,dif}\)) and direct (\(H_{s,dir}\)) components using the polynomial function of a clearness index, \(k_t\), developed by Erbs et al. (1982).
We estimate thermal radiative flux as the sum of radiation from the sky and ground. We assume that one half of the grasshopper’s body is subject to atmospheric radiation and the other half is subject to thermal radiation from the ground surface. Thermal radiation is calculated using the Stefan-Boltzmann law, which states that radiative flux is proportional to the fourth power of the absolute temperature of a body. Here \(T_{bK}\) is the body temperature (K), \(T_{gK}\) is the ground surface temperature (K), and \(T_{skyK}\) is the black body sky temperature (K) that corresponds to the air temperature [\(T_{skyK}= 0.0552T_{aK}^{1.5}\), (Swinbank 1963)]. The Stefan-Boltzmann constant (\(\sigma\)) characterizes the proportionality of this relationship. The thermal emissivity (\(\epsilon\)) accounts for incomplete absorption or emission of thermal radiation, but in this case, we assume that both the grasshopper and ground are perfect black bodies (\(\epsilon= 1\)). We account for the thermal radiative heat transfer surface area (\(A_t=A\)). The relationship is thus:
\[Q_t= 0.5A_t \epsilon \sigma (T_{bK}^4-T_{skyK}^4)+0.5A_t \epsilon \sigma (T_{bK}^4-T_{gK}^4).\]
Convective heat flux is estimated as the product of the convective heat transfer coefficient in turbulent air \(h_{cs} (Wm^{-2}C^{-1})\), the grasshopper’s surface area exposed to convective heat flux (\(A_c=A\)), and the temperature difference between the grasshopper’s body temperature \(T_{bK}\) and air temperature \(T_{aK}\):
\[Q_c= h_{cs} A_c (T_{bK}-T_{aK}).\]
We calculate \(h_{cs}\) from the convective heat transfer coefficient as \(h_{cs}=h_c(-0.007z/L+1.71)\) where \(z\)= 0.001m is the height above the ground. We use an empirically-derived relationship for grasshoppers to estimate the heat transfer coefficient, \(h_c (Wm^{-2}K^{-1})\) (Lactin and Johnson 1998): \(h_c = Nu K_f/L\), where the thermal conductivity of fluid, \(K_f(Wm^{-1}K^{-1})= 0.024+0.00007 T_{aK}\).
We use an empirical relationship from Anderson et al. (1998) to estimate the Nusselt number, \(Nu\), as \(nu = 0.41 Re^{0.5}\) where \(Re\) is the Reynolds number, \(Re =u L/v\), where \(u\) is wind speed (\(m/s\)) and \(v\) is the kinematic viscosity of air (\(m^2/s\)) (\(v = 15.68 \times 10^{-6}\) at \(300K\)). The rate of conduction is a function of the body area in contact with the substrate and the temperature differential between the body and the surface: \[Q_{cond}= h_{cut} A_{cond} (T_{bK}-T_{gK})/T, \] where \(h_{cut}\) is the thermal conductivity of the grasshopper cuticle [approximated as \(0.15 W m^{-1} K^{-1}\); value for hornets; (Galushko et al. 2005)]; \(A_{cond}\) is the surface area of the grasshopper in contact with the substrate; and \(T\) is the cuticle thickness [approximated as \(6 \times 10^{-5}m\); (Galushko et al. 2005)]. We only model conductance through the cuticle as we assume that the interior of the grasshopper is well mixed.
We provide the operative temperature model for grasshoppers as follows:
t.seq <- lapply(20:40,
FUN = Tb_grasshopper,
T_g = 25,
u = 0.4,
S = 800,
K_t = 0.7,
psi = 30,
l = 0.05,
Acondfact = 0.00,
z = 0.001,
abs = 0.7,
r_g = 0.3)
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)",
xlim = c(20,40),
ylim = c(20,40))
abline(a = 0,
b = 1,
col = "gray")
t.seq <- lapply(20:40,
FUN = Tb_grasshopper,
T_g = 25,
u = 0.4,
S = 400,
K_t = 0.7,
psi = 30,
l = 0.05,
Acondfact = 0.00,
z = 0.001,
abs = 0.7,
r_g = 0.3)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tb_grasshopper,
T_g = 25,
u = 0.4,
S = 600,
K_t = 0.7,
psi = 30,
l = 0.05,
Acondfact = 0.00,
z = 0.001,
abs = 0.7,
r_g = 0.3)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
legend(x = "bottomright",
title = expression("direct radiation, S"),
legend = c(200, 400, 600),
lty = c("dashed", "dotted", "solid"))
Mussels
We adapted a biophysical model for intertidal mussels that was developed by Helmuth (1998) to predict mussel body temperature \(T_b (C)\) based on environmental conditions, body size (height and length), and whether the mussel is solitary or aggregated. The amount of heat stored within a mussel body \(Q_{stored} (W)\) can be described as
\[Q_{stored} = Q_{sol} \pm Q_{rad,sky} \pm Q_{rad, ground} \pm Q_{cond} \pm Q_{conv} - Q_{evap} \tag{1}\]
where \(Q_{sol} (W)\) is the short-wave solar radiation to the mussel, \(Q_{rad,sky} (W)\) and \(Q_{rad,ground} (W)\) are the infrared radiation between the mussel and the sky and between the mussel and the surroundings, \(Q_{conv} (W)\) is the convective heat exchange, \(Q_{cond} (W)\) is the conduction to and from the ground, and \(Q_{evap} (W)\) is the evaporative cooling from water loss (B. S. T. Helmuth 1998). Since \(Q_{stored} (W)\) is a product of mass (\(m\)), specific heat (\(c\)), and body temperature (\(T_{bK}\), K), the rate of change of stored heat can be expressed as \(Q_{stored} = d(mcT_{bK})/dt\).
We calculate \(Q_{sol}\) from \(Q_{sol} = \frac{\alpha}{\sin\theta} A_{sol}S\) where \(\alpha\) is the solar absorptivity, \(\theta\) is the solar elevation angle (radian), \(A_{sol}\) is the mussel area in direction of the sun (\(m^2\)), \(S\) is the direct solar flux density on the area (B. S. T. Helmuth 1998). \(Q_{rad,sky}\) and \(Q_{rad,ground}\) can be calculated as follows (B. S. T. Helmuth 1998).
\[Q_{rad,sky} = \sigma A_{rad,sky}(\varepsilon_{org}T_{bK}^4-\varepsilon_{sky}T_{aK}^4) \tag{2}\] \[Q_{rad,ground} = \sigma A_{rad,ground}(\varepsilon_{org}T_{bK}^4-\varepsilon_{surr}T_{gK}^4) \tag{3}\]
Here, \(\sigma\) is the Stefan-Boltzmann constant (\(5.66 \times 10^{-8} W m^{-2} K^{-4}\)), \(T_{aK}\) and \(T_{gK}\) are air and ground temperatures (\(K\)), respectively, and \(A_{rad,sky}\) and \(A_{rad,ground}\) are the surface areas subject to long-wave radiation from sky and ground (\(m^2\)) respectively, both of which are simplified as half of the total surface area due to the nature of mussels. \(\varepsilon_{org}\) and \(\varepsilon_{surr}\) are the longwave thermal emissivities of the organism and its surroundings, which are roughly 0.97. \(\varepsilon_{sky}\) represents a functional infrared emissivity of the sky, which varies by air temperature and cloud cover.
In Environmental Biophysics (Campbell and Norman 2012), \(\varepsilon_{sky}\) is given by \(\varepsilon_{sky}=(1 - 0.84 c)(9.2 \times 10^{-6} T_{aK}^2) + 0.84c\) where \(c\) represents the fraction of the sky covered by cloud. Using Taylor series approximation, equations (2) and (3) become
\[Q_{rad,sky} = 4.0\varepsilon_{org}\varepsilon_{sky}^{3/4}\sigma A_{rad,sky}T_{aK}^3(T_{bK}-\varepsilon_{sky}^{1/4}T_{aK})\] \[Q_{rad,ground} = 4.0\varepsilon_{org}\sigma A_{rad,ground}T_g^3(T_{bK}-T_{gK})\]
Heat conduction \(Q_{cond}\) can be calculated as \(Q_{cond} = \frac{k_b}{0.5H}A_{cond}(T_{bK}-T_{gK})\) where \(k_b\) is the thermal conductivity of heat in the mussel body. Because most of the body of a mussel is composed of water, we approximated conductivity as equal to that of water, \(0.60 W m^{-2} K^{-1}\). \(H\) represents the height of the mussel and \(A_{cond}\) is the mussel area that is touching the ground.
Convective heat flux (\(Q_{conv}\)) can be calculated from \(Q_{conv}=h_cA_{conv}(T_{bK}-T_{aK})\) where \(h_c\) is the coefficient for forced convection and \(A_{conv}\) is the total surface area exposed to convective heat loss (B. S. T. Helmuth 1998). We calculate the value of \(h_c\) using three empirical formulas relating Nusselt numbers and Reynolds numbers: \(Nu = \frac{h_c L}{K_a}\), \(Re = \frac{uL}{v}\), and \(Nu = aRe^b\), where \(L\) is the length of the mussel, \(K_a\) is the conductivity of air \(K_a = 0.00501 + 7.2 \times10^{–5}T_a (W·m^{–1}·K^{–1})\), \(u\) is the wind speed (m/s), and \(v\) is the kinematic viscosity of air \(v = -1.25 \times 10^{-5} + 9.2 \times 10^{-8} T_{aK}\) (Denny and Harley 2006). \(A_{conv}\) can be assumed equal to the total surface area. The authors empirically acquired \(a\) and \(b\) in the formula to obtain \(Nu\), which varies depending on the position and the state of the mussels. For solitary mussels, \(a = 0.38\) and \(b = 0.51\) if the anterior or posterior end are facing upwind, and \(a = 0.63\) and \(b = 0.47\) if the valve is facing upwind. The values were measured to be \(a = 0.67\), \(b = 0.42\) for aggregated mussels (B. S. T. Helmuth 1998).
Lastly, \(Q_{evap} = \lambda \dot{m}\) where \(\lambda\) is the latent heat of vaporization of water (\(J/kg\)) and \(\dot{m}\) is the rate of water loss (B. S. T. Helmuth 1998). In this model, we assume the rate of water loss to be 5% of initial body mass per hour only when the mussels are gaping.
For simplicity, Tb_mussel function examines the steady-state (equilibrium) body temperature of an organism exposed over long periods to constant environmental conditions. See Helmuth (1998) for the transient model. Inserting the values in equation (1) and solving for \(T_{bK}\) gives us the following equation.
\[T_{bK,steady} = (k_1 A_{sol} S + k_2 A_{rad,sky} k_3 T_{aK}^4 + k_4 A_{rad,ground} T_{gK}^4 + k_5 A_{cond} T_{gK} + h_c A_{conv} T_{aK} - \lambda\dot{m}) / \] \[(k_2 A_{rad,sky} T_{aK}^3 + k_4 A_{rad, ground} T_{gK}^3 + k_5 A_{cond} + h_c A_{conv} + \dot{m}c)\]
Here, \(k_1\) to \(k_5\) are physical constants and measured coefficients. We provide the operative temperature model as follows.
t.seq <- lapply(20:40,
FUN = Tb_mussel,
l = 0.1,
h = 0.05,
T_g = 30,
S = 500,
k_d = 0.2,
u = 0.2,
psi = 30,
evap = TRUE,
cl = 0.5,
group = "solitary")
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)",
ylim = c(25, 45))
t.seq <- lapply(20:40,
FUN = Tb_mussel,
l = 0.1,
h = 0.05,
T_g = 30,
S = 500,
k_d = 0.2,
u = 1,
psi = 30,
evap = TRUE,
cl = 0.5,
group = "solitary")
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tb_mussel,
l = 0.1,
h = 0.05,
T_g = 30,
S = 500,
k_d = 0.2,
u = 2,
psi = 30,
evap = TRUE,
cl = 0.5,
group = "solitary")
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
abline(a = 0,
b = 1,
col = "gray")
legend(x = "bottomright",
title = expression("Wind speed, m/s"),
legend = c(0.2,1,2),
lty = c("solid", "dashed", "dotted"))
We also illustrate how the model can be applied to investigate the effect of mussel size and evaporation as in Helmuth (1998). The red line indicates the effect of evaporative heat loss.
t.seq <- lapply(seq(0.02, 0.14, 0.01),
FUN = Tb_mussel,
T_a = 25,
h = 0.05,
T_g = 30,
S = 500,
k_d = 0,
u = 0.5,
psi = 60,
evap = FALSE,
cl = 0,
group = "solitary")
plot(x = seq(0.02, 0.14, 0.01),
y = t.seq,
type = "l",
xlab = "mussel length (m)",
ylab = "body temperature (°C)",
ylim = c(25,32))
t.seq <- lapply(seq(0.02, 0.14, 0.01),
FUN = Tb_mussel,
T_a = 25,
h = 0.05,
T_g = 30,
S = 500,
k_d = 0,
u = 1,
psi = 60,
evap = FALSE,
cl = 0,
group = "solitary")
points(x = seq(0.02, 0.14, 0.01),
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(seq(0.02, 0.14, 0.01),
FUN = Tb_mussel,
T_a = 25,
h = 0.05,
T_g = 30,
S = 500,
k_d = 0,
u = 3,
psi = 60,
evap = FALSE,
cl = 0,
group = "solitary")
points(x = seq(0.02, 0.14, 0.01),
y = t.seq,
type = "l",
lty = "dotted")
t.seq <- lapply(seq(0.02, 0.14, 0.01),
FUN = Tb_mussel,
T_a = 25,
h = 0.05,
T_g = 30,
S = 500,
k_d = 0,
u = 1,
psi = 60,
evap = TRUE,
cl = 0,
group = "solitary")
points(x = seq(0.02, 0.14, 0.01),
y = t.seq,
type = "l",
col = "red")
legend(x = "bottomright",
title = expression("Wind speed, m/s"),
legend = c(0.5,1,3),
lty = c("solid", "dashed", "dotted"))
TrenchR also provides a version of the model modified for mussel beds (B. Helmuth 1999). The red line indicates the effect of evaporative heat loss.
t.seq <- lapply(20:40,
FUN = Tbed_mussel,
l = 0.1,
S = 500,
k_d = 0.2,
u = 0.2,
evap = FALSE)
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)",
ylim = c(20, 50))
t.seq <- lapply(20:40,
FUN = Tbed_mussel,
l = 0.1,
S = 500,
k_d = 0.2,
u = 1,
evap = FALSE)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tbed_mussel,
l = 0.1,
S = 500,
k_d = 0.2,
u = 2,
evap = FALSE)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
t.seq <- lapply(20:40,
FUN = Tbed_mussel,
l = 0.1,
S = 500,
k_d = 0.2,
u = 0.2,
evap = TRUE)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "solid",
col = "red")
abline(a = 0,
b = 1,
col = "gray")
legend(x = "bottomright",
title = expression("Wind speed, m/s"),
legend = c(0.2,1,2),
lty = c("solid", "dashed", "dotted"))
Additionally, land surface models have been modified to predict mussel bed temperatures (Mislan and Wethey 2015). We do not include these more complex models in TrenchR because they are available elsewhere (Mislan and Wethey 2014: https://zenodo.org/record/13380#.XpidHJNKjm0)
Limpets
We adapted a biophysical model Tb_Limpet() from Denny and Harley (2006) to predict the body temperature of limpets based on body size (length and height), environmental conditions, and the position of the limpet. We note that our implementation is simplified from the original. In this model, the shape of the limpet is approximated to a cone of radius \(r\) and height \(H\).
The rate at which energy is stored within the organism \(W_{stored} (W)\) is influenced by several heat exchange processes. The processes include the rate of heat absorbed from short-wave radiation \(W_{sw} (W)\), net heat transferred by long-wave radiation \(W_{lw} (W)\), and the rate of heat transferred by conduction and convection \(W_{cd} (W)\) and \(W_{cv} (W)\). We neglected heat gain by metabolism due to the limpets’ low metabolic rate. We likewise omitted evaporative heat loss assuming that the limpet has its shell firmly clamped against the substratum, which limits evaporation.
The rate at which solar radiation provides heat is defined as
\[W_{sw} = A_p\alpha_{sw}I_{sw}\]
where \(A_p\) is the area of the limpet’s shell (\(m^2\)) projected in the direction of the sunlight, \(\alpha_{sw}\) is the short-wave absorptivity of the shell, and \(I_{sw}\) is the solar irradiance (\(W m^{-2}\)). Under the assumption that limpet has a cone shape, where \(tan(\psi) > r / h\),
\[A_p = \pi r^2 cos(\psi).\]
Otherwise,
\[A_p = \pi r^2 cos(\psi) + hrsin(\psi) - \frac{\pi r^2}{2}cos(\psi)\]
where \(\psi\) is the solar zenith angle (Pennell and Deignan 1989). \(\alpha_{sw}\) was empirically acquired to be 0.68 for limpets (Denny and Harley 2006).
Net rate of long-wave energy transfer (\(W_{lw}\)) is the difference between the rate at which long-wave radiation from the sky is absorbed by the shell (\(W_{lw,a}\)) and the total effective rate at which energy is radiated from limpet to sky (\(W_{lw,s}\)). It is expressed as
\[W_{lw} = V_s A_l \epsilon_{lw,s}\sigma(\epsilon_{lw,a}\sigma T_{aK}^4 - T_{bK}^4)\]
where \(V_s\) is the view factor, \(A_l\) is the lateral area of limpet shell (\(m^2\)), \(\epsilon_{lw,s}\) is long-wave emissivity of the shell (= 0.97), and \(\epsilon_{lw,a}\) is the emissivity of air. \(V_s\) equals to \(cos(\psi)\cdot r / \sqrt{r^2 + H^2}\) by simulating limpet as a cone (Campbell and Norman 2012).
This can be linearized by Taylor expansion around \(T_{aK}\) to yield
\[W_{lw} = V_s A_l\epsilon_{w,s}\sigma T_{aK}^4(\epsilon_{w,a}-1) + 4V_sA_l\epsilon_{w,s}\sigma T_{aK}^3(T_{aK} - T_{bK}).\]
Convective heat flux (\(W_{cv}\)) can be estimated as
\[W_{cv} = h_cA_{cv}(T_{aK}-T_{bK}),\]
which is derived from the Newton’s law of cooling. \(A_{cv}\) is the area of the shell in contact with the air (\(m^2\)), which is approximated to the lateral area of the limpet (\(A_l\)). The convective heat transfer coefficient, \(h_c\) (\(W m^{-2} K^{-1}\)), is a component of the Nusselt number, \(Nu = \frac{h_c R}{K_a}\) where \(R\) is the limpet’s diameter and \(K_a\) is the conductivity of air (\(W m^{-1} K^{-1}\)). Denny (1993) found that conductivity varies with air temperature as follows: \(K_a = 0.00501 + 7.2 \times 10^{-5}T_a\). Additionally,
\[Re = \frac{uR}{v}\]
and
\[Nu = aRe^b\]
where \(u\) is the wind speed (\(m/s\)), \(v\) is the kinematic viscosity of air (\(m^2/s\)) that equals \(-1.25 \times 10^{-5} + 9.2 \times 10^{-8} T_{aK}\), and \(a\) and \(b\) are coefficients that need to be empirically obtained. When the anterior side of the limpet is facing upwind, \(a = 1.955\) and \(b= 0.371\). When the posterior side is facing upwind, \(a = 1.881\) and \(b= 0.376\). Lastly when the broadside is facing upwind, \(a = 1.304\) and \(b= 0.404\) (Denny and Harley 2006).
Conductive heat exchange is determined by the following equation:
\[W_{cd} = A_{cd}K_r\left(\frac{dT_{rK}}{dz}\right),\]
where \(A_{cd}\) is the area of contact between the organism’s foot and the rock (\(m^2\)), and \(K_r\) is the conductivity of the rock (\(3.06 W m^{-1} K^{-1}\)). \(T_{rK} (K)\) is the temperature of the rock as a function of distance \(z (m)\) from the surface.
Finally, \(T_{bK}\) can be calculated from the above equations and
\[W_{sw} \pm W_{lw} \pm W_{cd} \pm W_{cv} = 0.\]
We provide the operative temperature model as follows:
t.seq <- lapply(20:40,
FUN = Tb_limpet,
T_r = 30,
l = 0.0176,
h = 0.0122,
S = 1000,
u = 1,
psi = 30,
c = 1,
position = "anterior")
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)",
ylim = c(25, 50))
t.seq <- lapply(20:40,
FUN = Tb_limpet,
T_r = 30,
l = 0.0176,
h = 0.0122,
S = 1300,
u = 1,
psi = 30,
c = 1,
position = "anterior")
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tb_limpet,
T_r = 30,
l = 0.0176,
h = 0.0122,
S = 1600,
u = 1,
psi = 30,
c = 1,
position = "anterior")
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
abline(a = 0,
b = 1,
col = "gray")
legend(x = "bottomright",
title = expression("Solar irradiance, W/m^2"),
legend = c(1000, 1300, 1600),
lty = c("solid", "dashed", "dotted"))
TrenchR also provides two other versions of the limpet model modified by the Helmuth Lab. The first modifies radiation and convection in the function for limpets: Tb_limpetBH. The second modifies the function to predict snail body temperatures Tb_snail. We illustrate the functions below:
t.seq <- lapply(20:40,
FUN = Tb_limpet,
T_r = 30,
l = 0.0176,
h = 0.0122,
S = 800,
u = 1,
psi = 30,
c = 1,
position = "anterior")
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)",
ylab = "body temperature (°C)",
ylim = c(25, 50))
t.seq <- lapply(20:40,
FUN = Tb_limpetBH,
T_r = 30,
l = 0.0176,
h = 0.0122,
S = 800,
u = 1,
s_aspect = 90,
s_slope = 60,
c = 1)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dashed")
t.seq <- lapply(20:40,
FUN = Tb_snail,
l = 0.012,
S = 800,
u = 1,
CC = 0.5,
WL = 0,
WSH = 10)
points(x = 20:40,
y = t.seq,
type = "l",
lty = "dotted")
abline(a = 0,
b = 1,
col = "gray")
legend(x = "bottomright",
title = expression("function"),
legend = c("Tb_limpet", "Tb_limpetBH", "Tb_Snail"),
lty = c("solid", "dashed", "dotted"))
Nocturnal Salamanders (humid operative temperatures)
We present a function Tb_salamander_humid adapted from Riddell et al. (2018) and Campbell and Norman (2012) that calculates the humid operative temperature (\(T_{eh}\)) of a nocturnal salamander. The function estimates \(T_{eh}\) by accounting for evaporative cooling from the high water loss rates typically associated with amphibians. Rates of evaporative cooling are primarily a function of the humidity in the air and the resistance of the skin to water loss. The vapor pressure deficit (\(D\)) describes the drying power of the ambient air and thus the potential for evaporative cooling.
Similarly, the skin resistance to water loss (\(r_{i}\)) describes the capacity of the integument to resist the flux of water given the gradient of humidity in ambient air. High rates of water loss, for instance, can occur if humidity or skin resistance to water loss is low, resulting in greater amounts of evaporative cooling. The model also incorporates additional variables, such as nighttime absorbed radiation, radiative conductance, elevation, and boundary layer resistance. Direct and indirect solar radiation are not included because the code was developed for a nocturnal species. In addition to the cooling effect from evaporation, operative temperatures (\(T_{e}\)) at night are also generally lower than air temperatures due to emitted radiation being larger than absorbed (Campbell and Norman 2012).
The \(T_{eh}\) function is defined as \[T_{eh}=T_{a}+\frac{\gamma^{*}}{s+\gamma^{*}} \left( \frac{R_{abs} - \epsilon_{s}\sigma T^{4}_{a}}{c_{p}g_{Hr}} -\frac{D} {\gamma^{*}p_{a}} \right)\]where \(T_{a}\) is the air temperature in Celsius, \(\gamma^{*}\) is the apparent psychrometer constant (\(C^{-1}\)), \(s\) is the slope of the saturation mole fraction function (\(C^{-1}\)), \(R_{abs}\) is the absorbed radiation (\(W\:m^{-2}\)), \(\epsilon _{s}\) is the emissivity of the surface, \(\sigma\) is the Stefan-Boltzmann constant (\(W\:m^{−2}\:C^{−4}\)), \(c_{p}\) is specific heat of air (29.3 \(J\:mol^{-1}\:C^{-1}\)), \(g_{Hr}\) is the sum of boundary layer conductance for heat and radiative conductance (\(mol\:m^{-2}\:s^{-1}\)), \(D\) is the vapor pressure deficit (vpd, kPa), and \(p_{a}\) is the atmospheric pressure (kPa).
To calculate \(s\), we use \(s=\Delta/p_{a}\) where \(\Delta\) is the slope of saturation vapor pressure.
We calculate \(\Delta\) using
\[\Delta = \frac{bce_{s}(T)}{(c+T)^2}\]
where \(b = 17.502\), \(c = 240.97^{\circ}C\), \(e_s(T)\) is the saturation vapor pressure of the air (\(kPa\)), and T is temperature in Celsius. For \(p_a\) we use
\[p_{a}= 101.3exp\Big(\frac{-A}{8200}\Big)\]
where \(A\) is the altitude in meters above sea level and \(p_{a}\) is in \(kPa\). For \(e_s{s}\) we use \[e_{s}(T)=a\:exp\ Big(\frac{bT}{T+c}\Big)\] where \(a\) is 0.611 kPa and \(T\) is the temperature in Celsius.
For \(\gamma^{*}\), we use \[\gamma^{*}=\frac{\gamma(1/g_{vs} + 1/g_{vc}+1/g_{va})} {1/g_{Hc}+1/g_{Hr}}\] where \(\gamma\) is the thermodynamic psychrometer constant, which is \(6.66 \times 10^{-4}C^{-1}\) (Campbell and Norman 2012), \(g_{vs}\) is the conductance of the skin to vapor (\(mol\:m^{−2}\:s^{−1}\)), \(g_{vc}\) is the conductance of the coat to vapor (\(mol\:m^{−2}\:s^{−1}\)), \(g_{va}\) is the boundary layer conductance for vapor (\(mol\:m^{−2}\:s^{−1}\)), \(g_{Hc}\) is the conductance of the coat for heat (\(mol\:m^{−2}\:s^{−1}\)), and \(g_{Hr}\) is the sum of boundary layer conductance for heat and radiative conductance (\(mol\:m^{−2}\:s^{−1}\)).
Both \(1/g_{Hc}\) and \(1/g_{vc}\) are assumed to be zero since salamanders do not have insulation, such as fur, feathers, or clothing. Thus, \(\gamma^{*}\) can also be calculated as \(\gamma^{*}=\gamma H_{r}/g_{v}\) where \(g_{v}\) is the total vapor conductance (\(mol\:m^{−2}\:s^{−1}\)). We calculate \(g_{v}\) using
\[g_{v}=\frac{g_{vs} g_{va}}{g_{vs}+g_{va}}.\]
We convert skin and boundary layer resistance for vapor to conductance by multiplying the values by 100 to convert to \(s\:m^{-2}\) and then dividing by 41.1 \(mol\:m^{-3}\). The boundary layer conductance for heat is calculated as \(g_{Ha}= 1.4*0.135*sqrt(u/d)\), where \(u\) is the wind speed \((m/s)\) and \(d\) is the characteristic dimension \((m)\). We use a wind speed 0.1 m/s in our calculations of \(g_{Ha}\) to mimic still air near the ground (Campbell and Norman 2012), and for \(D\) the cross-sectional diameter of a cylinder approximating the size of a salamander. Boundary layer resistance can be estimated as follows.
par(mar = c(5, 5, 3, 2))
blr.seq= lapply(293:313, FUN=boundary_layer_resistance, e_s= 2.5, e_a = 2.3, elev = 500, D = 0.007, u = 2)
plot(x = 293:313,
y = blr.seq,
type = "l",
xlab = "air temperature (K)",
ylab = expression("boundary layer resistance" ~ (s ~ cm^{-1})))
TrenchR includes several support functions, often adapted from Campbell and Norman (2012), used for estimating humid operative temperatures. Functions estimate thermal radiation absorbed (thermal_radiation_absorbed) and emitted (Qemitted_thermal_radiation) as well as actual (actual_vapor_pressure) and saturation (saturation_vapor_pressure) vapor pressure. The functions also include a statistical approximation from Campbell and Norman (2012) for estimating soil temperatures at specified depths.
plot(x = 1:24,
y = Tsoil(T_g_max = 30,
T_g_min = 15,
hour = 1:24,
depth = 0),
type = "l",
xlab = "hour",
ylab = "soil temperature (°C)")
points(x = 1:24,
y = Tsoil(T_g_max = 30,
T_g_min = 15,
hour = 1:24,
depth = 5),
type = "l",
lty = "dashed")
points(x = 1:24,
y = Tsoil(T_g_max = 30,
T_g_min = 15,
hour = 1:24,
depth = 10),
type = "l",
lty = "dotted")
points(x = 1:24,
y = Tsoil(T_g_max = 30,
T_g_min = 15,
hour = 1:24,
depth = 20),
type = "l",
lty = "dotdash")
legend(x = "topleft",
title = "depth (cm)",
legend = c(0,5,10,20),
lty = c("solid", "dashed", "dotted", "dotdash"))
We illustrate the humid operative temperature model by comparing an example parameterization (solid lines) to a 1:1 line (dotted line) indicating when humid operative temperatures equal air temperatures. Evaporative cooling lowers humid operative temperatures below air temperatures for the nocturnal salamander. Higher skin resistance to evaporation reduced cooling (red) whereas higher vapor pressure deficit augments cooling (blue).
rad.seq <- sapply(20:40,
FUN = Qthermal_radiation_absorbed,
T_g = 30,
epsilon_ground = 0.97,
a_longwave = 0.965)
t.seq <- mapply(FUN = Tb_salamander_humid,
T_a = 20:40,
r_i = 4,
r_b = 1,
D = 0.01,
elev = 500,
e_a = 1.5,
e_s = 2.5,
Qabs = rad.seq,
epsilon = 0.96)
plot(x = 20:40,
y = t.seq,
type = "l",
xlab = "ambient temperature (°C)", ylab = "humid operative temperature (°C)")
abline(a = 0,
b = 1,
lty = "dotted")
# check higher skin resistance-> Te closer to air due to reduced cooling
t.seq <- mapply(FUN = Tb_salamander_humid,
T_a = 20:40,
r_i = 20,
r_b = 1,
D = 0.01,
elev = 500,
e_a = 1.5,
e_s = 2.5,
Qabs = rad.seq,
epsilon = 0.96)
points(x = 20:40,
y = t.seq,
type = "l",
col = "red")
# check higher vpd-> lower Te due to greater cooling
t.seq <- mapply(FUN = Tb_salamander_humid,
T_a = 20:40,
r_i = 4,
r_b = 1,
D = 0.01,
elev = 500,
e_a = 1.0,
e_s = 2.5,
Qabs = rad.seq,
epsilon = 0.96)
points(x = 20:40,
y = t.seq,
type = "l",
col = "blue")
legend(x = "topleft",
title = NULL,
legend = c("base scenario", "higher skin resistance", "higher vpd"),
lty = c("solid"),
col = c("black", "red", "blue"))
References
Anderson, R. V., C. R. Tracy, and Z. Abramsky. 1979. “Habitat Selection in Two Species of Short-Horned Grasshoppers.” Oecologia 38 (3): 359–74.
Bakken, G. S. 1992. “Measurement and Application of Operative and Standard Operative Temperatures in Ecology.” American Zoologist 32 (2): 194–216.
Barlett, P. N., and D. M. Gates. 1967. “The Energy Budget of a Lizard on a Tree Trunk.” Ecology 48 (2): 315–22.
Buckley, L. B. 2007. “Linking Traits to Energetics and Population Dynamics to Predict Lizard Ranges in Changing Environments.” The American Naturalist 171 (1): E1–19.
Buckley, L. B., and J. G. Kingsolver. 2012. “The Demographic Impacts of Shifts in Climate Means and Extremes on Alpine Butterflies.” Functional Ecology 26 (4): 969–77.
Buckley, L. B., C. R. Nufio, and J. G. Kingsolver. 2014. “Phenotypic Clines, Energy Balances and Ecological Responses to Climate Change.” Journal of Animal Ecology 83 (1): 41–50.
Campbell, G. S., and J. Norman. 2012. An Introduction to Environmental Biophysics. Springer Science & Business Media.
Denny, M. W. 1993. Air and Water: The Biology and Physics of Life’s Media. Princeton University Press.
Denny, M. W., and C. D. G. Harley. 2006. “Hot Limpets: Predicting Body Temperature in a Conductance-Mediated Thermal System.” Journal of Experimental Biology 209 (13): 2409–19.
Erbs, D. G., S. A. Klein, and J. A. Duffie. 1982. “Estimation of the Diffuse Radiation Fraction for Hourly, Daily and Monthly-Average Global Radiation.” Solar Energy 28 (4): 293–302.
Fei, T., A.w K. Skidmore, V. Venus, T. Wang, M. Schlerf, B. Toxopeus, S. Van Overjijk, M. Bian, and Y. Liu. 2012. “A Body Temperature Model for Lizards as Estimated from the Thermal Environment.” Journal of Thermal Biology 37 (1): 56–64.
Forsythe, W. C., E. J. Rykiel Jr., R. S. Stahl, H. Wu, and R. M. Schoolfield. 1995. “A Model Comparison for Daylength as a Function of Latitude and Day of Year.” Ecological Modeling 80 (1): 87–95.
Galushko, D., N. Ermakov, M. Karpovski, A. Palevski, J. S. Ishay, and D. J. Bergman. 2005. “Electrical, Thermoelectric and Thermophysical Properties of Hornet Cuticle.” Semiconductor Science and Technology 20 (3): 286.
Gates, D. M. 2012. Biophysical Ecology. Courier Corporation.
Gillooly, J. F., J. H. Brown, G. B. West, V. M. Savage, and E. L. Charnov. 2001. “Effects of Size and Temperature on Metabolic Rate.” Science 293 (5538): 2248–51.
Helmuth, B. 1999. “Thermal Biology of Rocky Intertidal Mussels: Quantifying Body Temperatures Using Climatological Data.” Ecology 80 (1): 15–34.
Helmuth, B. S. T. 1998. “Intertidal Mussel Microclimates: Predicting the Body Temperature of a Sessile Invertebrate.” Ecological Monographs 68 (1): 51–74.
Kingsolver, J. G. 1983. “Thermoregulation and Flight in Colias Butterflies: Elevational Patterns and Mechanistic Limitations.” Ecology 64 (3): 534–45.
Lactin, D. J., and D. L. Johnson. 1998. “Environmental, Physical, and Behavioural Determinants of Body Temperature in Grasshopper Nymphs (Orthoptera: Acrididae).” The Canadian Entomologist 130 (5): 551–77.
Liu, B. Y. H., and R. C. Jordan. 1960. “The Interrelationship and Characteristic Distribution of Direct, Diffuse and Total Solar Radiation.” Solar Energy 4 (3): 1–19.
Mislan, K. A. S., and D. S. Wethey. 2015. “A Biophysical Basis for Patchy Mortality During Heat Waves.” Ecology 96 (4): 902–7.
Mitchell, John W. 1976. “Heat Transfer from Spheres and Other Animal Forms.” Biophysical Journal 16 (6): 561–69.
Muth, A. 1977. “Thermoregulatory Postures and Orientation to the Sun: A Mechanistic Evaluation for the Zebra-Tailed Lizard, Callisaurus Draconoides.” Copeia, 710–20.
Norris, K. S. 1967. “Color Adaptation in Desert Reptiles and Its Thermal Relationships.” In Lizard Ecology: A Symposium, 162–229. University of Missouri Press Columbia.
Pennell, S., and J. Deignan. 1989. “Computing the Projected Area of a Cone.” SIAM Review 31 (2): 299–302.
Porter, W. P., and F. C. James. 1979. “Behavioral Implications of Mechanistic Ecology II: The African Rainbow Lizard, Agama Agama.” Copeia, 594–619.
Porter, W. P., J. W. Mitchell, W. A. Beckman, and C. B. DeWitt. 1973. “Behavioral Implications of Mechanistic Ecology.” Oecologia 13 (1): 1–54.
Riddell, E. A., J. P. Odom, J. D. Damm, and M. W. Sears. 2018. “Plasticity Reveals Hidden Resistance to Extinction Under Climate Change in the Global Hotspot of Salamander Diversity.” Science Advances 4 (7): eaar5471.
Roughgarden, J., W. Porter, and D. Heckel. 1981. “Resource Partitioning of Space and Its Relationship to Body Temperature in Anolis Lizard Populations.” Oecologia 50 (2): 256–64.
Samietz, J., M. A. Salser, and H. Dingle. 2005. “Altitudinal Variation in Behavioural Thermoregulation: Local Adaptation Vs. Plasticity in California Grasshoppers.” Journal of Evolutionary Biology 18 (4): 1087–96.
Spotila, J. R., M. P. O’Connor, and G. S. Bakken. 1992. “Biophysics of Heat and Mass Transfer.” Environmental Physiology of the Amphibians. University of Chicago Press, Chicago, 59–80.
Swinbank, W. C. 1963. “Long-Wave Radiation from Clear Skies.” Quarterly Journal of the Royal Meteorological Society 89 (381): 339–48.
