thermod-Package

thermod is a simple two-layer heat transport model that assumes that the lake is divided into two completely mixed volumes: the epilimnion and the hypolimnion. Both layers are divided by a thermocline zone. The entrainment over the thermocline depends on a diffusion coefficient which is a function of the diffusion at neutral stability to the Richardson number. All equations and derivations are based on or from Steven C. Chapra (2008) ‘Surface Water-Quality Modeling’ Waveland Press, Inc.

thermod allows you to easily configure a lake setup from a LakeEnsemblR-configuration file:

# use a LakeEnsemblR setup
config_file <- 'LakeEnsemblR.yaml'
folder = '.'
parameters <- configure_from_ler(config_file <- config_file, folder = folder)

# load in the meteorological boundary data
bound <-read_delim(paste0(folder,'/meteo.txt'), delim = '\t')
bound <- bound[, -c(1)]
colnames(bound) <- c('Day','Jsw','Tair','Dew','vW')

# function to calculate wind shear stress (and transforming wind speed from km/h to m/s)
bound$Uw <- 19.0 + 0.95 * (bound$vW * 1000/3600)^2 
bound$vW <- bound$vW * 1000/3600

boundary <- bound

thermod then solves its ordnary differential equations on a vertical grid with two elements and on a daily time step using the 4th-order Runge-Kutta method: \(c_{i+1}= c_i + [\frac{1}{6}(k_1+2k_2+k_3+k_4)]h\). Further, the depth of the thermocline is fixed and either provided by the user or estimated based on maximum length or width of the basin using the empirical eqation from Hanna M. (1990) Evaluation of Models predicting Mixing depth. Can. J. Fish. Aquat. Sci. 47: 940-947: \(z_{therm} = 10^{0.336 log10(max(L,W))-0.245}\).

Temperature simulation

thermod has three main modes: (1) temperature, (2) temperature and dissolved oxygen, and (3) temperature, dissolved oxygen and nutrient-food web model (NPZ). Let us first simulate only water temperatures in both layers for Lake Mendota, and compare this with observed data.

Here, thermod solves two ordinary differential equations for water temperature: \(\frac{dT_{epi}}{dt} = \frac{Q}{V_{epi}}T_{in} - \frac{Q}{V_{epi}}T_{epi} + {\upsilon }_t \frac{A_t}{V_{epi}}(T_{hypo} - T_{epi}) \pm \frac{A_{surf}}{V_{epi} \rho c_p} J\) and \(\frac{dT_{hypo}}{dt} = {\upsilon}_t \frac{A_t}{V_{hypo}} (T_{epi} - T_{hypo})\)

The atmospheric heat fluxes, \(J\), are separated into the five sub-heat fluxes shortwave radiation, incominglongwave radiation, outgoing longwace radiation, latent heat flux and sensible heat flux: \(J = J_{sw} + J_{in} * J_{out} + J_{E} + J_{H}\). Here, \(J_{sw}\) is a measured bondary condition, whereas the other heat terms are estimated from meteorological data: \(J_{in} = \sigma (T_{air} + 273)^4 (K+0.031 \sqrt{e_{air}})(1-R_L)\), \(J_{out}=\epsilon \sigma (T_{epi} +273)^4\), \(J_E = f(Uw)(e_{sat}-e_{air})\) with \(f(Uw)= 19.0 + 0.95 {u_{wind}}^2\), and \(J_H = c_1 f(Uw)(T_{epi} - T_{air})\).

thermod aims to get a dynamic representation of the vertical diffusivity in dependence of the water column stability (in contrast to telling the model to use high or low values on specific days). Therefore, the model states the diffusion coefficient as being dependent on the diffusion at neutral stability to the Richardson number: \({\upsilon}_t = \frac{E_0}{(1+\alpha {R_i})^{3/2}}\). The diffusion coefficient at neutral stability then depends on the ratio of wind shear stress to water density: \(E_0 = c \sqrt{\frac{\tau}{\rho_{water}}}\) with wind shear stress as: \({\tau} = {\rho}_{air} C_d {u_{wind}}^2\). Finally, we can estimate the Richardson number, which generally describes the work of buoyancy against wind-induced turbulence. If it is below 1/4, the system will experience shear-induced turbulence: \(Ri = \frac{-\frac{g}{\rho}\frac{d\rho}{dz}}{\sqrt{\frac{\tau}{\rho_{water}}}H^{-2}}\) which acts over a thickness \(H\) which we assume to be the thermocline depth.

