Integrating with the gasfluxes package
As noted in the introductory vignette, there are many reasons why greenhouse gas concentration observations might not be well fit by a linear model: diffusion effects, saturation, mixing problems, etc.
The fluxfinder package provides extensive outputs based
on linear model fits, and a few QA/QC indicators based on a polynomial
model, but there are cases where we might want to take advantage of more
sophisticated model-fitting routines.
Some problematic sample data
library(fluxfinder)
# Load the data
# Data from a LI-7810
f <- system.file("extdata/TG10-01087-curvature.data", package = "fluxfinder")
dat <- ffi_read_LI7810(f)
#> TG10-01087-curvature.data: read 300 rows of TG10-01087 data, 2022-07-12 10:07:00 to 2022-07-12 10:11:59 EST
# Fit a linear flux and QA/QC it
flux <- ffi_compute_fluxes(dat, group_column = NULL, time_column = "TIMESTAMP", gas_column = "CO2")
#> NOTE: HM81_flux.estimate is not NA, implying nonlinear data
print(flux$lin_flux.estimate)
#> [1] 0.1020207
print(flux$lin_r.squared)
#> [1] 0.9321557
print(flux$poly_r.squared)
#> [1] 0.994207
print(flux$HM81_flux.estimate)
#> [1] 0.2448754There’s a fairly large jump from the R2 of the linear
model (0.93) to that of a polynomial (0.99+), and the flux computed by
the Hutchinson
and Mosier (1981) nonlinear approach is numeric (i.e.,
non-NA).
This implies nonlinearity in our data:
library(ggplot2)
dat$SECONDS <- dat$SECONDS-min(dat$SECONDS)
ggplot(dat, aes(SECONDS, CO2)) + geom_point() +
# naive linear model
geom_smooth(method = "lm", formula = 'y ~ x') +
# HM1981 flux line
geom_abline(slope = flux$HM81_flux.estimate, intercept = min(dat$CO2), linetype = 2)
There’s definitely a problem! Depending on what we think is happening
here, one option would be to change the observation length so that the
flux is computed based on only the first ~75 seconds, which looks
linear. A second option would be to use the flux_HM1981
field as our flux.
A third option would be fit more sophisticated model(s).
Using the gasfluxes package
The gasfluxes package also provides routines to calculate greenhouse gas fluxes from chamber measurements, and includes code to fit the HMR model as well as several model-selection routines.
library(gasfluxes)
# Add some columns that gasfluxes expects
dat$ID <- "test"
dat$V <- 0.1
dat$A <- 0.16
gasfluxes(dat, .times = "SECONDS", .C = "CO2", plot = FALSE)
#> test: lm fit successful
#> test: rlm fit successful
#> test: HMR fit successful
#> test: NDFE fit successful
#> ID linear.f0 linear.f0.se linear.f0.p linear.C0 linear.AIC
#> <char> <num> <num> <num> <num> <num>
#> 1: test 0.06376292 0.0009964887 3.758958e-176 449.7366 1378.524
#> linear.AICc linear.RSE linear.r linear.diagnostics robust.linear.f0
#> <num> <num> <num> <char> <num>
#> 1: 1378.605 2.39156 0.9654821 0.06306838
#> robust.linear.f0.se robust.linear.f0.p robust.linear.C0
#> <num> <num> <num>
#> 1: 0.001022341 4.611554e-170 449.9836
#> robust.linear.weights
#> <char>
#> 1: 0.65|0.69|0.75|0.61|0.73|0.68|0.66|0.82|0.78|0.79|0.9|1|0.86|0.87|0.68|0.78|1|1|0.97|0.87|0.93|0.95|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|0.93|1|1|1|1|1|1|0.98|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|0.75|1|0.92|1|1|1|0.93|1|1|0.91|0.82|1|0.89
#> robust.linear.diagnostics HMR.f0 HMR.f0.se HMR.f0.p HMR.kappa
#> <char> <num> <num> <num> <num>
#> 1: 0.157757 0.002322147 5.338988e-183 0.006730131
#> HMR.phi HMR.AIC HMR.AICc HMR.RSE HMR.diagnostics NDFE.f0 NDFE.f0.se
#> <num> <num> <num> <num> <char> <num> <num>
#> 1: 481.152 646.9649 647.1005 0.705413 0.855136 0.2842011
#> NDFE.f0.p NDFE.tau NDFE.C0 NDFE.AIC NDFE.AICc NDFE.RSE NDFE.diagnostics
#> <num> <num> <num> <num> <num> <num> <char>
#> 1: 0.002846498 2.131125 441.6831 906.8762 907.0118 1.087861gasfluxes will compute on multiple groups (measurements)
via its .id parameter, but we can also use use
fluxfinder::ffi_compute_fluxes() to handle the grouping and
let it call gasfluxes:
# Define a small wrapper function to put things into the format
# that gasfluxes expects
f <- function(time, conc) {
d <- data.frame(time = time, C = conc, ID = 1, A = 1, V = 1)
gasfluxes(d, plot = FALSE)
}
ffi_compute_fluxes(dat, NULL, "SECONDS", "CO2", fit_function = f)
#> 1: lm fit successful
#> 1: rlm fit successful
#> 1: HMR fit successful
#> 1: NDFE fit successful
#> SECONDS ID linear.f0 linear.f0.se linear.f0.p linear.C0 linear.AIC
#> <num> <num> <num> <num> <num> <num> <num>
#> 1: 149.5 1 0.1020207 0.001594382 3.758958e-176 449.7366 1378.524
#> linear.AICc linear.RSE linear.r linear.diagnostics robust.linear.f0
#> <num> <num> <num> <char> <num>
#> 1: 1378.605 2.39156 0.9654821 0.1009094
#> robust.linear.f0.se robust.linear.f0.p robust.linear.C0
#> <num> <num> <num>
#> 1: 0.001635746 4.611554e-170 449.9836
#> robust.linear.weights
#> <char>
#> 1: 0.65|0.69|0.75|0.61|0.73|0.68|0.66|0.82|0.78|0.79|0.9|1|0.86|0.87|0.68|0.78|1|1|0.97|0.87|0.93|0.95|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|0.93|1|1|1|1|1|1|0.98|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|1|0.75|1|0.92|1|1|1|0.93|1|1|0.91|0.82|1|0.89
#> robust.linear.diagnostics HMR.f0 HMR.f0.se HMR.f0.p HMR.kappa
#> <char> <num> <num> <num> <num>
#> 1: 0.2524112 0.003715436 5.338988e-183 0.006730131
#> HMR.phi HMR.AIC HMR.AICc HMR.RSE HMR.diagnostics NDFE.f0 NDFE.f0.se
#> <num> <num> <num> <num> <char> <num> <num>
#> 1: 481.152 646.9649 647.1005 0.705413 1.368217 0.4547209
#> NDFE.f0.p NDFE.tau NDFE.C0 NDFE.AIC NDFE.AICc NDFE.RSE NDFE.diagnostics
#> <num> <num> <num> <num> <num> <num> <char>
#> 1: 0.002846472 2.131128 441.6831 906.8762 907.0118 1.087861
#> SECONDS_min SECONDS_max
#> <int> <int>
#> 1: 0 299Conclusion
This vignette covered how to integrate gasfluxes with
fluxfinder, if you want to take advantage of the former
package’s HMR or NDFE routines. More information on these models can be
found in:
- Hüppi, R., Felber, R., Krauss, M., Six, J., Leifeld, J., and Fuß, R.: Restricting the nonlinearity parameter in soil greenhouse gas flux calculation for more reliable flux estimates, PLoS One, 13, e0200876, 2018. http://dx.doi.org/10.1371/journal.pone.0200876
- Leiber-Sauheitl, K., Fuß, R., Voigt, C., and Freibauer, A.: High CO2 fluxes from grassland on histic Gleysol along soil carbon and drainage gradients, Biogeosciences, 11, 749–761, 2014. http://dx.doi.org/10.5194/bg-11-749-2014
- Pedersen, A. R., Petersen, S. O., and Schelde, K.: A comprehensive approach to soil-atmosphere trace-gas flux estimation with static chambers, Eur. J. Soil Sci., 61, 888–902, 2010. http://dx.doi.org/10.1111/j.1365-2389.2010.01291.x