How to calculate gasfluxes

Roman Hüppi, Roland Fuß

2018-08-04

This vignette introduces the gasfluxes package to calculate soil greenhouse gas fluxes from static chamber measurements. The package offers different models for measured concentration-time relationships from static chambers within the function gasfluxes

and offers the function selectfluxes to use selection algorithms that combine the different models appropriately

Data import

The input data.frame or data.table can be imported easily from a CSV file (e.g., as exported by Excel):

library(data.table)
#adjust path (see help("setwd")) and file name as appropriate
fluxMeas <- fread("fluxmeas.csv") 
#here we use two flux measurements from the file as an example
fluxMeas <- fluxMeas[ID %in% c("ID6","ID11")]

fluxMeas
#>        ID        V     A      time         C
#>    <char>    <num> <int>     <num>     <num>
#> 1:    ID6 0.535125     1 0.0000000 0.3241982
#> 2:    ID6 0.535125     1 0.3333333 0.3838311
#> 3:    ID6 0.535125     1 0.6666667 0.3623941
#> 4:    ID6 0.535125     1 1.0000000 0.5067303
#> 5:   ID11 0.550000     1 0.0000000 0.3337777
#> 6:   ID11 0.550000     1 0.3333333 0.4410218
#> 7:   ID11 0.550000     1 0.6666667 0.5080499
#> 8:   ID11 0.550000     1 1.0000000 0.5417371

The ID column must be a unique identifier for each flux measurement, hence the same value for all the concentration-time points that belong to the same flux. Column V stands for the chamber volume, A for the chamber area, time for the elapsed time after chamber closure and C is the concentration of the measured species (CO2, N2O, CH4 etc.)

Units of the output fluxes depend on input units. It’s recommended to use [V] = m3, [A] = m2, [time] = h, [C] = [mass or mol]/m3, which results in [f0] = [mass or mol]/m2/h. Since all algorithms only use V/A, A can be input as 1 and V as the chamber height.

The following features of the input will be checked for sanity:

Calculating fluxes with the gasfluxes function

If the input dataframe has the appropriate format, plug it into the function to calculate the fluxes:

library(gasfluxes)
flux.results <- gasfluxes(fluxMeas, method = c("linear","robust linear", "HMR"), plot = TRUE)
#> ID6: lm fit successful
#> ID6: rlm fit successful
#> ID6: HMR fit failed
#> ID11: lm fit successful
#> ID11: rlm fit successful
#> ID11: HMR fit successful
flux.results
#>        ID linear.f0 linear.f0.se linear.f0.p linear.C0 linear.AIC linear.AICc
#>    <char>     <num>        <num>       <num>     <num>      <num>       <num>
#> 1:    ID6 0.0844683   0.03531910  0.13923269 0.3153645  -9.516839         Inf
#> 2:   ID11 0.1139995   0.01920685  0.02723196 0.3525107 -14.609437         Inf
#>    linear.RSE  linear.r linear.diagnostics robust.linear.f0 robust.linear.f0.se
#>         <num>     <num>             <char>            <num>               <num>
#> 1: 0.04919468 0.8607673                          0.09719292         0.002805431
#> 2: 0.02602898 0.9727680                          0.11399952         0.019206850
#>    robust.linear.f0.p robust.linear.C0 robust.linear.weights
#>                 <num>            <num>                <char>
#> 1:       0.0006741157        0.3232908           1|1|0.033|1
#> 2:       0.0272319550        0.3525107               1|1|1|1
#>    robust.linear.diagnostics    HMR.f0  HMR.f0.se   HMR.f0.p HMR.kappa
#>                       <char>     <num>      <num>      <num>     <num>
#> 1:                                  NA         NA         NA        NA
#> 2:                           0.2332008 0.01194134 0.03257046  1.628278
#>      HMR.phi   HMR.AIC HMR.AICc     HMR.RSE
#>        <num>     <num>    <num>       <num>
#> 1:        NA        NA       NA          NA
#> 2: 0.5938413 -33.15831       NA 0.002821235
#>                                                    HMR.diagnostics
#>                                                             <char>
#> 1: Missing value or an infinity produced when evaluating the model
#> 2:

The output of gasfluxes contains the estimate of the initial flux (f0) for each method chosen (linear, robust linear and HMR in this case) including statistical parameters of the fits:

For the improved HMR fitting function the starting value k_HMR = log(1.5) is used by default. This is appropriate for hourly values. If you change the input i.e. by providing your measurement time in seconds use k_HMR = log(1.5/3600).

It is possible to have more than one ID column in the input file. For instance, it is often convenient to use a plot_ID and a date column by specifying the column names, e.g., gasfluxes(fluxMeas, .id = c("plot_ID", "date"), method = c("linear","robust linear", "HMR")).

###Graphical output

By default and with plot = TRUE the flux calcuation function plots a figure for each chamber flux that is stored in the subfolder /pics in the working directory. The following example on the left shows a flux where HMR could not be fitted and the 3rd concentration is regarded as outlier by the robust linear estimate. The second examples shows a flux that is well fitted by HMR which increases the flux estimate twofold compared to the linear fit (see numbers in the output above; HMR.f0/linear.f0).

Applying selection algorithm with the selectfluxes function

