2  Simulate Habitat

if any changes are made, need to re-run qmd_to_r_script() function at end, in console

Purpose: This script simulates habitat conditions (temperature and food availability) for a simple, spatially-implicit habitat network consisting of two patches.

Notes:

• In theory, habitat network is generalizable and can be expanded to include as many patches as needed. But in doing so we will need to update the simulation to allow more than 2 patches.
• We can specify food available in one of two ways, ration size in g/g/d or proportion of maximum consumption (P_C_max). Fullerton’s IBM in Armstrong et al. (2021) sets a constant ration size (g/g/d) in each patch, and varies this from upstream to downstream. The problem with this is that is that it forces small fish to feed way below C_max, but as fish grow and allometric C_max declines, fish feed much closer to or at C_max. Therefore, growth rates are low and constraining for small fish, but increase as fish grow larger, such that fish are less constrained by food availability.
  • Many researchers using bioenergetics models simply set a constant P_C_max, then estimate growth from that. But the problem with this, as Railsback points out, is that under warmer conditions fish are able to attain higher ration sizes simply because of the temperature dependence function: C = Cmax * p * f(T). Thus, there is no way to separate ration size from habitat temperature using a constant P_C_max.
  • Instead, we can simply modify the temperature dependence function to set some threshold P_C_max. Basically, temperature acts on ration/consumption up to some threshold value of P_C_max, above which body size alone limits consumption (via allometric C_max).

Current configuration:

• Temp: simple warm and cold patches…cold regime is simply 60% of daily warm regime values. No noise.
• Ration: constant food availability, same in each patch. No noise.
• P_Cmax: constant proportion of maximum consumption, same in each patch. No noise.

2.1 Temperature

Simulate thermal regimes as annual sine waves. Cold vs. Warm using simple multiplier.

# Set seed for reproducibility
set.seed(123)

# 1. Create a date range (1 year)
dates <- seq(as.Date("2024-05-01"), as.Date("2034-12-31"), by="day")
n <- length(dates)
doy <- as.numeric(format(dates, "%j")) # Day of year

# 2. Simulate Stream Temperature (Sine wave + Noise)
# Parameters: Midpoint, Amplitude, Offset
base_temp <- 11        # Average annual temperature
amplitude <- 14         # Seasonal fluctuation amplitude
# Peak temperature occurs late in the year (July/Aug)
stream_temp <- base_temp + amplitude * sin(2 * pi * (doy - 100) / 365)

# 3. Add random noise (daily weather variation)
noise <- rnorm(n, mean = 0, sd = 0) # no variation for simplicity
simulated_temp <- stream_temp + noise

# 4. Force values near 0 to small postive, max = 25
simulated_temp[simulated_temp < 0.5] <- 0.5
simulated_temp[simulated_temp > 25] <- 25

# 4. Simulate cold regime using a multiplier
simulated_temp_cold <- simulated_temp*0.6

# 4. Create Data Frame
habitat_df <- tibble(date = dates, 
                     doy = doy, 
                     dayofsim = 1:length(dates),
                     temp_warm = round(simulated_temp, digits = 1), 
                     temp_cold = round(simulated_temp_cold, digits = 1))

# 5. Plot the simulation
summer_shading <- data.frame(
  xmin = as.Date(paste0(unique(format(habitat_df$date, "%Y")), "-06-01")),
  xmax = as.Date(paste0(unique(format(habitat_df$date, "%Y")), "-08-31")),
  ymin = -Inf, ymax = Inf
)

habitat_df %>% ggplot() +
  geom_rect(data = summer_shading,
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = "grey", alpha = 0.3, inherit.aes = FALSE) +
  geom_line(aes(x = date, y = temp_warm), color = "red") +
  geom_line(aes(x = date, y = temp_cold), color = "blue") +
  theme_bw() + 
  xlab("Date") + 
  ylab("Simulated daily stream temperature (deg C)") + 
  labs(title = "Daily stream temperature with Jun-Aug shaded")

# function to retrieve water temperature in patch = patch and at doy/time = t
get_patchtemp <- function(t, patch) {
  if (patch == "warm") {
    habitat_df$temp_warm[habitat_df$dayofsim == t]
  } else {
    habitat_df$temp_cold[habitat_df$dayofsim == t]
  }
}

2.2 Ration

Simulate food availability in each patch, over time.

# define habitat specific ration size
ration_cold <- 0.10
ration_warm <- 0.10

