---
title: "Using energy balances to estimate body temperatures"
author: "Lauren Buckley and the TrEnCh project, Department of Biology, University of Washington"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: bibliography.bib
vignette: >
  %\VignetteIndexEntry{Using energy balances to estimate body temperatures}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown:  
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}

# Load TrenchR package

 library(TrenchR)

```

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.

# Operative environmental temperatures

The body temperatures, $T_b$, of organisms can depart dramatically from air temperatures due to energy exchange with the environment. 
Heat is gained by absorbing solar and thermal radiation and is generated from metabolic reactions. 
Heat is lost through the organism's emission of radiation and the evaporation of water. 
The organism exchanges heat with the surrounding air or water via convection and with the ground and other solid surfaces in contact via conduction. 
The balance of these heat exchanges can be estimated and is often referred to as operative environmental temperature, $T_e$. 
The $T_e$ indicates the steady-state temperature of an organism with given physical properties in a particular environment, omitting heat gain via metabolism and heat loss via evaporation [@bakken1992measurement].

The following accounting of heat gains and losses at the surface of the organism can be used to estimate $T_b$ (or $T_e$ if $Q_{met}=0$ and $Q_{evap}=0$)  [@gates2012biophysical]:

$$Q_{net} = Q_{abs} - Q_{emit} - Q_{conv} - Q_{cond} - Q_{met} - Q_{evap},$$

where $Q_{net}$ is the net energy exchange with the environment (W), $Q_{abs}$ is the solar radiation absorbed (W), $Q_{emit}$ is the net thermal radiation emitted (W), $Q_{conv}$ is energy exchange due to convection (W), $Q_{cond}$ is energy exchange due to conduction (W), $Q_{met}$ is the energy emitted due to metabolism (W), and $Q_{evap}$ is the energy emitted due to evaporative water loss (W).

The energy balance is available in TrenchR as follows:

```{r}
Qnet_Gates(Qabs  = 500, 
           Qemit = 10, 
           Qconv = 100, 
           Qcond = 100,
           Qmet  = 10, 
           Qevap = 5)
```

We briefly review each component of the energy balance below. 
Details are available in biophysical ecology texts [@gates2012biophysical; @campbell2012introduction]. We note that we have standardized most functions to accept and return temperatures in degrees Celsius to promote accessibility. However, the heat budget calculations generally use Kelvin since it is the natural unit, particularly for radiation, and metrics such as convective coefficients standardly use kelvin. TrenchR includes tables of parameters such as the absorptivity of animal surfaces to solar and thermal radiation in the `inst/extdata` folder.

# 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:

```{r, fig.height = 4, fig.width = 6}

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.

# Thermal radiation, Qemit

The net rate of emission of thermal radiation from the surface of an animal, $Q_{emit} (W)$, is determined by the difference between the surface temperature of the animal $T_b (K)$ and the temperatures of the air $T_a (K)$ and of the ground $T_g (K)$. Exchange is thermal radiation between the animal and the air is estimated using sky temperature $T_{sky} (K)$, the effective radiant temperature of the sky. TrenchR includes a function for estimating $T_{sky} (K)$ using two commonly approaches: the Brunt or Swinbank formulas [@gates2012biophysical]. 
```{r}
T_sky(T_a=293, formula="Swinbank")
```


$T_a$ is additionally used to estimate sky temperature 
The following expressions can be used to estimate $Q_{emit} (W)$ for animals in enclosed and open environments, respectively. The functions uses the Brunt equation to estimate $T_{sky} if a value is not provided: $T_{sky}= 1.22\times (T_a-273.15)-20.4+273.15$ [@gates2012biophysical]. 
$$
enclosed: Q_{emit}= A_r \epsilon \sigma (T_b^4 - T_a^4)\\
open: Q_{emit}= \epsilon \sigma (A_s (T_b^4 - T_{sky}^4)+A_r (T_b^4 - T_g^4)),
$$ 