It is not apriori obvious which of the flux estimates should be selected for an individual chamber measurement. In case of small fluxes and large relativ measurement uncertainty linear estimates are the most appropriate choice. However, if non-linear effects (like decreasing diffusion gradient, lateral gas flow and chamber leakage) impact the increase in concentration within the chamber, the non-inear estimate of the HMR model is a better choice. The selectfluxes function offers some well defined decision algorithms that avoid the need for manual decisions (like suggested in the original HMR package). The most extensively tested decision ruleset is kappa.max, which optimises the balance between bias and uncertainty by taking the minimal detectable flux of the current system (f.detect see its suggested calculation below) and the size of the flux into account (see Hüppi et al. (2018) for details). t.meas is the time of the final concentration time point of the individual chamber measurement. To apply the kappa.max selection to your calculated fluxes use the selectfluxes function the following way:

selectfluxes(flux.results, select = "kappa.max", f.detect = 0.031, t.meas = 1)
#>        ID linear.f0 linear.f0.se linear.f0.p linear.C0 linear.AIC linear.AICc
#>    <char>     <num>        <num>       <num>     <num>      <num>       <num>
#> 1:    ID6 0.0844683   0.03531910  0.13923269 0.3153645  -9.516839         Inf
#> 2:   ID11 0.1139995   0.01920685  0.02723196 0.3525107 -14.609437         Inf
#>    linear.RSE  linear.r linear.diagnostics robust.linear.f0 robust.linear.f0.se
#>         <num>     <num>             <char>            <num>               <num>
#> 1: 0.04919468 0.8607673                          0.09719292         0.002805431
#> 2: 0.02602898 0.9727680                          0.11399952         0.019206850
#>    robust.linear.f0.p robust.linear.C0 robust.linear.weights
#>                 <num>            <num>                <char>
#> 1:       0.0006741157        0.3232908           1|1|0.033|1
#> 2:       0.0272319550        0.3525107               1|1|1|1
#>    robust.linear.diagnostics    HMR.f0  HMR.f0.se   HMR.f0.p HMR.kappa
#>                       <char>     <num>      <num>      <num>     <num>
#> 1:                                  NA         NA         NA        NA
#> 2:                           0.2332008 0.01194134 0.03257046  1.628278
#>      HMR.phi   HMR.AIC HMR.AICc     HMR.RSE
#>        <num>     <num>    <num>       <num>
#> 1:        NA        NA       NA          NA
#> 2: 0.5938413 -33.15831       NA 0.002821235
#>                                                    HMR.diagnostics kappa.max
#>                                                             <char>     <num>
#> 1: Missing value or an infinity produced when evaluating the model  2.724784
#> 2:                                                                  3.677404
#>          flux     flux.se       flux.p        method
#>         <num>       <num>        <num>        <char>
#> 1: 0.09719292 0.002805431 0.0006741157 robust linear
#> 2: 0.23320077 0.011941336 0.0325704607           HMR
flux.results[,c(1,26:30)]
#>        ID kappa.max       flux     flux.se       flux.p        method
#>    <char>     <num>      <num>       <num>        <num>        <char>
#> 1:    ID6  2.724784 0.09719292 0.002805431 0.0006741157 robust linear
#> 2:   ID11  3.677404 0.23320077 0.011941336 0.0325704607           HMR

This appends the columns shown above with the selected and recommanded flux estimate and its method used.

Estimate the minimal detectable flux f.detect

The minimal detectable flux relevant for kappa.max selection algorithm depends on the measurement precision of the GC (GC.sd), the chamber size (area A and volume V) and the timing of the sampling scheme (t i.e. at 0, 20, 40 and 60 minutes).

### estimate f.detect by simulation 
C0    <- 325   #ambient concentration, here in [ppm]
GC.sd <- 5 #uncertainty of GC measurement, here in [ppm]
#create simulated concentrations corresponding to flux measurements with zero fluxes:
set.seed(42)
sim <- data.frame(t = seq(0, 1, length.out = 4), C = rnorm(4e3, mean = C0, sd = GC.sd),
                  id = rep(1:1e3, each = 4), A = 1, V = 0.535125)   # specify your sampling scheme t (here in [h]) and chamber volume (V) and area (A)
#fit HMR model:                  
simflux <- gasfluxes(sim, .id = "id", .times = "t", methods = c("HMR", "linear"), plot = FALSE, verbose = F) 
simflux[, f0 := HMR.f0]
simflux[is.na(f0), f0 := linear.f0] # use linear estimates where HMR could not be fitted
#dection limit as 97.5 % quantile (95 % confidence):
f.detect <- simflux[, quantile(f0, 0.975)]
f.detect # here in [ppm/h/m^2], use same unit as your flux estimates, 
#>    97.5% 
#> 20.45867
#e.g., convert to mass flux assuming chamber temperature of 15 degree celcius and standard air pressure
f.detect / 1000 * 28 * 273.15 / 22.4 / (273.15 + 15)
#>      97.5% 
#> 0.02424209

Visualise the flux calculation procedure for a specific GHG chamber measurement system

Simulations are useful for understanding the behaviour of flux calculation algorithms. The performance of a flux calculation scheme depends on the precision of the GC (GC.sd), height of the chamber (V/A), the sampling time scheme (especially the chamber closure time, refered as t.meas), the tightness of the chambers, the reliability of the sampling vial handling etc. This is why it makes sense to look at each specific measurement system and how well it copes with non linear flux estimates (i.e. HMR).

With a simulation we can estimate the impact of the calculation scheme on the flux itself and hence its bias:

and the uncertainty associated with it:

Once having created a simulated flux dataset one can check, where the fluxes are placed in the parameter space of the estimated kappa vs. flux f0.

The code for the example simulation is available as a Gist and further details are described in Hüppi et al. (2018).