# simulation maximum length
times <- seq(from = 1, to = max(boundary$Day), by = 1)

# initial water temperatures
yini <- c(3,3) 

# there is calibration parameter for the vertical diffusivity over the thermocline, which we
# will set slightly higher to 3.2
parameters[19] = 3.2 

# you can also set ice_on to TRUE to reduce atmospheric exchanges during idealized ice conditions
ice_on = TRUE 

out <- run_model(bc = boundary, params = parameters, ini = yini, times = times, ice = ice_on)

result <- data.frame('Time' = out[,1],
                     'WT_epi' = out[,2], 'WT_hyp' = out[,3])
head(result)
##   Time   WT_epi   WT_hyp
## 1    1 3.000000 3.000000
## 2    2 2.903372 2.412467
## 3    3 2.499375 2.234035
## 4    4 2.174156 1.831215
## 5    5 1.749387 1.488693
## 6    6 1.344808 1.221217
# load observed data
obs <-read_delim(paste0('obs.txt'), delim = ',')
simple_therm_depth = parameters[18]
start_date <- get_yaml_value(config_file, "time", "start")
stop_date <- get_yaml_value(config_file, "time", "stop")

obs_sfc <- obs %>%
  filter(Depth_meter <= simple_therm_depth) %>%
  mutate('date' = yday(datetime),
         'sfc' = Water_Temperature_celsius) %>%
  group_by(datetime) %>%
  summarise('wtr_avg' = mean(sfc, na.rm = TRUE))
obs_sfc$time = match(as.Date(obs_sfc$datetime), seq(as.Date(start_date), as.Date(stop_date), by = 'day'))
obs_btm <- obs %>%
  filter(Depth_meter >= simple_therm_depth) %>%
  mutate('date' = yday(datetime),
         'btm' = Water_Temperature_celsius) %>%
  group_by(datetime) %>%
  summarise('wtr_avg' = mean(btm, na.rm = TRUE))
obs_btm$time = match(as.Date(obs_btm$datetime), seq(as.Date(start_date), as.Date(stop_date), by = 'day'))

ggplot(result) +
  geom_line(aes(x=Time, y=WT_epi, col='Surface Mixed Layer (model)')) +
  geom_line(aes(x=(Time), y=WT_hyp, col='Bottom Layer (model)')) +
  geom_point(data = obs_sfc, aes(x=time, y=wtr_avg, col='Surface Mixed Layer (obs)'),linetype = "dashed") + # sfc
  geom_point(data = obs_btm, aes(x=(time), y=wtr_avg, col='Bottom Layer (obs)'), linetype = "dashed") + # btm
  labs(x = 'Simulated Time', y = 'WT in deg C')  +
  scale_color_manual(values = c('blue','blue','red','red')) +
  theme_bw()+
  guides(col=guide_legend(title="Layer")) +
  theme(legend.position="bottom")

Dissolved oxygen simulation

We can now also add dissolved oxygen equations to the mix, where net ecosystem production and sediment oxygen demand are temperature-dependet but also fitting parameters (here using \(P_{NEP}\) and \(P_{SED}\)): \(\frac{dDO_{epi}}{dt}= F_{NEP} + F_{ATM} + {\upsilon}_{t,DO}*A_T *(\frac{DO_{hypo}}{V_{hypo}}-\frac{DO_{epi}}{V_{epi}})\) with \(F_{ATM} = k(DO_{sat}-\frac{DO_{epi}}{V_{epi}})* A_{surf}\) (gas exchange after the method according to Cole (1998)), and \(F_{NEP}= \theta^{T_{epi}-20} * P_{NEP} *V_{epi}\); and \(\frac{dDO_{hypo}}{dt}= {\upsilon}_{t,DO}*A_T *(\frac{DO_{epi}}{V_{epi}}-\frac{DO_{hypo}}{V_{hypo}})\) with \(F_{SED}= P_{SED} * A_{SED}* \theta^{T_{hypo}-20}(\frac{\frac{DO_{hypo}}{V_{hypo}}}{k_{1/2}+\frac{DO_{hypo}}{V_{hypo}}})\). The parameter \({\upsilon}_{t,DO}\) depends on the estimated vertical diffusivity but is set several magnitudes lower. Note that thermod itself tracks the mass of dissolved oxygen per layer.