where $A_s$ and $A_r$ are the areas ($m^2$) exposed to the sky (or enclosure) and the ground, respectively; $\epsilon$ is the longwave infrared emissivity of skin [proportion, 0.95 to 1 for most animals [@gates2012biophysical]]; and $\sigma$ is the Stefan-Boltzmann constant ($5.673\times10^{-8} W m^{-2} K^{-4}$).

The function is available in R as follows:

```{r, fig.height = 4, fig.width = 6}

plot(x    = 293:313, 
     y    = Qemitted_thermal_radiation(epsilon  = 0.96, 
                                       A        = 1, 
                                       psa_dir  = 0.4, 
                                       psa_ref  = 0.5,  
                                       T_b      = 293:313,  
                                       T_g      = 293,  
                                       T_a      = 298,  
                                       enclosed = FALSE), 
     type = "l", 
     xlab = "body surface temperature, Tb (K)", 
     ylab = "emitted thermal radiation, Qemit (W)")
points(x    = 293:313, 
       y    = Qemitted_thermal_radiation(epsilon  = 0.96, 
                                         A        = 1, 
                                         psa_dir  = 0.2, 
                                         psa_ref  = 0.5, 
                                         T_b      = 293:313, 
                                         T_g      = 293, 
                                         T_a      = 298, 
                                         enclosed = TRUE), 
       type = "l", 
       lty  = "dashed")
points(x    = 293:313, 
       y    = Qemitted_thermal_radiation(epsilon  = 0.96,
                                         A        = 1,
                                         psa_dir  = 0.4, 
                                         psa_ref  = 0.5, 
                                         T_b      = 293:313, 
                                         T_g      = 283, 
                                         T_a      = 298, 
                                         enclosed = FALSE), 
       type = "l", 
       lty  = "dotted")
points(x    = 293:313, 
       y    = Qemitted_thermal_radiation(epsilon= 0.96, 
                                         A = 1, 
                                         psa_dir = 0.2, 
                                         psa_ref = 0.5, 
                                         T_b      = 293:313, 
                                         T_g      = 283, 
                                         T_a      = 298, 
                                         enclosed = FALSE), 
       type = "l", 
       lty  = "dotdash")

legend(x      = "topleft", 
       title  = "parameters", 
       legend = c("psa_dir= 0.4, T_g = 293", "psa_dir= 0.2, T_g = 293", "psa_dir= 0.4, T_g = 283", "psa_dir= 0.2, T_g = 283"), 
       lty    = c("solid", "dashed", "dotted", "dotdash"))

```

# 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:

```{r, fig.height = 4, fig.width = 6}

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 [@mitchell1976heat].

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:

```{r, fig.height = 4, fig.width = 6}

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"))

```

```{r, fig.height = 4, fig.width = 6}

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"))

```

```{r, fig.height = 4, fig.width = 6}

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"))


```

```{r, fig.height = 4, fig.width = 6}

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:

```{r, fig.height = 4, fig.width = 6}

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" ))

```

```{r}

Nusselt_from_Grashof(Gr = 5)

```

```{r}

free_or_forced_convection(Gr = 100,  
                          Re = 5)

```

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 [@mitchell1976heat]. 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.

```{r, fig.height = 4, fig.width = 6}

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"))

```

```{r, fig.height = 3, fig.width = 6}

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"))

```

```{r, fig.height = 3, fig.width = 5}

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 [@mitchell1976heat].

The functions are available in R:

```{r, fig.height = 4, fig.width = 6}

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"))

```

```{r, fig.height = 4, fig.width = 6}

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"))

```

# Metabolism, Qmet

Rates of heat production associated with metabolism are generally estimated using scaling relationships based on generalizing empirical data. Metabolic rates scale as a power law function of mass with an exponent less than 1 such that the per-mass metabolic rate declines with increasing mass. Metabolic rates also increase exponentially with temperature.  The general form of the relationship is $B=M^x \exp(-E_i/{kT})$, where $M$ is mass, $x$ is an exponent approximating 0.75, $E_i$ is the average activation energy for the rate-limiting enzyme-catalyzed biochemical reactions of metabolism, $k$ is the Boltzmann constant, and $T$ is body temperature [@gillooly2001effects]. We offer functions that both do not [`Qmetabolism_from_mass()`] and do [`Qmetabolism_from_mass_temp()`] include temperature dependence for a variety of taxa.

```{r, fig.height = 4, fig.width = 4}

