Forward model comparison
This document provides code and analysis to compare metabolism
results from the Odum open-water method, BASEmetab, and EBASE with known
results from the forward model included in the Fwoxy package. The Odum
open-water method is available from the sine branch of
the WtRegDO package, which includes a new option to output instantaneous
metabolism flux estimates at each time step. The BASEmetab results were
created using this
script. For BASEmetab, the prior distribution for the gas transfer
velocity parameter (K) was based on a mean value from Fwoxy with an
appropriate standard deviation in proportion to the value used in EBASE.
Also note that the BASEmetab model requires a PAR time series, which was
created using the fun_par_sin_model from Fwoxy. The EBASE
results were created using this
script and is similar to BASEmetab with the exception of a
parameterization for gas exchange following Wanninkhof 2014, i.e., gas
exchange is estimated directly in the model and not calculated
empirically outside of the Bayesian context.
Setup and analysis
Load the libraries, BASEmetab results, EBASE results, and setup model parameters and constants for the simulated dissolved oxygen time series from Fwoxy.
library(fwoxy)
library(tidyverse)
library(lubridate)
library(WtRegDO)
library(here)
load(file = here('data/fwoxybasemetab.RData'))
load(file = here('data/fwoxyebase.RData'))
# plot colors
colors <- c(gasexd = "red3", gppd = "orange", erd = "purple4", dcdtd = "steelblue3")
# plot breaks, labels
brks <- seq(1, 518400,by = 43200)
labs <- rep(c(0, 12), length.out = length(brks))
# Set model parameters
oxy_ic <- 250 # (mmol/m^3), initial oxygen concentration
a_param <- 0.2 # ((mmol/m^3)/day)/(W/m^2), light efficiency
er_param <- 20 # (mmol/m^3/day), ecosystem respiration
# Constant Forcings
ht_const <- 3 # m, height of the water column
salt_const <- 25 # ppt, salinity
temp_const <- 25 # deg C, water temperature
wspd_const <- 3 # m/s, wind speed at 10 mCreate Fwoxy model using the parameters above. The
example object creates a plot showing the oxygen time
series and resulting metabolism estimates as flux of dissolved oxygen in
mmol/m3/d. These are volumetric rates.
example <- fwoxy(oxy_ic = oxy_ic, a_param = a_param, er_param = er_param,
ht_in = ht_const, salt_in = salt_const, temp_in = temp_const,
wspd_in = wspd_const)
head(example)## time, sec oxy, mmol/m3 troc, mmol/m3/d gasex, mmol/m3/d gpp, mmol/m3/d
## 1 0 250.0000 NA NA NA
## 2 900 249.7305 -25.86951 5.869511 0
## 3 1800 249.4617 -25.80931 5.809309 0
## 4 2700 249.1935 -25.74925 5.749247 0
## 5 3600 248.9259 -25.68932 5.689325 0
## 6 4500 248.6589 -25.62954 5.629542 0
## er, mmol/m3/d oxysu, mmol/m3 wspd2, m2/s2 sc kw, m/s
## 1 NA NA NA NA NA
## 2 20 26.27292 9 431.8863 7.757129e-06
## 3 20 26.00345 9 431.8863 7.757129e-06
## 4 20 25.73460 9 431.8863 7.757129e-06
## 5 20 25.46638 9 431.8863 7.757129e-06
## 6 20 25.19878 9 431.8863 7.757129e-06
## oxysat, mmol/m3
## 1 NA
## 2 223.7271
## 3 223.7271
## 4 223.7271
## 5 223.7271
## 6 223.7271The oxygen time series can be used from the Fwoxy output with the Odum open-water method in WtRegDO. The time series must be formatted to be used with WtRegDO, which includes converting the seconds column to a date/time object (arbitrarily set to the system date), changing the units of oxygen from mmol/m3 to mg/L, and including the constants defined for the Fwoxy input as input to WtRegDO.
tomod <- example %>%
mutate(
DateTimeStamp = force_tz(as.POSIXct(`time, sec`, origin = Sys.Date(), tz = 'UTC'), tzone = 'America/Jamaica'),
DO_obs = `oxy, mmol/m3` * 0.032, # to mg/L, same as g/m3
Temp = temp_const,
WSpd = wspd_const,
Sal = salt_const,
ATemp = NA,
BP = 1013.25,
Tide = NA
) %>%
select(DateTimeStamp, Temp, Sal, DO_obs, WSpd, ATemp, BP, Tide)
head(tomod)## DateTimeStamp Temp Sal DO_obs WSpd ATemp BP Tide
## 1 2023-03-15 00:00:00 25 25 8.000000 3 NA 1013.25 NA
## 2 2023-03-15 00:15:00 25 25 7.991377 3 NA 1013.25 NA
## 3 2023-03-15 00:30:00 25 25 7.982774 3 NA 1013.25 NA
## 4 2023-03-15 00:45:00 25 25 7.974191 3 NA 1013.25 NA
## 5 2023-03-15 01:00:00 25 25 7.965628 3 NA 1013.25 NA
## 6 2023-03-15 01:15:00 25 25 7.957084 3 NA 1013.25 NAThe Odum open-water method is then used. Note that the location and time zone is required. This will impose some slight differences in the metabolic estimates between Fwoxy and the Odum method because the latter uses a “metabolic day” as the time between sunrises based on the date and location. Also note the chosen method for estimating gas exchange. The results return areal (m2) and volumetric (m3) estimates for metabolism. The first several rows do not include estimates because of internal checks in WtRegDO that verify a “complete” metabolic day is used for the calculations. For the Fwoxy oxygen output, the WtRegDO package treats the first few hours of data as the tail end of an incomplete metabolic day and results are not estimated.
# estimate instantaneous Odum metabolism from Fwoxy output
opmetab <- ecometab(tomod, DO_var = 'DO_obs', tz = 'America/Jamaica', lat = 29.75, long = -85, depth_val = NULL,
depth_vec = ht_const, gasex = 'Wanninkhof', gasave = 'instant', metab_units = 'mmol', instant = T)
head(opmetab)## DateTimeStamp Temp Sal WSpd BP DO metab_date solar_period
## 1 2023-03-15 00:07:30 25 25 3 1013.25 249.8653 2023-03-14 sunset
## 2 2023-03-15 00:22:30 25 25 3 1013.25 249.5961 2023-03-14 sunset
## 3 2023-03-15 00:37:30 25 25 3 1013.25 249.3276 2023-03-14 sunset
## 4 2023-03-15 00:52:30 25 25 3 1013.25 249.0597 2023-03-14 sunset
## 5 2023-03-15 01:07:30 25 25 3 1013.25 248.7924 2023-03-14 sunset
## 6 2023-03-15 01:22:30 25 25 3 1013.25 248.5257 2023-03-14 sunset
## solar_time day_hrs DOsat dDO H D KL
## 1 2023-03-14 19:47:36 11.94053 1.115411 -25.86951 3 -5.775812 0.6702159
## 2 2023-03-14 19:47:36 11.94053 1.114210 -25.80931 3 -5.715680 0.6702159
## 3 2023-03-14 19:47:36 11.94053 1.113011 -25.74925 3 -5.655688 0.6702159
## 4 2023-03-14 19:47:36 11.94053 1.111815 -25.68932 3 -5.595836 0.6702159
## 5 2023-03-14 19:47:36 11.94053 1.110622 -25.62954 3 -5.536122 0.6702159
## 6 2023-03-14 19:47:36 11.94053 1.109432 -25.56990 3 -5.476548 0.6702159
## Ka Pg Rt NEM Pg_vol Rt_vol
## 1 0.009308554 NA NA NA NA NA
## 2 0.009308554 NA NA NA NA NA
## 3 0.009308554 NA NA NA NA NA
## 4 0.009308554 NA NA NA NA NA
## 5 0.009308554 NA NA NA NA NA
## 6 0.009308554 NA NA NA NA NAComparison
The metabolic output from WtRegDO, BASEmetab, EBASE, and Fwoxy are then combined for a direct comparison. This requires a bit of additional wrangling, the details of which are not important. The output shows the combined time series used to create the comparison plots.
# function to make seconds in opmetab match fwoxy
getsec <- function(x){
sec <- seconds(x)
sec <- sec - (7.5 * 60)
sec <- sec - min(sec)
sec <- as.numeric(sec)
return(sec)
}
# format opmetab to match fwoxy
opmetab <- opmetab %>%
select(DateTimeStamp, dDO, D, Pg_vol, Rt_vol) %>%
mutate(
secs = getsec(DateTimeStamp),
D = -1 * D,
Rt_vol = -1 * Rt_vol,
typ = 'Odum'
) %>%
select(-DateTimeStamp) %>%
gather(var, val, -secs, -typ)
# format BASEmetab to match fwoxy
bsmetab <- fwoxybasemetab %>%
mutate(
secs = getsec(DateTimeStamp),
typ = 'BASEmetab'
) %>%
select(secs, typ, Pg_vol, Rt_vol, D, dDO) %>%
gather(var, val, -secs, -typ) %>%
na.omit()
# format EBASE to match fwoxy
ebmetab <- fwoxyebase %>%
mutate(
secs = getsec(DateTimeStamp),
D = -1 * D / H,
Pg_vol = P / H,
Rt_vol = R / H,
typ = 'EBASE'
) %>%
select(secs, typ, Pg_vol, Rt_vol, D, dDO) %>%
gather(var, val, -secs, -typ) %>%
na.omit()
# combine fwoxy and opmetab
toplo1 <- example %>%
select(
secs = `time, sec`,
dDO = `troc, mmol/m3/d`,
D = `gasex, mmol/m3/d`,
Pg_vol = `gpp, mmol/m3/d`,
Rt_vol = `er, mmol/m3/d`
) %>%
mutate(
typ = 'Fwoxy'
) %>%
gather(var, val, -secs, -typ) %>%
bind_rows(opmetab, bsmetab, ebmetab) %>%
mutate(
var = factor(var, levels = c('D', 'Pg_vol', 'Rt_vol', 'dDO'), labels = c('gasexd', 'gppd', 'erd', 'dcdtd')),
typ = factor(typ, levels = c('Fwoxy', 'Odum', 'BASEmetab', 'EBASE'))
)
head(toplo1)## secs typ var val
## 1 0 Fwoxy dcdtd NA
## 2 900 Fwoxy dcdtd -25.86951
## 3 1800 Fwoxy dcdtd -25.80931
## 4 2700 Fwoxy dcdtd -25.74925
## 5 3600 Fwoxy dcdtd -25.68932
## 6 4500 Fwoxy dcdtd -25.62954This shows the time series for the metabolic estimates provided by each method. The plot format is similar to that provided by the Fwoxy package.
ggplot(toplo1, aes(x = secs, y = val, color = var)) +
geom_line() +
facet_wrap(~typ, ncol = 1) +
scale_color_manual(values = colors) +
scale_x_continuous(labels = labs, breaks = brks) +
theme_bw() +
labs(
y = 'Flux, mmol/m3/d',
x = 'Hour of day'
) +
theme(
legend.title = element_blank(),
strip.background = element_blank()
)
This shows a direct comparison of the estimates provided by each method. The y-axis shows the Odum results and the x-axis shows the Fwoxy results. There is a slight difference in phasing of the time series that creates an offset from the 1:1 line. Also note that the difference in respiration is minor but related to the daily estimate for the Odum method versus an overall estimate for the entire time series for the Fwoxy method.
toplo2 <- toplo1 %>%
spread(typ, val)
ggplot(toplo2, aes(x = Fwoxy, y = Odum, color = var)) +
geom_point() +
facet_wrap(~var, ncol = 2, scales = 'free') +
# coord_equal() +
scale_color_manual(values = colors) +
theme_bw() +
geom_abline(intercept = 0, slope = 1) +
labs(
y = 'Odum, Flux, mmol/m3/d',
x = 'Fwoxy, flux, mmol/m3/d'
) +
theme(
legend.title = element_blank(),
strip.background = element_blank()
)
The same comparison as above is made for the BASEmetab results with Fwoxy.
ggplot(toplo2, aes(x = Fwoxy, y = BASEmetab, color = var)) +
geom_point() +
facet_wrap(~var, ncol = 2, scales = 'free') +
scale_color_manual(values = colors) +
theme_bw() +
geom_abline(intercept = 0, slope = 1) +
labs(
y = 'BASEmetab, Flux, mmol/m3/d',
x = 'Fwoxy, Flux, mmol/m3/d'
) +
theme(
legend.title = element_blank(),
strip.background = element_blank()
)
The same comparison as above is made for the EBASE results with Fwoxy.
ggplot(toplo2, aes(x = Fwoxy, y = EBASE, color = var)) +
geom_point() +
facet_wrap(~var, ncol = 2, scales = 'free') +
scale_color_manual(values = colors) +
theme_bw() +
geom_abline(intercept = 0, slope = 1) +
labs(
y = 'EBASE, Flux, mmol/m3/d',
x = 'Fwoxy, Flux, mmol/m3/d'
) +
theme(
legend.title = element_blank(),
strip.background = element_blank()
)
This shows a direct comparison between the methods for each metabolic result as time series.
ggplot(toplo1, aes(x = secs, y = val, color = typ)) +
geom_line() +
facet_wrap(~var, ncol = 1, scales = 'free') +
scale_x_continuous(labels = labs, breaks = brks) +
theme_bw() +
labs(
y = 'Flux, mmol/m3/d',
x = 'Hour of day'
) +
theme(
legend.title = element_blank(),
strip.background = element_blank()
)
Summary
All methods produce similar results for the instantaneous metabolic
estimates, although EBASE over-estimates gas exchange and
under-estimates respiration. Fixing the b parameter in the
EBASE method produces comparable estimates for the other methods. Slight
differences between the Odum and Fwoxy results relate to how the time
series is divided to estimate the rates in the WtRegDO package.
Specifically, the Odum method uses the metabolic day, whereas the Fwoxy
approach uses a fixed day. Additionally, the Odum method calculates
rates based on the midpoint difference of DO between two time steps.
This creates a slight phase offset between the methods. The BASEmetab
results are more similar to the Fwoxy results than the Odum results.
For EBASE, we can compare the individual parameter estimates from
values used to create the Fwoxy input. The a parameter is
light efficiency and was set at 0.2 (mmol/m3/d)/(W/m2) in the initial
Fwoxy input. The estimate from EBASE is:
a <- fwoxyebase %>%
pull(a) %>%
na.omit() %>%
unique %>%
mean()
a## [1] 0.2013696This is very close to the input used for Fwoxy. A second parameter
output from Fwoxy is b as part of the gas exchange
estimate. A reasonable estimate for this value is 0.251 (cm/hr)/(m2/s2).
The estimate from EBASE is:
b <- fwoxyebase %>%
pull(b) %>%
na.omit() %>%
unique %>%
mean()
b## [1] 0.4961419