goAquaFlux

Diffusive and ebullitive GHG flux partitioning for aquatic floating chambers

Overview

goAquaFlux extends the goFlux R package (Rheault et al., 2024) to aquatic floating-chamber measurements, where methane reaches the atmosphere by two pathways with very different behaviour:

Diffusion — a smooth, near-linear accumulation in the chamber headspace.
Ebullition — bubbles, appearing as abrupt steps in concentration.

Standard chamber-flux software models a single smooth signal, so a naive fit either discards bubble-affected incubations or silently folds ebullition into a biased “diffusive” estimate. goAquaFlux detects the steps, fits the diffusive flux on a bubble-free window using goFlux’s validated LM/HM models and QA/QC, and returns a partitioned total.

How it works

The workflow is a pipeline of four steps, one function each, coordinated by goAquaFlux() across all incubations in a dataset (split by UniqueID):

Detect the bubbles. find.bubbles() scans the concentration series (typically CH4dry_ppb) with a rolling dispersion metric — MAD, variance, or the MAD of first differences ("diff", robust to a diffusive trend) — and flags contiguous supra-threshold periods as bubbling events. Each event’s size is then estimated as a local step on top of the diffusive trend.
Fit the diffusive flux. goAquaFlux.diffusive() fits the diffusive component on the bubble-free window (truncated at the first event) using goFlux’s linear (LM) and non-linear (HM) models and the model-selection criteria. For non-bubble gases, the window is only cut at a bubble that coincides with an abrupt slope change.
Compute the ebullitive flux. goAquaFlux.ebullition() converts the summed bubble magnitudes into a flux, \(F_E = \frac{\sum \Delta C_{\text{bubble}}}{t_{inc}} \times K\), with standard errors propagated from the per-event magnitudes, and returns an independent two-point endpoint estimate for diagnosis.
Combine and check. goAquaFlux.total() adds the two components, \(F_T = F_E + F_D\), propagates the error assuming independence (\(SE_T = \sqrt{SE_E^2 + SE_D^2}\)), and flags results that are suspiciously larger than the endpoint estimate (e.g. a diffusive ramp mistaken for a bubble step).

flux.plot.aqua() then produces one diagnostic plot per incubation: the concentration series, the shaded diffusive window, the detected bubbling events, and a legend of the total / diffusive / ebullitive fluxes.

Installation

# install.packages("remotes")
remotes::install_github("Qepanna/goFlux")          # the parent package
remotes::install_github("Qepanna/goFlux@aquaghg")    # this fork + aqua extension

Quick start

library(goFlux)

fluxres <- goAquaFlux(
  dataframe            = IDed,          # output of goFlux::autoID()
  gastype              = "CH4dry_ppb",
  use_bubble_detection = TRUE
)

flux.plot.aqua(fluxres, IDed, gastype = "CH4dry_ppb",
               plot.display = c("diffusive.window", "ebullition.events"))

Functions

`goAquaFlux` for computing gas fluxes from aquatic chamber measurements

Main entry point: partitions aquatic GHG fluxes into diffusive and ebullitive components across many incubations.

Usage

Note

Code chunks under Usage sections are not part of the demonstration. They are meant to show you how to use the arguments in the function.

goAquaFlux(
  dataframe,
  gastype,
  H2O_col = "H2O_ppm",
  prec = NULL,
  criteria = c("MAE", "RMSE", "AICc", "SE", "g.factor", "kappa", "MDF", "nb.obs",
    "intercept", "p-value"),
  Area = NULL,
  offset = NULL,
  Vtot = NULL,
  Vcham = NULL,
  Pcham = NULL,
  Tcham = NULL,
  use_bubble_detection = TRUE,
  bubble.window.size = 30,
  bubble_gas = "CH4dry_ppb",
  bubble.method = "mad",
  bubble.args = list(),
  ebullition.final_window_min = 30,
  ebullition.window_C0Cf = 10,
  diffusion.minimum_window = 30,
  return_df = TRUE
)

Arguments