plot(x    = 1:100, 
     y    = Qmetabolism_from_mass(m     = 1:100,
                                  taxon = "bird"), 
     type = "l", 
     xlab = "mass (g)", 
     ylab = "metabolism (W)", 
     log  = "xy", 
     ylim = c(0.1, 2))

points(x    = 1:100, 
       y    = Qmetabolism_from_mass(m     = 1:100,
                                    taxon = "mammal"), 
       type = "l", 
       lty  = "dashed")

points(x    = 1:100, 
       y    = Qmetabolism_from_mass(m     = 1:100,
                                    taxon = "reptile"), 
       type = "l", 
       lty  = "dotted")

legend(x      = "topleft", 
       title  = "taxa", 
       legend = c("bird", "mammal", "reptile"), 
       lty    = c("solid", "dashed", "dotted"))

```

```{r, fig.height = 4, fig.width = 4}

plot(x    = 1:100, 
     y    = Qmetabolism_from_mass_temp(m     = 1:100, 
                                       T_b   = 30,
                                       taxon = "bird"), 
     type = "l", 
     xlab = "mass (g)", 
     ylab = "metabolism (W)", 
     ylim = c(0, 0.1))

points(x    = 1:100, 
       y    = Qmetabolism_from_mass_temp(m     = 1:100, 
                                         T_b   = 30,
                                         taxon = "reptile"), 
       type = "l", 
       lty  = "dotted")

points(x    = 1:100, 
       y    = Qmetabolism_from_mass_temp(m     = 1:100, 
                                         T_b   = 30,
                                         taxon = "amphibian"), 
       type = "l", 
       lty  = "dotdash")

points(x    = 1:100, 
       y    = Qmetabolism_from_mass_temp(m     = 1:100,
                                         T_b   = 30,
                                         taxon = "invertebrate"), 
       type = "l", 
       lty  = "twodash")

legend(x      = "top", 
       title  = "taxa", 
       legend = c("bird,mammal", "reptile", "amphibian", "invertebrate"), 
       lty    = c("solid", "dotted", "dotdash", "twodash"))

plot(x    = 20:40, 
     y    = Qmetabolism_from_mass_temp(m     = 5, 
                                       T_b   = 20:40,
                                       taxon = "bird"), 
     type = "l", 
     xlab = "temperature (C)", 
     ylab = "metabolism (W)", 
     ylim = c(0,0.1))

points(x    = 20:40, 
       y    = Qmetabolism_from_mass_temp(m     = 5, 
                                         T_b   = 20:40,
                                         taxon = "reptile"), 
       type = "l", 
       lty  = "dotted")

points(x    = 20:40, 
       y    = Qmetabolism_from_mass_temp(m     = 5, 
                                         T_b   = 20:40,
                                         taxon = "amphibian"), 
       type = "l", 
       lty  = "dotdash")

points(x    = 20:40, 
       y    = Qmetabolism_from_mass_temp(m     = 5, 
                                         T_b   = 20:40,
                                         taxon = "invertebrate"), 
       type = "l", 
       lty  = "twodash")

legend(x      = "topleft", 
       title  = "taxa", 
       legend = c("bird,mammal", "reptile", "amphibian", "invertebrate"), 
       lty    = c("solid", "dotted", "dotdash", "twodash"))

```

# Evaporation, Qevap

We provide functions to estimate heat loss associated with evaporation. We offer an empirically-derived relationship for a lizard [@porter1973behavioral; in @gates2012biophysical] and a function based on resistances for amphibians [@spotila1992biophysics]. 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:

```{r, fig.height = 4, fig.width = 4}

# 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 [@spotila1992biophysics]. We provide the relationship for $r_e$ as a function:

