Biomass Script

Time series needed

Necessary
column 1: t
column 2: water temperature
column 3: solar irradiance
column 4: nitrogen concentration
column 5: phosphorous concentration

Potential further columns
column 6: silica concentration
column 7: carbon concentration
column 8: secchi disk depth

Global object definitions

Ecological constants
PN = N half-saturation constant
PP = P half-satruation constant
PSi = Si half-saturation constant
PC = C half-saturation constant
Is = saturating light intensity
Tmult, Kmult = T(1-4), K(1-4)
Gmax = Maximum growth rate
Rmax = Maximum respiration rate
max = Maximum excretion rate
Emax = Maximum mortality rate

Reading in data

# Input data frame containing necessary time series (currently using fake data that encompasses values above and below optimal)
df <- data.frame(t = c(1:50),
                 temp = c(1:50),
                 irradiance = seq(30, 128, 2),
                 nitrogen = seq(0, 0.147, 0.003),
                 phosphorous = seq(0,0.0245, 0.0005))

Setting constant values

# Initiate biomass list and set initial biomass (units: g/m^3)
density <- 45

# Set timestep (units: day)
timestep <- 1

# Set a space to experiment with different constant values. Replacing NA with a number will automatically make the algorithm use that number
customs <- list(PN = NA,
                PP = NA,
                Is = NA,
                Ti = c(NA, NA, NA, NA),
                Ki = c(NA, NA, NA, NA),
                Gmax = NA,
                Rmax = NA,
                Emax = NA,
                Mmax = NA)

# Set default constant values (currently the same as CE-QUAL-W2 defaults, found on pg.190 and pg.203 of user manual part 3)
defaults <- list(PN = 0.014,
            PP = 0.003,
            Is = 100,
            Ti = c(5, 25, 35, 40),
            Ki = c(0.1, 0.99, 0.99, 0.1),
            Gmax = 2,
            Rmax = 0.04,
            Emax = 0.04,
            Mmax = 0.1)

# Set constants list. Use custom values where present; otherwise use defaults.
constants <- defaults
for (i in 1:length(customs)) {
  non_na <- !is.na(customs[[i]])
  constants[[i]][non_na] <- customs[[i]][non_na]
}

# Make each column in the constants list its own object
list2env(constants, envir = .GlobalEnv)

## <environment: R_GlobalEnv>

Create functions

# Function to get temperature multiplier
TempMultCalc <- function(Temp, Ti, Ki){
  Y1 <- (1 / (Ti[2] - Ti[1])) * log( (Ki[2]/Ki[1]) * ((1-Ki[1]) / (1-Ki[2])) )
  Y2 <- (1 / (Ti[4] - Ti[3])) * log( (Ki[3]/Ki[4]) * ((1-Ki[4]) / (1-Ki[3])) )
  
  Yar <- (Ki[1] * exp(Y1*Temp - Y1*Ti[1])) / (1 + Ki[1]*exp(Y1*Temp - Y1*Ti[1]) - Ki[1])
  Yaf <- (Ki[4] * exp(Y2*Ti[4] - Y2*Temp)) / (1 + Ki[4]*exp(Y2*Ti[4] - Y2*Temp) - Ki[4])
  
  Mult <- Yar * Yaf
  
  return(Mult)
}

# Function to get nutrient multipliers. Based on "Monod Relationship"
Monod <- function(conc, halfsat){
  Mult <- conc / (halfsat + conc)
  return(Mult)
}

Algal biomass algorithm

# Initiate lists to be used in algorithm
TempMult <- NMult <- PMult <- LMult <- MinMult <- Kag <-
  Kar <- Kae <- Kam <- Sa <- vector("double", length = length(df$t))

# Initiate for loop to iterate over timesteps t
for(t in 1:length(df$t)){

# Temperature Multiplier
Temp <- df$temp[t]
  
if (Temp <= Ti[1] | Temp >= Ti[4]){
  TempMult[t] <- 0
} else {
  TempMult[t] <- TempMultCalc(Temp, Ti, Ki)
}

# Nutrient Mutlipliers
N_Conc <- df$nitrogen[t]
NMult[t] <- Monod(N_Conc, PN)

P_Conc <- df$phosphorous[t]
PMult[t] <- Monod(P_Conc, PP)

# Light Multiplier
I <- df$irradiance[t]
LMult[t] <- (I/Is) * exp((-I/Is) + 1)

# Set minimum multiplier
MinMult[t] <- min(NMult[t], PMult[t], LMult[t])

# Growth rate calculations
Kag[t] <- TempMult[t] * MinMult[t] * Gmax

Kar[t] <- TempMult[t] * Rmax

Kae[t] <- TempMult[t] * (1 - LMult[t]) * Emax

Kam[t] <- TempMult[t] * Mmax

# Update biomass time-series

Sa[t] <- density[t] * (Kag[t] - Kar[t] - Kae[t] - Kam[t])
density[t+1] <- Sa[t]*timestep + density[t]

# Finish iteration
}