Argument	Description
`dataframe`	A `data.frame` containing time series observations for one or more chamber incubations.
`gastype`	Character string indicating the gas concentration column to use for flux calculation. Allowed values include: `"CO2dry_ppm"`, `"CH4dry_ppb"`, `"COdry_ppb"`, `"N2Odry_ppb"`, `"NH3dry_ppb"`, `"NO2dry_ppb"`, `"NOdry_ppb"`, and `"H2O_ppm"`.
`H2O_col`	Character string specifying the column containing water vapor concentration used for dilution correction. Default is `"H2O_ppm"`. If `NULL`, water vapor correction is disabled.
`prec`	Numeric scalar specifying instrument precision for the gas analyzer. If `NULL`, precision is retrieved from the corresponding precision column in `dataframe` (e.g. `CH4_prec`).
`criteria`	Character vector specifying model selection criteria used by `goFlux`. Default criteria include `"MAE"`, `"RMSE"`, `"AICc"`, `"SE"`, `"g.factor"`, `"kappa"`, `"MDF"`, `"nb.obs"`, `"intercept"`, and `"p-value"`.
`Area`	Numeric scalar; chamber base area in cm² (as in `goFlux`; note the internal factor of 10,000 to m²). If `NULL`, the value is retrieved from the `Area` column in `dataframe`.
`offset`	Numeric scalar; height between the water surface and the chamber top in cm, used to compute `Vtot` as `Vcham + Area * offset / 1000` when `Vtot` is not provided.
`Vtot`	Numeric scalar specifying total chamber volume (L). If `NULL`, it is calculated as: Vtot = Vcham + (Area * offset) / 1000
`Vcham`	Numeric scalar representing chamber headspace volume (L). Used only when `Vtot` is not provided.
`Pcham`	Numeric scalar representing chamber pressure (kPa). If `NULL`, atmospheric pressure (101.325 kPa) is assumed.
`Tcham`	Numeric scalar representing chamber temperature (°C). If `NULL`, a default temperature of 15 °C is assumed.
`use_bubble_detection`	Logical; if `TRUE` (default) bubbling events are detected using `bubble_gas` and, when `gastype == bubble_gas`, the flux is partitioned into diffusive and ebullitive components. If `FALSE`, no bubble detection is performed and only a diffusive flux is returned.
`bubble.window.size`	Integer specifying the rolling window size (number of observations) used for bubble detection.
`bubble_gas`	Character string specifying the gas used to detect bubbling events. Default is `"CH4dry_ppb"`.
`bubble.method`	Character; dispersion metric passed to `find.bubbles`: `"mad"` (default), `"variance"` or `"diff"` (rolling variance of increments, robust to a diffusive trend).
`bubble.args`	Named list of additional arguments forwarded to `find.bubbles` (e.g. `list(k = 5, min_magnitude = 10)`). Values here override the defaults for advanced tuning.
`ebullition.final_window_min`	Minimum time window (seconds) required after the last bubble event to define the incubation end for ebullition flux calculation.
`ebullition.window_C0Cf`	Time window (seconds) used to compute initial and final concentrations for the endpoint flux estimate.
`diffusion.minimum_window`	Minimum number of observations required to compute diffusive flux before the first bubble event.
`return_df`	Logical. If `TRUE` (default) the function returns a tidy list of three data frames (see Value). If `FALSE`, the raw per-incubation results list is returned instead, which is convenient for advanced users who want the untidied intermediate objects.

Details

The workflow includes: (1) data validation and cleaning; (2) calculation of auxiliary variables (flux conversion term and MDF); (3) bubble detection (typically using CH4 concentration); (4) estimation of ebullition flux based on detected bubble magnitudes; (5) estimation of diffusive flux using model selection; and (6) combination of the diffusive and ebullitive components into a total flux. The function operates on datasets containing multiple chamber incubations, which are automatically split using the UniqueID column. Bubble detection is performed using a rolling window approach applied to concentration variability of the selected gas (typically CH4). Detected bubbling events are used to estimate ebullition fluxes and restrict the time window used for diffusive flux estimation. Diffusive fluxes are calculated using the goFlux framework, which evaluates multiple regression models and selects the best model according to user-defined criteria.