```{r, fig.height = 3, fig.width = 4}

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:

```{r, fig.height = 3, fig.width = 4}

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)$.

```{r, fig.height = 4, fig.width = 4}

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 [-@campbell2012introduction] 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, [@gates2012biophysical]] 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 [-@campbell2012introduction] 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 [@mitchell1976heat]:

$$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}$):

```{r, fig.height = 4, fig.width = 4}

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 [-@buckley2007linking]. The function follows the modeling approach contained in Campbell and Norman [-@campbell2012introduction].

### 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 [@forsythe1995model]: 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 [@mitchell1976heat]. A factor of 1.4 is introduced to account for the increased convection due to environmental turbulence [@mitchell1976heat].  

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 [@liu1960interrelationship], $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 [@gates2012biophysical]. 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 [-@barlett1967energy].

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* [@muth1977thermoregulatory]. 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. [-@roughgarden1981resource] provide an expression for total area based on empirical lizard data from Norris [-@norris1967color] and Porter and James [-@porter1979behavioral]: $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 [-@barlett1967energy]. 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:

```{r, fig.height = 4, fig.width = 6}

# 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. [-@fei2012body] 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>.

```{r, fig.height = 4, fig.width = 4}

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 [-@kingsolver1983thermoregulation] 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 [-@buckley2012demographic]. 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:

```{r, fig.height = 4, fig.width = 4}

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 [@buckley2014phenotypic], partially adapted from [@samietz2005altitudinal]. 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 [@anderson1979habitat].

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 [@kingsolver1983thermoregulation]:

$$ 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$, [@anderson1979habitat]], 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 [@samietz2005altitudinal]. 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 [-@lactin1998environmental]. 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. [-@erbs1982estimation].

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}$, [@swinbank1963long]]. 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})$ [@lactin1998environmental]: $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. [-@lactin1998environmental] 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; [@galushko2005electrical]]; $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$; [@galushko2005electrical]]. 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:

```{r, fig.height = 4, fig.width = 4}

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 [-@helmuth1998intertidal] 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 [@helmuth1998intertidal]. 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 [@helmuth1998intertidal]. $Q_{rad,sky}$ and $Q_{rad,ground}$ can be calculated as follows [@helmuth1998intertidal]. 

$$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 [@campbell2012introduction], $\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 [@helmuth1998intertidal]. 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}$ [@denny2006hot]. $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 [@helmuth1998intertidal].

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 [@helmuth1998intertidal]. 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 [-@helmuth1998intertidal] 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.

```{r, fig.height = 4, fig.width = 4}

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 [-@helmuth1998intertidal]. The red line indicates the effect of evaporative heat loss.

```{r, fig.height = 4, fig.width = 4}

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 [@helmuth1999thermal]. The red line indicates the effect of evaporative heat loss.

```{r, fig.height = 4, fig.width = 4}

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 [@mislan2015biophysical]. 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 [-@denny2006hot] 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 [@pennell1989computing]. $\alpha_{sw}$ was empirically acquired to be 0.68 for limpets [@denny2006hot].

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 [@campbell2012introduction].

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 [-@denny1993air] 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$ [@denny2006hot].

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:

```{r, fig.height = 4, fig.width = 4}

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:

```{r, fig.height = 4, fig.width = 4}

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. [-@riddell2018resistance] and Campbell and Norman [-@campbell2012introduction] 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 [@campbell2012introduction].

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}$ [@campbell2012introduction], $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 [@campbell2012introduction], and for $D$ the cross-sectional diameter of a cylinder approximating the size of a salamander. Boundary layer resistance can be estimated as follows.

```{r, fig.height = 4, fig.width = 4}
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  [-@campbell2012introduction], 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 [-@campbell2012introduction] for estimating soil temperatures at specified depths.

```{r, fig.height = 4, fig.width = 6}

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).  

```{r, fig.height = 4, fig.width = 4}

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"))

