TrenchR_TeExample.Rmd

---
title: 'TrenchR: an R package for modular and accessible microclimate and biophysical ecology'
author: 
  - "Lauren B. Buckley, Bryan A. Briones Ortiz, Isaac Caruso, Aji John, Ofir Levy,"
  - "Abigail V. Meyer, Eric A. Riddell, Yutaro Sakairi, and Juniper Simonis"
output:
  word_document: default
  html_document: default
  pdf_document: 
    pandoc_args: --listings
    includes:
      in_header: preamble.tex
bibliography: bibliography.bib
csl: plos.csl
editor_options: 
  markdown: 
    wrap: 72
---

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

library(TrenchR)

```

We illustrate the use of the TrenchR package by estimating an energy budget for a *Sceloporus* lizard on June 1, 2021 in Santa Fe, New Mexico, USA (35.69$^\circ$N, -105.944$^\circ$W, elevation: 2121 m). 
The simplified example, which is designed to be self-contained, is also incorporated in the *Using energy balances to estimate body temperatures* vignette. 
We start by generating environmental inputs .
Using these inputs, we estimate the energy budget with component functions. 
Finally, we use an integrated biophysical model to estimate operative environmental temperatures. 

Let us assume the lizard is in an unshaded location where a weather station at standard height (2 meters) reports that the daily air temperature varies from a minimum of 10 $^\circ$C to a maximum of 25 $^\circ$C and the wind speed averages 1 m/s. 
The soil surface temperature varies from a minimum of 15 $^\circ$C to a maximum of 30 $^\circ$C. 
We assume that atmospheric transmissivity $\tau = 0.7$ and albedo $\rho = 0.6$.

### Environmental data

At the first stage, we prepare the environmental data for analysis. 
We will estimate hourly air and soil temperatures and radiation using a function describing diurnal temperature variation. 
We start by estimating the day of year and the timing of sunrise and sunset:

```{r}
# Set up input data as variables

lat  <- 35.69    # Latitude (degrees)
lon  <- -105.944 # Longitude (degrees)
elev <- 2121     # Elevation (meters)

Tmin   <- 10 # Minimum air temperature (C)
Tmax   <- 25 # Maximum air temperature (C)
Tmin_s <- 15 # Minimum soil temperature (C)
Tmax_s <- 30 # Maximum soil temperature (C)
u      <- 1  # Wind speed (m/s)

# Assumptions

tau <- 0.7 # Atmospheric transmissivity
rho <- 0.6 # Albedo
Tb0 <- 25  # Initial assumption of body temperature (C)

doy   <- day_of_year("2021-06-01", format = "%Y-%m-%d") # Day of year

snoon <- solar_noon(lon = lon, doy = doy) # Estimate solar noon
dayl  <- daylength(lat = lat, doy = doy)  # Estimate day length

tr    <- snoon - dayl / 2 # Time of sunrise
ts    <- snoon + dayl / 2 # Time of sunset
```

Although measured solar radiation is preferable if available, we can estimate hourly solar radiation by discounting incoming solar radiation as it moves through the atmosphere as follows. We use the approach from Cambpell and Norman [-@campbell2000introduction], which uses an empirical relation to partition radiation into direct, diffuse, and reflected components. The `partition_solar_radiation()` function includes 8 empirical relationships for, and the proportion_diffuse_solar_radiation() includes a more complex numerical approximation for, partitioning radiation components as described in the *Estimating microclimates* vignette.

```{r}
# Estimate zenith angle (degrees)

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

# Convert to radians

psi_rad <- degrees_to_radians(psi_deg) 

# Estimate radiation

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

# Separate solar radiation into direct, diffuse, and reflected components

Sdir <- Srad[1,] # Direct solar radiation (W/m2)
Sdif <- Srad[2,] # Diffuse solar radiation  (W/m2)
Sref <- Srad[3,] # Reflected solar radiation  (W/m2)

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 then calculate hourly air and soil surface temperatures based on daily minimum and maximum temperatures. 
We select the sine-exponential model for air temperature and the sine model for surface temperature [@wann1985evaluation]:

```{r}
# Air temperature (C)

Ta <- sapply(1:24, 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) 

# Soil surface temperature (C)

Ts <- sapply(1:24, diurnal_temp_variation_sine, T_max = Tmax_s, T_min = Tmin_s) 
```

At the second stage, we use microclimate models to scale air temperature ($T_r$) and wind speed ($u_r$) from weather station height (reference height $z_r$= 2 m) to lizard height (organism height $z$= 0.02 m). 
We assume a surface roughness of $z_0 = 0.2$ m, which corresponds to bare sand and determines the turbulence of airflow. 
We implement free air temperature and wind speed profiles driven by density differences but profiles forced by wind speed are also available.