# we need to configure additional parameters for oxygen dynamics first
wq_parameters <- append(parameters, c(0.001 / 1000, 
                                       100, 15000 * 1e4, 100))

wq_parameters[19] = parameters[19]

# simulation maximum length
times <- seq(from = 1, to = max(boundary$Day), by = 1)

# initial water temperatures and oxygen conditions
yini <- c(3,3, 10 * 1000/1e6  * wq_parameters[1], 10 * 1000/1e6  * wq_parameters[2]) 

ice_on = TRUE 
out <- run_oxygen_model(bc = boundary, params = wq_parameters, ini = yini, times = times, ice = ice_on)

result <- data.frame('Time' = out[,1],
                     'WT_epi' = out[,2], 'WT_hyp' = out[,3],
                     'DO_epi' = out[,4], 'DO_hyp' = out[,5])
head(result)
##   Time   WT_epi   WT_hyp       DO_epi       DO_hyp
## 1    1 3.000000 3.000000 2.743250e+12 2.125000e+12
## 2    2 2.903372 2.412467 2.789142e+12 2.116976e+12
## 3    3 2.499375 2.234035 2.832751e+12 2.109784e+12
## 4    4 2.174156 1.831215 2.872690e+12 2.106131e+12
## 5    5 1.749387 1.488693 2.912628e+12 2.102348e+12
## 6    6 1.344808 1.221217 2.952660e+12 2.097875e+12
# load observed oxygen data
obs <-read_delim(paste0('obs_oxygen.txt'), delim = ',')
simple_therm_depth = parameters[18]
start_date <- get_yaml_value(config_file, "time", "start")
stop_date <- get_yaml_value(config_file, "time", "stop")

obs_sfc <- obs %>%
  filter(Depth_meter <= simple_therm_depth) %>%
  mutate('date' = yday(datetime),
         'sfc' = Dissolved_Oxygen_gPerCubicMeter ) %>%
  group_by(datetime) %>%
  summarise('do_avg' = mean(sfc, na.rm = TRUE))
obs_sfc$time = match(as.Date(obs_sfc$datetime), seq(as.Date(start_date), as.Date(stop_date), by = 'day'))
obs_btm <- obs %>%
  filter(Depth_meter >= simple_therm_depth) %>%
  mutate('date' = yday(datetime),
         'btm' = Dissolved_Oxygen_gPerCubicMeter ) %>%
  group_by(datetime) %>%
  summarise('do_avg' = mean(btm, na.rm = TRUE))
obs_btm$time = match(as.Date(obs_btm$datetime), seq(as.Date(start_date), as.Date(stop_date), by = 'day'))

ggplot(result) +
  geom_line(aes(x=Time, y=DO_epi / 1000 /  wq_parameters[1] * 1e6, col='Surface Mixed Layer (model)')) +
  geom_line(aes(x=(Time), y=DO_hyp / 1000 /  wq_parameters[2] * 1e6, col='Bottom Layer (model)')) +
  geom_point(data = obs_sfc, aes(x=time, y=do_avg, col='Surface Mixed Layer (obs)'),linetype = "dashed") +
  geom_point(data = obs_btm, aes(x=(time), y=do_avg, col='Bottom Layer (obs)'),linetype = "dashed") + 
  labs(x = 'Simulated Time', y = 'DO in g/m3')  +
  scale_color_manual(values = c('blue','blue','red','red')) +
  theme_bw()+
  guides(col=guide_legend(title="Layer")) +
  theme(legend.position="bottom")