```

# Example operative environmental temperature calculations

Suppose we want to estimate the hourly body temperature (operative environmental temperature, $T_e$) of a *Sceloporus* lizard on June 1, 2021, in Santa Fe, New Mexico, USA (35.69, -105.944, elevation: 2121m). Let's assume the lizard is in an unshaded location where the daily air temperature varies from a minimum of 10°C to a maximum of 25°C, the soil surface temperature varies from a minimum of 15°C to a maximum of 30°C, and the windspeed is 0.5 m/s. We assume that atmospheric transmissivity $\tau = 0.7$ and albedo $\rho= 0.6$.

We will initially estimate $T_b$ by balancing an energy balance for steady-state conditions. We then use the $T_b$ estimates to explore the components of an energy balance.

## Environmental data

Let us first prepare the environmental data for analysis. We will need to estimate hourly air and soil temperatures and radiation. We start by estimating day of year and the timing of sunrise and sunset:

```{r}

# set up input data as variables
# latitude and longitude (degrees) and elevation (m)

  lat  <- 35.69
  lon  <- -105.944
  elev <- 2121 

# minimum and maximum of air and soil temperatures, respectively (C)

  Tmin   <- 10
  Tmax   <- 25 
  Tmin_s <- 15 
  Tmax_s <- 30 

# wind speed (m/s)

  u <- 0.5 

# assumptions
#  atmospheric transmissivity and albedo

  tau <- 0.7
  rho <- 0.6 
 
# initial assumption of body temperature (C)

  Tb0 <- 25

  doy   <- day_of_year("2021-06-01", 
                       format= "%Y-%m-%d") 
  snoon <- solar_noon(lon = lon, 
                      doy = doy) 
  dayl  <- daylength(lat = lat,
                     doy = doy) 

# time of sunrise

  tr <- snoon - dayl / 2 

# time of sunset

  ts <- snoon + dayl / 2 

```

While measured solar radiation is preferable if available, we can estimate hourly solar radiation by discounting incoming solar radiation as it moved through the atmosphere as follows:

```{r, fig.height = 4, fig.width = 6}

psi_deg <- sapply(0:23, 
                  FUN = zenith_angle, 
                  doy = doy, 
                  lat = lat, 
                  lon = lon)

psi_rad <- degrees_to_radians(psi_deg)

Srad <- sapply(psi_rad, 
               FUN  = solar_radiation, 
               doy  = doy, 
               tau  = tau, 
               elev = elev, 
               rho  = rho)

# Separate into direct, diffuse, and reflected solar radiation

  Sdir <- Srad[1,]
  Sdif <- Srad[2,]
  Sref <- Srad[3,]
 
plot(x    = 0:23, 
     y    = Sdir, 
     type = "l", 
     xlab = "hour", 
     ylab = expression(radiation ~ (W/m^{2})), 
     ylim = c(0,1200))

points(x    = 0:23, 
       y    = Sdif, 
       type = "l", 
       lty  = "dotted")

points(x    = 0:23, 
       y    = Sref, 
       type = "l", 
       lty  = "dashed")

legend(x      = "topright", 
       title  = "Solar radiation component", 
       legend = c("direct", "diffuse", "reflected"), 
       lty    = c("solid", "dotted", "dashed"))

```

We can then calculate hourly air and soil temperatures. We use the sine-exponential model for air temperature and the sine model for surface temperature:

```{r}

Ta <- sapply(0:23, 
             FUN   = diurnal_temp_variation_sineexp, 
             T_max = Tmax,  
             T_min = Tmin,  
             t_r   = tr,  
             t_s   = ts,  
             alpha = 2.59,  
             beta  = 1.55,  
             gamma = 2.2) #air temperature (C)

Ts <- sapply(0:23,  
             FUN   = diurnal_temp_variation_sine,  
             T_max = Tmax_s,  
             T_min = Tmin_s) #soil temperature (C)

```

## Operative temperature estimation

We will first use functions based on energy balances to estimate body temperature given the environmental conditions. We first need to estimate several energy balance parameters. We first estimate the
surface area exposed to radiation We model a 10g *Sceloporus* lizard with solar absorptivity $a = 0.9$ [@gates2012biophysical].

```{r}