# define random variation in ration size
ration_sd <- 0.00 # no variation for simplicity

# function to generate patch food availability at time t (stochastic)
patch_food <- function(patch) {
  if (patch == "warm") {
    max(0, min(1, rnorm(1, ration_warm, ration_sd)))
  } else {
    max(0, min(1, rnorm(1, ration_cold, ration_sd)))
  }
}

# define food in habitat tibble
habitat_df <- habitat_df %>% 
  mutate(ration_warm = replicate(n = dim(.)[1], patch_food("warm")),
         ration_cold = replicate(n = dim(.)[1], patch_food("cold")))

# plot ration size
habitat_df %>% ggplot() + 
  geom_line(aes(x = date, y = ration_warm), color = "red") +
  geom_line(aes(x = date, y = ration_cold), color = "blue") +
  theme_bw() + 
  xlab("Date") + 
  ylab("Simulated daily ration size")

# function to retrieve ration in patch = patch and at doy/time = t
get_patchration <- function(t, patch) {
  if (patch == "warm") {
    habitat_df$ration_warm[habitat_df$dayofsim == t]
  } else {
    habitat_df$ration_cold[habitat_df$dayofsim == t]
  }
}

2.3 Proportion C_max

Simulate food availability as proportion of C_max (p-value) in each patch, over time.

# define habitat specific ration size
pcmax_cold <- 0.5
pcmax_warm <- 0.5

# define random variation in ration size
pcmax_sd <- 0.00 # no variation for simplicity

# function to generate patch food availability at time t (stochastic)
patch_pcmax <- function(patch) {
  if (patch == "warm") {
    max(0, min(1, rnorm(1, pcmax_warm, ration_sd)))
  } else {
    max(0, min(1, rnorm(1, pcmax_cold, ration_sd)))
  }
}

# define food in habitat tibble
habitat_df <- habitat_df %>% 
  mutate(pcmax_warm = replicate(n = dim(.)[1], patch_pcmax("warm")),
         pcmax_cold = replicate(n = dim(.)[1], patch_pcmax("cold")))

# plot ration size
habitat_df %>% ggplot() + 
  geom_line(aes(x = date, y = pcmax_warm), color = "red") +
  geom_line(aes(x = date, y = pcmax_cold), color = "blue") +
  theme_bw() + 
  xlab("Date") + 
  ylab("Simulated daily P_Cmax")

# function to retrieve ration in patch = patch and at doy/time = t
get_patchpcmax <- function(t, patch) {
  if (patch == "warm") {
    habitat_df$pcmax_warm[habitat_df$dayofsim == t]
  } else {
    habitat_df$pcmax_cold[habitat_df$dayofsim == t]
  }
}

2.4 View data

str(habitat_df)

tibble [3,897 × 9] (S3: tbl_df/tbl/data.frame)
 $ date       : Date[1:3897], format: "2024-05-01" "2024-05-02" ...
 $ doy        : num [1:3897] 122 123 124 125 126 127 128 129 130 131 ...
 $ dayofsim   : int [1:3897] 1 2 3 4 5 6 7 8 9 10 ...
 $ temp_warm  : num [1:3897] 16.2 16.4 16.6 16.8 17.1 17.3 17.5 17.7 17.9 18.1 ...
 $ temp_cold  : num [1:3897] 9.7 9.8 10 10.1 10.2 10.4 10.5 10.6 10.7 10.9 ...
 $ ration_warm: num [1:3897] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
 $ ration_cold: num [1:3897] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...
 $ pcmax_warm : num [1:3897] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...
 $ pcmax_cold : num [1:3897] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

head(habitat_df)

# A tibble: 6 × 9
  date         doy dayofsim temp_warm temp_cold ration_warm ration_cold
  <date>     <dbl>    <int>     <dbl>     <dbl>       <dbl>       <dbl>
1 2024-05-01   122        1      16.2       9.7         0.1         0.1
2 2024-05-02   123        2      16.4       9.8         0.1         0.1
3 2024-05-03   124        3      16.6      10           0.1         0.1
4 2024-05-04   125        4      16.8      10.1         0.1         0.1
5 2024-05-05   126        5      17.1      10.2         0.1         0.1
6 2024-05-06   127        6      17.3      10.4         0.1         0.1
# ℹ 2 more variables: pcmax_warm <dbl>, pcmax_cold <dbl>

2.5 Write R file

For sourcing later

#quarto::qmd_to_r_script("Habitat.qmd")