```{r}
# Neutral air temperature profile

Ta_liz  <- air_temp_profile_neutral(T_r = Ta, zr = 2, z0 = 0.2, z = 0.02, T_s = Ts) 

# Neutral wind speed profile

u_liz <- wind_speed_profile_neutral(u_r = u, zr = 2, z0 = 0.2 , z = 0.02) 

plot(x    = 0:23, 
     y    = Ta, 
     type = "l", 
     xlab = "hour", 
     ylab = "temperature", 
     ylim = c(8,32))

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

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

legend(x      = "topright", 
       title  = "Temperature", 
       legend = c("Ta", "Ts", "Ta liz"), 
       lty    = c("solid", "dotted", "dashed"))
```

### Energy balance

Finally, we will use our microclimates estimates to solve the following energy balance to estimate $T_e$:

$$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 generated by metabolism (W), and $Q_{evap}$ is the energy generated by evaporative water loss (W). 
We will estimate each term on the right side of the equation in turn. 
Estimating $Q_{abs}$ requires the surface area exposed to radiation and the solar absorptivity of the animal surface ($a$ proportion). We use zenith angle $psi$ to estimate the projected (silhouette) area as a portion of the surface area of the organism, which allows estimating absorbed solar radiation.
We model a 10 gram *Sceloporus* lizard with solar absorptivity $a = 0.9$ [@gates1980biophysical]. 
We will initially assume $T_b = T_a + 10$ to illustrate the calculations before solving for $T_b$ given the environmental conditions.

```{r}
mass <- 10  # Mass (g)
svl  <- 0.006 # Snout vent length (meters)
a    <- 0.9 # Solar absorptivity (proportion)

#assume 1/3 of surface area is in contact with surface
psa_g <- 0.33

# Estimate surface area (m^2) and the proportion sihouette area

A   <- surface_area_from_mass(mass, "lizard")  
psa <- sapply(psi_deg, proportion_silhouette_area, taxon = "lizard", posture = "prostrate") 

# Change negative values to zero
psa[psa < 0] = 0
```

We calculate the hourly solar and thermal radiation absorbed (W) as follows:

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

}
```

We estimate thermal radiation $Q_{emit}$ (W) for the lizard outdoors, where $psa_{dir}$ and $psa_{ref}$ are the view factors, also refered to as configuration factors, that indicate the proportions of surface area $A$ ($m^2$) exposed to the sky and ground, respectively. We assume the surface emissivity of lizards, $epsilon_s = 0.965$ [@barlett1967].

```{r}
epsilon_s <- 0.965 # Surface emissivity of lizards

Qemit <- rep(NA, 24)
for (hour in 1:24) {
  
  Qemit[hour] <- Qemitted_thermal_radiation(epsilon = epsilon_s, A = A, psa_dir = 0.5, psa_ref = 0.5, T_b = Ta_liz[hour] + 273.15, T_g = Ts[hour] + 273.15, T_a = Ta_liz[hour] + 273.15, enclosed = FALSE)
  
}
```

We next estimate convection $Q_{conv}$ (W) and conduction $Q_{cond}$ (W). 
We will estimate the lizard's heat transfer coefficient, $H_L$ ($Wm^{-2}K^{-1}$) using an empirical relationship for lizards (`heat_transfer_coefficient()`). We average thermal conductivity and kinematic viscosity across the day for simplicity and since there is not substantial diurnal variation. We also illustrate a function estimating $H_L$ using a spherical approximation (`heat_transfer_coefficient_approximation()`) and a simplified approximation (`heat_transfer_coefficient_simple()`) for cases when taxon specific relationships for estimating heat transfer coefficients are not available. We estimate the characteristic dimension, which determines exposure to convective heat exchange as the cube root of volume, assuming the animal density approximates that of water [@mitchell1976heat]. 

These coefficients assume convection is forced by the wind. TrenchR includes approaches for free convection and a function (`free_or_forced_convection()`) that evaluates whether free or forced convection is appropriate. The function uses dimensionless numbers, which have been developed to describe heat transfer coefficients associated with convection over different geometries and can be estimated using TrenchR (e.g., Grashof, Nusselt, and Reynolds numbers).

```{r}
# Use DRYAIR from NicheMapR to estimate the thermal conductivity of air and kinematic viscosity. 

ap <- airpressure_from_elev(elev) * 1000 # Barometric pressure (pascal)

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

K  <- mean(DRYAIRout$thcond) # Thermal conductivity (Wm^-2K^-1) 
nu <- mean(DRYAIRout$viskin) # Kinematic viscosity (m2 s-1)

# Estimate the characteristic dimension as cube root of volume, assuming density of water as 1000kg/m^3
    
D <- ((mass / 1000) / 1000) ^ (1 / 3)

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

H_L <- heat_transfer_coefficient(u = u_liz, D = D, K = K, nu = nu, taxon = "lizard_surface")