mass <- 10 #grams

# solar absorptivity
 
  a <- 0.9 

# estimate surface area (m^2) from mass (g)

  A <- surface_area_from_mass(m     = mass, 
                              taxon = "lizard")

# estimate projected (silhouette) area as a portion of surface area

  psa <- sapply(psi_deg, 
                FUN     = proportion_silhouette_area, 
                taxon   = "lizard", 
                posture = "elevated")

# change negative values to zero

  psa[psa < 0] <- 0

# proportion of surface area in contact with ground, assume 0.25   

  psa_g= 0.25  

# Total radiation 

  Qabs <- psa * Sdir + 0.5 * Sdif + (0.5-psa_g) * Sref
 
```

We then estimate the lizard's heat transfer coefficient, $H_L$ using an empirical relationship for lizards (`heat_transfer_coefficient()`), but also illustrate a function estimating $H_L$ using a spherical approximation (`heat_transfer_coefficient_approximation()`) and a simplified approximation (`heat_transfer_coefficient_simple()`).

```{r}

# Use DRYAIR from NicheMapR to estimate the thermal conductivity of air and kinematic viscosity
# estimate Barometric pressure (pascal)

  ap <- airpressure_from_elev(elev) * 1000 
 
  DRYAIRout <- DRYAIR(db  = Ta, 
                      bp  = ap, 
                      alt = elev)

  # thermal conductivity (Wm^-2K^-1) 
  
    K   <- DRYAIRout$thcond 
  
  # kinematic viscosity (m2 s-1)  
    
    nu  <- DRYAIRout$viskin 

  # approximate snout vent length (m) for Sceloporus

    svl <- 0.006 

  # approximate characteristic dimension as cube root of volume, assuming density of water as 1000kg/m^3
    
    D <- ((mass / 1000) / 1000) ^ (1 / 3)
    
  # We will use the average of K and nu across the day for simplicity and since there's not a substantial variation

    K  <-  mean(K)
    nu <- mean(nu)

  # Estimate the heat transfer coefficient using an empirical relationship for lizards

    H_L <- heat_transfer_coefficient(u     = u,
                                     D     = D,
                                     K     = 25.7 * 10^(-3), 
                                     nu    = 15.3 * 10^(-6), 
                                     taxon = "lizard_surface")

  # Also illustrate estimations using a spherical approximation and a simplified version of the approximation.

    heat_transfer_coefficient_approximation(u     = u,
                                            D     = D,
                                            K     =  25.7 * 10^(-3), 
                                            nu    = 15.3 * 10^(-6), 
                                            taxon = "lizard")

    heat_transfer_coefficient_simple(u    = 0.5,
                                     D    = D, 
                                     type = "Spotila")

```

We now implement the first energy balance.

```{r}

epsilon_s <- 0.965

TeGates <- rep(NA, 24)

for (hr in 1:24) {
  TeGates[hr] <- Tb_Gates(A       = A, 
                          D       = D, 
                          psa_dir = psa[hr], 
                          psa_ref = 0.5 -psa_g, 
                          psa_air = 0.5, 
                          psa_g   = psa_g, 
                          T_g     = Ts[hr], 
                          T_a     = Ta[hr], 
                          Qabs    = Qabs[hr]*A, 
                          epsilon = epsilon_s, 
                          H_L     = H_L, 
                          ef      = 1.3, 
                          K       = 0.5)
}

```

We also implement a similar but simplified energy balance. The energy balance omits conduction with the ground and thus underestimates temperatures, but we implement it for comparison purposes:

```{r}

TeCN <- rep(NA, 24)

for (hr in 1:24) {
  TeCN[hr] <- Tb_CampbellNorman(T_a     = Ta[hr], 
                                T_g     = Ts[hr], 
                                S       = Qabs[hr], 
                                a_l     = 0.96, 
                                epsilon = epsilon_s, 
                                c_p     = 29.3, 
                                D       = D, 
                                u       = u)
}

#S is solar radiation flux (W m^-2), so we divide by surface area, A

