Integrating with the gasfluxes package

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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.2448754

There’s a fairly large jump from the R² 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)
#> Warning: package 'ggplot2' was built under R version 4.3.3
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)
#> Warning: package 'gasfluxes' was built under R version 4.3.3

# 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 linear.AICc
#> 1: test 0.06376292 0.0009964887 3.758958e-176  449.7366   1378.524    1378.605
#>    linear.RSE  linear.r linear.diagnostics robust.linear.f0 robust.linear.f0.se
#> 1:    2.39156 0.9654821                          0.06306838         0.001022341
#>    robust.linear.f0.p robust.linear.C0
#> 1:      4.611554e-170         449.9836
#>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      robust.linear.weights
#> 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
#> 1:                           0.157757 0.002322147 5.338986e-183 0.006730131
#>    HMR.phi  HMR.AIC HMR.AICc  HMR.RSE HMR.diagnostics   NDFE.f0 NDFE.f0.se
#> 1: 481.152 646.9649 647.1005 0.705413                 0.8551356  0.2842008
#>      NDFE.f0.p NDFE.tau  NDFE.C0 NDFE.AIC NDFE.AICc NDFE.RSE NDFE.diagnostics
#> 1: 0.002846487 2.131127 441.6831 906.8762  907.0118 1.087861

gasfluxes 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
#> 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
#> 1:    1378.605    2.39156 0.9654821                           0.1009094
#>    robust.linear.f0.se robust.linear.f0.p robust.linear.C0
#> 1:         0.001635746      4.611554e-170         449.9836
#>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      robust.linear.weights
#> 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
#> 1:                           0.2524112 0.003715436 5.338985e-183 0.006730131
#>    HMR.phi  HMR.AIC HMR.AICc  HMR.RSE HMR.diagnostics  NDFE.f0 NDFE.f0.se
#> 1: 481.152 646.9649 647.1005 0.705413                 1.368218  0.4547217
#>    NDFE.f0.p NDFE.tau  NDFE.C0 NDFE.AIC NDFE.AICc NDFE.RSE NDFE.diagnostics
#> 1: 0.0028465 2.131125 441.6831 906.8762  907.0118 1.087861                 
#>    SECONDS_min SECONDS_max
#> 1:           0         299

Conclusion

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.