# Estimate the heat transfer coefficient using a spherical approximation

H_L2 <- heat_transfer_coefficient_approximation(u = u_liz, D = D, K = K, nu = nu, taxon = "lizard")

# Estimate the heat transfer coefficient using a simplified version of the approximation

H_L3 <- heat_transfer_coefficient_simple(u = u_liz, D = svl, type = "Gates")
```

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_a-T_b),$$ 

where an enhancement factor, $ef$, multiplier can be incorporated to account for increases in heat exchange resulting from air turbulence in field conditions. 
We implement the function in R assuming that 2/3 of the lizard's surface area is exchanging heat through convection. 

```{r}
Qconv <- rep(NA, 24)
for (hour in 1:24) {
  
  Qconv[hour] <- Qconvection(T_a = Ta_liz[hour] + 273.15, T_b = Ta_liz[hour] + 10 + 273.15, H = H_L, A = A, proportion = 0.67, ef = 1.3)
  
}
```

We estimate conductive heat flow (W) from the lizard to the surface assuming conductance through the animal tissue is the rate limiting step as follows:

$$Q_{cond} = A \cdot \mbox{proportion} \cdot K_{skin}(T_g-T_b)/d$$ 
where $K_{skin}$ in the thermal conductivity of lizard skin ($Wm^-2K^-1$). We implement the estimate 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.15, T_b = Ta_liz[hr] + 10 + 273.15, d = 0.025, K = 0.5, A = A, proportion = psa_g)
  
}
```

We assume, as is generally done for lizards, that heat exchange associated with metabolism and evaporation is negligible. However, functions for estimating both forms of heat exchange available in TrenchR.

```{r}
Qmet  <- 0 
Qevap <- 0
```

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

```{r}
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"))

```

We now use a function based on the Gates energy balance above to estimate body temperature given the environmental conditions:

```{r}
Te <- rep(NA, 24)
for (hour in 1:24) {
  
  Te[hour] <- Tb_Gates(A = A, D = svl, psa_dir = psa[hour], psa_ref = 1 - psa[hour], psa_air = 0.67, psa_g = 0.25, T_g = Ts[hour], T_a = Ta_liz[hour], Qabs = Qabs[hour], epsilon = epsilon_s, H_L = H_L, ef = 1.3, K = K)
  
}
```

We also implement a similar but simplified energy balance [@campbell2000introduction]. 
The energy balance omits conduction with the ground:

```{r}
Te2 <- rep(NA, 24)
for (hr in 1:24) {
  
  # S is solar radiation flux (W m^-2), so we divide by surface area, A

  Te2[hr] <- Tb_CampbellNorman(T_a = Ta_liz[hr], T_g = Ts[hr], S = Qabs[hr] / A, a_l = 0.96, epsilon = epsilon_s, c_p = 29.3, D = D, u = u_liz)
  
}

```

We additionally estimate $T_b$ using a specialized function for lizards [@buckley2008link], where $F_d$, $F_r$, $F_a$, and $F_g$ are the view factors between the surface of the lizard and diffuse solar radiation, reflected solar radiation, atmospheric thermal radiation, and ground thermal radiation, respectively:

```{r}
Te3 <- rep(NA, 24)
for (hour in 1:24) {
  
  Te3[hour] <- Tb_lizard(T_a = Ta_liz[hour], T_g = Ts[hour], u = u_liz, svl = svl * 1000, m = mass, psi = psi_deg[hour], 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)

}
```

We then plot a comparison of operative temperature estimates. 

```{r, echo=FALSE}
par(mar = c(5, 4, 4, 8) + 0.1)

plot(1:24, Te, type = "l", xlab = "Hour", ylab = expression("Temperature ("*~degree*C*")"), col = "#5DC863FF", ylim = c(10,70), las = 1) # Gates
points(1:24, Te2, type ="l", col = "#5DC863FF", lty = "dotted") # Campbell and Norman
points(1:24, Te3, type = "l", col = "#5DC863FF", lty = "dashed") # Lizard, Buckley 2008
points(1:24, Ta, type = "l", col = "#440154FF")
points(1:24, Ta_liz, type = "l", col = "#3B528BFF")
points(1:24, Ts, type = "l", col = "#21908CFF")

# add additional axis with radiation

par(new = TRUE)
plot(1:24, Sdir, pch = 16, axes = FALSE, xlab = NA, ylab = NA, type = "l", col = "#FDE725FF")
axis(side = 4, las = 1)
mtext(side = 4, line = 3, 'Radiation (W/m^2)')
legend("topleft", bty = "n",
       legend = c("Ta (2m)", "Ta (0.02m)", "Ts", "Te", "Radiation"),
       lty = 1, pch = NA, col = c("#440154FF", "#3B528BFF", "#21908CFF", "#5DC863FF", "#FDE725FF"))
```