```

We additionally estimate $T_e$ using the specialized function for lizards:

```{r}

TeLiz <- rep(NA, 24)

for (hr in 1:24) {
  TeLiz[hr] <- Tb_lizard(T_a       = Ta[hr], 
                         T_g       = Ts[hr], 
                         u         = u, 
                         svl       = svl * 1000, 
                         m         = mass, 
                         psi       = psi_deg[hr], 
                         rho_s     = rho, 
                         elev      = elev, 
                         doy       = doy, 
                         sun       = TRUE, 
                         surface   = TRUE, 
                         a_s       = a, 
                         a_l       = 0.965, 
                         epsilon_s = epsilon_s, 
                         F_d       = 0.8, 
                         F_r       = 0.5, 
                         F_a       = 0.5, 
                         F_g       = 0.5)
}

TeFei <- rep(NA, 24)

for (hr in 1:24) {
  TeFei[hr] <- Tb_lizard_Fei(T_a       = Ta[hr], 
                             T_g       = Ts[hr], 
                             S         = Qabs[hr], 
                             lw        = 30, 
                             shade     = 0.0, 
                             m         = 10.5,
                             Acondfact = 0.05,
                             Agradfact = 0.4)
}

```

We next plot the diurnal trends in environmental data and the operative data.

```{r, fig.height = 4, fig.width = 6}

par(mar = c(4, 4, 1, 4), 
    mgp = c(2, 1, 0), 
    las = 0)

# Gates

  plot(x    = 1:24, 
       y    = TeGates, 
       type = "l", 
       xlab = "Hour", 
       ylab = "Temperature (°C)", 
       col = "blue", 
       ylim = c(10,50)) 

# Campbell

  points(x    = 1:24,  
         y    = TeCN, 
         type = "l", 
         col  = "blue", 
         lty  = "dotted") 

# Buckley 2008

  points(x    = 1:24,  
         y    = TeLiz, 
         type = "l", 
         col  = "blue", 
         lty  = "dashed") 

# Fei et al.

  points(x    = 1:24,  
         y    = TeFei, 
         type = "l", 
         col  = "blue", 
         lty  = "dotdash") 

points(x    = 1:24,  
       y    = Ta, 
       type = "l", 
       col  = "orange")

points(x    = 1:24,  
       y    = Ts, 
       type = "l", 
       col  = "purple")
 
# add additional axis with radiation

  par(new = TRUE)

  plot(x    = 1:24, 
       y    = Sdir, 
       pch  = 16, 
       axes = FALSE, 
       xlab = NA, 
       ylab = NA, 
       type = "l", 
       col  = "green")

  axis(side = 4)

  mtext(side = 4, 
        line = 2, 
        text = "Radiation", 
        col  = "green")

legend(x      = "topleft",
       legend = c("Ta", "Ts", "Tb_Gates", "Tb_CampbellNorman", "Tb_lizard", "Tb_lizard_Fei"),
       lty    = c("solid", "solid", "solid", "dotted", "dashed", "dotdash"), 
       pch    = NA, 
       col    = c("orange", "purple", "blue", "blue", "blue", "blue"))

```

## Energy balance

We will now use body temperature estimates to examine the components of the energy balance. We will be solving the following energy balance:

$$Q_{net} = Q_{abs} - Q_{emit} - Q_{conv} - Q_{cond} + Q_{met} - Q_{evap}.$$

We will estimate each term on the right side of the equation in turn. We calculate the hourly radiation absorbed as follows:

```{r}

Qabs <- rep(NA, 24)

for (hr in 1:24) {
  Qabs[hr] <- Qradiation_absorbed(a       = a, 
                                  A       = A, 
                                  psa_dir = psa[hr],
                                  psa_dif = 0.5,
                                  psa_ref = 0.5-psa_g, 
                                  S_dir   = Sdir[hr], 
                                  S_dif   = Sdif[hr], 
                                  rho     = rho)
}