Examples

# Not run:

results <- goAquaFlux(
  dataframe = chamber_data,
  gastype = "CH4dry_ppb"
)

References

Rheault, K., Christiansen, J. R., & Larsen, K. S. (2024). goFlux: A user-friendly way to calculate GHG fluxes yourself, regardless of user experience. Journal of Open Source Software, 9(96), 6393. } \seealso{ \code{

`find.bubbles` for detecting ebullition events

Detects bubbling (ebullition) events as abrupt steps in a concentration series.

Usage

Note

Code chunks under Usage sections are not part of the demonstration. They are meant to show you how to use the arguments in the function.

find.bubbles(
  df,
  bubble_source,
  window.size = 15,
  dt = 1,
  method = c("mad", "variance", "diff"),
  var.quantile = 0.7,
  k = 4,
  min_ratio = 3,
  min_sd = NULL,
  min_gap = 10,
  min_length = 5,
  max_reg_window = 120,
  reg.min.obs = 10,
  min_magnitude = 5,
  min_snr = NULL
)

Arguments

Argument	Description
`df`	A data.frame containing the incubation time series. Must include an `Etime` column (elapsed time, seconds) and the concentration column named by `bubble_source`.
`bubble_source`	Character; name of the concentration column used for detection (e.g. `"CH4dry_ppb"`).
`window.size`	Integer; width of the moving window (in interpolated time steps) used to compute rolling dispersion. Default `15`.
`dt`	Numeric; temporal resolution (seconds) of the regular grid the signal is interpolated onto before rolling statistics. Default `1`.
`method`	Character; dispersion metric. One of `"mad"` (default), `"variance"` or `"diff"`. `"diff"` computes rolling variance of the first differences and is largely insensitive to a linear diffusive trend, so it is recommended when diffusion is strong (see Details).
`var.quantile`	Numeric in (0, 1); empirical quantile of the rolling dispersion distribution used in the adaptive threshold. Default `0.7`.
`k`	Numeric; multiplier on the MAD of the rolling dispersion in the robust threshold (median + k MAD). Larger values are more conservative. Default `4`.
`min_ratio`	Numeric; minimum ratio of the maximum to the median rolling dispersion required before any detection is attempted. Guards against detections in flat, low-variance incubations. Default `3`.
`min_sd`	Numeric or `NULL`; if set, incubations whose overall standard deviation is below this value return `NULL` (no detection).
`min_gap`	Numeric; events separated by less than this gap (seconds) are merged. Default `10`.
`min_length`	Numeric; minimum event duration (seconds) to retain. Default `5`.
`max_reg_window`	Numeric; maximum half-window (seconds) used for the local magnitude regression around each event. Default `120`.
`reg.min.obs`	Integer; minimum observations required for the magnitude regression. Default `10`.
`min_magnitude`	Numeric or `NULL`; minimum absolute step magnitude (same units as `bubble_source`) to retain an event. Default `5`.
`min_snr`	Numeric or `NULL`; if set, minimum signal-to-noise ratio (`\|magnitude\| / SE`) to retain an event.

Details

The signal is first robustly standardised (median / MAD) and interpolated onto a regular dt-second grid. Rolling dispersion is then computed and compared to an adaptive threshold defined as the maximum of an empirical quantile (var.quantile) and a robust bound (median + k MAD). Contiguous supra-threshold runs are merged (min_gap) and filtered by duration (min_length). A steep but purely diffusive rise inflates the rolling variance/MAD even in the absence of bubbles, which can produce false positives. The "diff" method mitigates this by operating on the increments of the signal: a genuine ebullition step produces a large positive spike in the first differences, whereas a linear diffusive trend produces roughly constant increments and therefore low differenced dispersion. Event magnitude is estimated with the local model C_t = _0 + _1 t + _2 I(t t_b) where I(t t_b) is a step dummy at the event start; _2 is the estimated magnitude. Only events with a positive magnitude are retained (ebullition adds gas to the headspace).

Examples

# Not run:

bubbles <- find.bubbles(df = incubation_data,
                        bubble_source = "CH4dry_ppb",
                        window.size = 15,
                        method = "diff")

`goAquaFlux.diffusive` for estimating diffusive gas flux

Diffusive flux on a bubble-free window, via goFlux’s LM/HM models.

Usage

Note

Code chunks under Usage sections are not part of the demonstration. They are meant to show you how to use the arguments in the function.

goAquaFlux.diffusive(
  df,
  gastype,
  criteria = c("MAE", "RMSE", "AICc", "SE", "g.factor", "kappa", "MDF", "nb.obs",
    "intercept", "p-value"),
  bubble_gas = "CH4dry_ppb",
  bubbles = NULL,
  minimum_window = 30,
  abrupt.window = 30,
  abrupt.min.points = 10,
  abrupt.threshold = 0.5
)

Arguments

Argument	Description
`df`	A data.frame for a single incubation, as prepared by `goAquaFlux`. Must contain `Etime` (starting at 0), `flag`, the gas column named by `gastype`, the instrument precision column (`*_prec`), `Vtot`, `Area`, `Pcham`, `Tcham`, and (for non-water gases) `H2O_mol` and/or `H2O_ppm`.
`gastype`	Character; name of the gas concentration column.
`criteria`	Character vector of model-selection criteria passed to `best.flux`. Defaults to the full goFlux criteria set.
`bubble_gas`	Character; the gas used to detect bubbles (typically `"CH4dry_ppb"`). When `gastype == bubble_gas` the diffusive window is always truncated at the first detected bubble. For other gases the window is only truncated if a bubble coincides with an abrupt change in that gas (see Details).
`bubbles`	Data.frame of bubbling events from `find.bubbles` (or `NULL`). Must contain a `start` column.
`minimum_window`	Integer; minimum number of observations required in the diffusive window. If fewer are available the function returns `NA`.
`abrupt.window`	Numeric; half-width (seconds) of the window used on each side of a candidate bubble time when testing for an abrupt slope change (non-bubble gases only). Default `30`.
`abrupt.min.points`	Integer; minimum points required on each side of the candidate time to run the abrupt-change test. Default `10`.
`abrupt.threshold`	Numeric; relative slope-change threshold above which a bubble is deemed to perturb the gas. Default `0.5` (i.e. 50%).

Details

For the bubble gas itself, ebullition by definition perturbs the signal, so the diffusive window ends at the first detected bubble. For other gases a bubble does not necessarily perturb the concentration; the window is only cut at the first bubble time that is accompanied by an abrupt change in local slope (the ratio of the post- to pre-bubble slope exceeds abrupt.threshold). If no such change is found, the full series is used.

`goAquaFlux.ebullition` for estimating ebullitive gas flux

Ebullitive flux from detected bubble magnitudes.

Usage

Note

Code chunks under Usage sections are not part of the demonstration. They are meant to show you how to use the arguments in the function.

goAquaFlux.ebullition(
  df,
  gastype,
  bubbles,
  flux.term,
  final_window.min = 30,
  window_C0Cf = 10
)

Arguments

Argument	Description
`df`	A data.frame for a single incubation. Must contain an `Etime` column (elapsed time, seconds, starting at 0) and the gas concentration column named by `gastype`. Rows are assumed to be sorted by `Etime`.
`gastype`	Character; name of the gas concentration column (e.g. `"CH4dry_ppb"`).
`bubbles`	Data.frame of bubbling events from `find.bubbles` (or `NULL`). Must contain `start`, `end` and `magnitude`; an optional `SE` column enables error propagation.
`flux.term`	Numeric; conversion factor turning a concentration change per unit time into flux units (from `flux.term`).
`final_window.min`	Numeric; minimum time (seconds) that must remain after the last bubble for the full deployment length to be used as the incubation time. If less time remains, the incubation time is truncated to the start of the last bubbling event. Default `30`.
`window_C0Cf`	Numeric; duration (seconds) of the initial and final windows used to compute mean start and end concentrations for the two-point endpoint estimate. Default `10`.

Details

Ebullitive flux: F_E = K with K the flux conversion term and t_{inc} the effective incubation time. Standard errors are propagated assuming independent bubble magnitude estimates: SE(F_E) =

`goAquaFlux.total` for combining diffusive and ebullitive flux

Combine components with error propagation and consistency checks.

Usage

Note

Code chunks under Usage sections are not part of the demonstration. They are meant to show you how to use the arguments in the function.

goAquaFlux.total(ebullition_flux, diffusive_flux, tolerance = 1.2)

Arguments

Argument	Description
`ebullition_flux`	A list returned by `goAquaFlux.ebullition`. Must contain at least `flux` (ebullitive flux), `SE` (its standard error) and `F_tot2pts` (the endpoint total-flux estimate).
`diffusive_flux`	A list returned by `goAquaFlux.diffusive`. Must contain at least `flux` (diffusive flux) and `SE`.
`tolerance`	Numeric > 1. Ratio above which the reconstructed total flux is flagged as suspiciously larger than the endpoint estimate. Default `1.2` (i.e. 20% larger).

Details

Total flux: F_T = F_E + F_D Error propagation (independent errors): SE_T =

`flux.plot.aqua` for plotting aquatic chamber incubations

Per-incubation diagnostic plots of the diffusive window and ebullition events.

Usage

Note

Code chunks under Usage sections are not part of the demonstration. They are meant to show you how to use the arguments in the function.

flux.plot.aqua(
  flux.results.ls,
  dataframe,
  gastype,
  shoulder = 30,
  plot.display = c("diffusive.window", "ebullition.events"),
  flux.unit = NULL,
  quality.check = FALSE,
  conversion.factor = 1
)

Arguments

Argument	Description
`flux.results.ls`	The list returned by `goAquaFlux` (with `return_df = TRUE`), i.e. a list with `flux_summary`, `bubbles` and `diffusive`. For backwards compatibility, a plain `best.flux`-style data.frame may also be supplied, in which case the function delegates to `flux.plot` (diffusion only).
`dataframe`	The data.frame of measurements used by `goAquaFlux` (with `Etime`, `flag` and the gas column).
`gastype`	Character; gas column to plot (e.g. `"CH4dry_ppb"`).
`shoulder`	Numeric; padding (seconds) added before/after the measurement on the x-axis. Default `30`.
`plot.display`	Character vector of elements to overlay. Currently supported: `"diffusive.window"` and `"ebullition.events"`. Default `c("diffusive.window", "ebullition.events")`.
`flux.unit`	Character or `NULL`; flux-unit label. If `NULL`, a sensible default is chosen from `gastype`.
`quality.check`	Logical; reserved for quality-check annotations. Default `FALSE`.
`conversion.factor`	Numeric > 0; multiplies displayed flux values. Default `1`.

Overview

How it works

Installation

Quick start

Functions

goAquaFlux for computing gas fluxes from aquatic chamber measurements

Usage

Arguments

Details

Examples

References

find.bubbles for detecting ebullition events

Usage

Arguments

Details

Examples

goAquaFlux.diffusive for estimating diffusive gas flux

Usage

Arguments

Details

goAquaFlux.ebullition for estimating ebullitive gas flux

Usage

Arguments

Details

goAquaFlux.total for combining diffusive and ebullitive flux

Usage

Arguments

Details

flux.plot.aqua for plotting aquatic chamber incubations

Usage

Arguments

`goAquaFlux` for computing gas fluxes from aquatic chamber measurements

`find.bubbles` for detecting ebullition events

`goAquaFlux.diffusive` for estimating diffusive gas flux

`goAquaFlux.ebullition` for estimating ebullitive gas flux

`goAquaFlux.total` for combining diffusive and ebullitive flux

`flux.plot.aqua` for plotting aquatic chamber incubations