```

We estimate thermal radiation $Q_{emit}$ for the lizard outdoors. We assume the surface emissivity of lizards, $\epsilon_s = 0.965$ [@barlett1967energy].

```{r}

# Use Gates model as Te estimate

  Te <- TeGates

Qemit <- rep(NA, 24)

for (hr in 1:24) {
  Qemit[hr] <- Qemitted_thermal_radiation(epsilon  = epsilon_s, 
                                          A        = A, 
                                          psa_dir  = 0.5, 
                                          psa_ref  = 0.5-psa_g, 
                                          T_b      = Te[hr] + 273, 
                                          T_g      = Ts[hr] + 273, 
                                          T_a      = Ta[hr] + 273, 
                                          enclosed = FALSE)
}

```

We next estimate convection $Q_{conv}$ and conduction $Q_{cond}$. We estimate convective heat exchange between the animal and surrounding air using the following relationship:

$$Q_{conv}= ef\cdot H_L(A\cdot \mbox{proportion})(T_{aK}-T_{bK}),$$ 

where an enhancement factor, $ef$, multiplier can be incorporated to account for increases in heat exchange resulting from air turbulence in field conditions.

The function is available in R assuming that 2/3 of the lizard's surface area is exchanging heat through convection.

```{r}

Qconv <- rep(NA, 24)

for (hr in 1:24) {
  Qconv[hr] <- Qconvection(T_a        = Ta[hr] + 273.15,
                           T_b        = Te[hr] + 273.15,
                           H          = H_L,
                           A          = A, 
                           proportion = 0.67, 
                           ef         = 1.3)
}

```

We estimate conductance between the lizard and surface assuming conductance through the animal tissue is the rate limiting step as follows: 

$$Q_{cond}= A\cdot \mbox{proportion}\cdot K(T_{gK}-T_{bK})/d.$$ 

We conduct the estimate in R assuming that conductive heat exchange occurs down to a soil depth of 2.5cm. We use this value rather than skin thickness, which results in rapid conduction and does not readily reach steady state conditions.

```{r}

Qcond <- rep(NA, 24)

for (hr in 1:24) {
  Qcond[hr] <- Qconduction_animal(T_g        = Ts[hr] + 273,
                                  T_b        = Te[hr] + 273,
                                  d          = 0.025, #assume for conductive heat exchange occurs down to 2.5cm
                                  K          = 0.5,
                                  A          = A, 
                                  proportion = psa_g)
}

```

We assume, as in generally done for lizards, that heat exchange associated with metabolism and evaporation is negligible, but we note that functions for estimating each heat exchange are described above.

```{r}

Qmet  <- 0
Qevap <- 0

```

The full heat budget can be calculated as follows. We plot out the energy balance components.

```{r, fig.height = 4, fig.width = 7}

Qnet <- Qnet_Gates(Qabs  = Qabs, 
                   Qemit = Qemit, 
                   Qconv = Qconv, 
                   Qcond = Qcond, 
                   Qmet  = Qmet, 
                   Qevap = Qevap)

par(mar = c(5, 5, 3, 5))

plot(x    = 1:24, 
     y    = Qnet, 
     type = "l", 
     xlab = "hour", 
     ylab = expression("heat flux" ~ (W/m^{2})), 
     ylim = c(-1,4))

points(x    = 1:24, 
       y    = Qabs, 
       type = "l", 
       lty  = "dotted")

points(x    = 1:24, 
       y    = Qemit, 
       type = "l", 
       lty  = "dashed")

points(x    = 1:24, 
       y    = Qconv, 
       type = "l", 
       lty  = "dotdash")

points(x    = 1:24, 
       y    = Qcond, 
       type = "l", 
       lty  = "twodash")

legend(x      = "topright", 
       title  = "component", 
       legend = c("net radition, Qnet", "solar radiation, Qabs", "thermal radiation, Qemit", "convection, Qconv", "conduction, Qcond"), 
       lty    = c("solid", "dotted", "dashed", "dotdash", "twodash"))

```

```{r include=FALSE}

suppressWarnings(par(oldpar))

```

# References
