It is generally recognized that when chambers are placed on the soil surface, conditions are altered so that the flux of gas is affected. Buildup of gas in the chamber headspace and soil pores reduces the diffusive flux of gas from the soil surface when the chamber is in place, resulting in chamber headspace gas concentrations vs. time data that are nonlinear (Hutchinson and Mosier, 1981; Livingston and Hutchinson, 1995). Application of linear regression to such data results in fluxes that underestimate the actual predeployment flux (Hutchinson and Mosier, 1981; Livingston and Hutchinson, 1995; Healy et al., 1996; Livingston et al., 2006; Rochette and Bertrand, 2007; Venterea and Baker, 2008; Venterea, 2010).
There have been several mathematical algorithms applied to correct for this diffusion effect. Hutchinson and Mosier (1981) developed an equation that uses three sample points collected at equal time intervals (Hutchinson/Mosier [H/M] method). The quadratic procedure described by Wagner et al. (1997) involves fitting a quadratic equation to the concentration vs. time data (Quad method). The flux is then computed as the first derivative of the quadratic equation at time zero. Pedersen et al. (2001) developed a stochastic diffusion model that is an extension of the H/M method and does not require data points determined at equal time intervals and can accommodate more than three data points. The non-steady-state diffusive flux estimator (NDFE) developed by Livingston et al. (2006) is a three-parameter model in which the predeployment flux can be derived from concentration vs. time data by nonlinear regression. Recently, Pedersen et al. (2010) developed a technique designated as the revised Hutchinson/Mosier (HMR) model, which is a modification of the H/M technique to account for horizontal gas diffusion and chamber leaks. Similar to the Pedersen et al. (2001) stochastic model, the Quad method, and the NDFE model, the HMR technique can be used with data sets of three or more points. The common goal of all of these techniques is elimination of the bias associated with the assumption that the concentration vs. time data are linear. However, there are statistical properties other than bias associated with estimators.
The variance associated with different flux estimation methods is influenced by the variability in the concentration vs. time data. This variance affects the minimum detection limit associated with a given flux calculation technique. In this study, we used Monte Carlo sampling to determine the positive and negative flux detection limits associated with different flux calculation procedures at a Type I error rate of 0.05. These assessments were done over a range of analytical precisions and chamber deployment times and for two sampling intensities (three or four points per flux data set). For the purposes of illustration, in this work we use N2O as an example. We also scaled the N2O results so they would be gas species independent. Furthermore, the minimum detectable fluxes of each flux calculation method for other gas species can be determined if analytical precision, chamber DT, and ambient concentration of the gas are known.
Materials and Methods
Sampling Precision of Gas Chromatography
Sampling and analytical precisions of the gas chromatographic measurements of N2O, CH4, and CO2 at ambient concentrations were determined by calculating the SDs and CVs from 35 air samples. Air samples (11.4 mL) were collected in a 10-mL polypropylene syringe and injected into evacuated glass serum vials (empty volume ∼10.5 mL), which were capped with gray butyl rubber stoppers (Voigt Global). In the laboratory, the samples were analyzed for N2O, CH4, and CO2 using a gas chromatograph (GC) (model 8610C, SRI Instruments). An autosampler similar in design to that described by Arnold et al. (2001) was connected to the GC to facilitate sample injection via a sample valve with a 1.0-mL sample loop. Gas species separation was accomplished with stainless steel columns (0.3175 cm diameter × 74.54 cm long) packed with Haysep D and contained in the GC column oven operated at 50°C. Nitrous oxide was detected with an electron capture detector operated at 325°C. Methane and CO2 were measured with a methanizer interfaced with a flame ionization detector operated at 350°C. The carrier gas (N2) flow rate through the column was 20 mL min−1. Certified standard gases (±5%) were obtained from Scott Specialty Gas and used to generate the relationships between detector voltage output and gas concentrations. Precisions of the gas chromatographic analyses of N2O, CH4, and CO2 were determined by computing the means and standard deviations of the 35 measurements of each gas species in air. Precision of measurement of each gas species is expressed as its CV as proscribed by American Public Health Association (1985). Normality of the distributions of the concentrations of the gas species and of the distributions of calculated fluxes was determined by the Kolmogorov Smirnov test. Symmetry of the distributions was assessed by calculating skewness and by examining the trends exhibited by the midsummaries of the ordered fluxes as described by Emerson and Stoto (1983). It was important to determine the distributional properties of the GHGs of interest to guide construction of the populations for the Monte Carlo samplings.
Monte Carlo Simulations: Limit of Detection of Gas Fluxes
Flux determinations using non–steady-state soil chamber techniques typically rely on discrete samples collected from a chamber headspace over a fixed time interval. The flux is then calculated by determining the change in gas concentration vs. time relationship by a curve fitting procedure. However, sampling and analytical error contribute to uncertainty in the gas concentration measurements at each point in time. For example, Fig. 1 shows hypothetical N2O flux determinations based on measurement of chamber headspace gas concentrations at three points in time. The data points in Fig. 1 show the measured concentrations at each time, and the bell-shape curves represent the fact that each discrete measurement is drawn from a population defined by the mean concentration and the measurement variability. In these examples, although the mean chamber headspace N2O concentration does not change over time (it remains constant at 322 nL L-1 [ppb] as represented by the dashed horizontal line), sampling and analytical variability (as represented by the bell-shaped normal distribution curves) result in potential gas measurements that could indicate an apparent positive or negative flux. Reporting these apparent positive or negative fluxes as significant would be committing a Type I error (i.e., rejecting the null hypothesis [H0; flux = 0]) when in fact the flux does equal 0).
The impacts of sampling and analytical variability on the flux threshold value associated with a Type I error rate of 0.05 (i.e., the flux detection limit) for the different methods were evaluated by constructing gas concentration distributions with different SDs (i.e., different sampling and analytical precisions) and by sampling these distributions using Monte Carlo analysis. For these evaluations, scenarios were established whereby distributions corresponding to trace gas concentration measurements at distinct points in time were generated (Fig. 2). In these analyses, the means of the distributions of each gas remained constant (i.e., no change in trace gas concentration over time); therefore, the underlying flux is 0.
The Monte Carlo simulations were performed by generating variates from the unit normal distribution (mean, 0; SD, 1) using the Box-Muller algorithm (Box and Muller, 1958) and the Microsoft Excel @RAND function. Normal variates from distributions corresponding to trace gas concentrations were then generated by multiplying the unit normal variate by the SD of the target population and adding the mean of the target population.
Samples were selected from the distributions at each time point, and a flux was calculated using several calculation techniques (described below). For all the scenarios, the means of the distributions were the same, but different scenarios had different SDs to represent different sampling and analytical precisions. For each scenario, 100,000 Monte Carlo simulations were run, and 100,000 flux estimates were calculated using each of the calculation methods (except the for HRM method). Due to the computationally intensive nature of the nonlinear regression used to evaluate the HMR model, only 10,000 Monte Carlo simulations were performed for the HMR technique over the ranges of analytical precisions and sampling intensities. In some cases, the populations of fluxes generated by the Monte Carlo simulations were significantly different from the normal distribution, so empirical cumulative probability density functions were constructed. These were done by ranking the fluxes in ascending order. The probability associated with each flux measurement was then calculated by dividing its rank number by the total number of fluxes represented in the cumulative distribution. From each cumulative probability density function, the flux corresponding to the 95th percentile was deemed to be the detection limit of the positive fluxes, and the flux corresponding to the fifth percentile was deemed to be the detection limit of negative fluxes (the cumulative probability density functions were approximately symmetrical around zero). These 95th and fifth percentile fluxes represent the 5% Type I error threshold fluxes at the positive and negative tails of the distributions.
The scenarios were run at sampling and analytical precisions (CVs) of 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.08, 0.10, and 0.12. Because chamber deployment time and chamber sampling frequency also affect the calculated flux, simulations were run at simulated chamber deployment times of 0.1, 0.25, 0.5, 0.75, 1.0, and 2.0 h using both three and four equally spaced time points for each CV level. For illustrative purposes, the populations used in each Monte Carlo sampling scenario were generated to represent N2O emissions in that they had means of 320 nL L−1, which correspond to the ambient atmospheric air N2O concentration. However, the results obtained were scaled and thus are applicable to other gas species with differing ambient concentrations.
Flux Calculation Procedures
The flux calculation methods evaluated were linear regression (LR), the Hutchinson and Mosier (1981) (H/M) technique, the quadratic method (Quad) (Wagner et al., 1997), and the HMR procedure (Pedersen et al., 2010). Linear regression and Quad fluxes were calculated using the Microsoft Excel LINEST function (Venterea et al., 2009). The H/M method was implemented as described in Eq. 1 when three gas time points were used.where fo is the calculated flux (units of nL L−1 h−1); C0, C1, and C2 are the chamber headspace gas concentrations (nL L−1) at times 0, 1, and 2, respectively; and t1 is the interval between gas sampling points (hours).
Implementation of the H/M technique requires that (i) only three time point gas concentrations are used (C0, C1, and C2) and (ii) the time interval between sample C0 and C1 must equal the time interval between C1 and C2. The gas concentration units are nL L-1. When fluxes were calculated from four gas time points (C0, C1, C2, and C3), the H/M technique was modified such that the average of gas concentrations of the two intermediate time points (C1, C2) was used (Eq. ).where fo is the calculated flux (units of nL L−1 h−1), C0 is the headspace concentration at time 0, CA1,2 is the average of the headspace concentrations at time C1 and C2, and C3 is the chamber headspace gas concentration at time 3. The term “tA1,2” is the time interval corresponding to the average of time 1 and time 2 (or one half of the total chamber deployment time).
The H/M model can only be evaluated if the quantity [(C1 − C0)/(C2 − C1)] > 1 or if [(C1 − C0)/(C2 − C1)] has a value between 0 and 1 in the case of three time points. In the case of four time points, a H/M flux can only be calculated if [(CA1,2 − C0)/(C3 − CA1,2)] > 1 or if [(CA1,2 − C0)/(C3 − CA1,2)] has a value between 0 and 1. In our simulations, these conditions occurred approximately 40% of the time. In calculating the empirical cumulative frequency distributions and corresponding fifth and 95th percentiles, only the nonfailure fluxes were used.
The HMR procedure was implemented using the R algorithm (R version 18.104.22.168, R, 2010) described by Pedersen et al. (2010) using the following command line: HMR(‘filename’, FollowHRM = TRUE, LR.always = TRUE). The FollowHMR = TRUE statement bypasses manual screening and selection of the appropriate flux by the user. It is recommended that, in application of the HMR procedure, manual selection of the flux be performed by the user, and the FollowHMR = TRUE option should be used in model verification analyses, such as the one presented in this study (A.R. Pedersen, personal communication). The HMR software failed to produce a flux calculated by the HMR model in approximately 7200 out of the 10,000 Monte Carlo samplings for each time point/precision/deployment time scenario. In these cases, the HMR software provided a flux computed by linear regression, or a “no flux” value of zero. Only the HMR-calculated (as identified by the software package) fluxes were used in computing the empirical cumulative frequency distributions and the corresponding fifth and 95th percentiles associated with the HMR model. Also, the HMR procedure was only analyzed with four-time-point data series.
Restricted Quadratic and Hutchinson/Mosier Procedures
Chamber gas concentration data often exhibit a convex upward pattern (for gas production processes) or a concave downward pattern (for gas consumption processes). These are the Type 2 and Type 5 data patterns discussed by Anthony et al. (1995) and may indicate that processes other than gas diffusion may be affecting the observed chamber headspace gas concentrations. In the absence of sampling and analytical error, changes in biological activity (production or consumption) over the course of the chamber deployment period could be responsible for such data patterns. In these cases, application of the H/M, Quad, or HMR procedures underestimate the flux, whereas LR provides a flux estimate that is more representative of the average soil gas flux during the chamber deployment period. The criterion used by Venterea et al. (2009) identifies these occurrences. For three-point data sets, this rejection criterion is: 0 ≤ [(C1 − C0)/(C2 − C1)] ≥ 1, where C0, C1, and C2 are the chamber headspace concentrations at times 0, 1, and 2, respectively. For four-point data sets, this rejection criterion is: 0 ≤ [(C1.2 − C0)/(C3 − C1,2)] ≥ 1, where C0, is the concentration at time 0, C1,2 is the mean of the gas concentrations at times 1 and 2, and C3 is the chamber headspace concentration at time 3. This rejection criterion was applied in additional evaluations of the Quad procedure and the H/M method (designated restricted quadratic [rQuad] and restricted Hutchinson/Mosier [rH/M], respectively). Table 1 summarizes the six flux calculation procedures evaluated in this study.
|Method no.||Method description†||Application description|
|1||Linear regression (LR)||use all 100,000 calculated fluxes|
|2||Quadratic method (Quad)||use all 100,000 calculated fluxes|
|3||Hutchinson and Mosier (H/M)||use only noncomputational-failure fluxes (∼40,000 fluxes for each Monte Carlo sampling)|
|4||Pedersen HMR (HMR)||use only fluxes that the software indicated were calculated by the HMR model (∼2800 fluxes for each Monte Carlo sampling scenario)|
|5||restricted quadratic (rQuad)||use only Quad-calculated fluxes when data did not meet the rejection criteria (∼50,000 fluxes for each Monte Carlo sampling scenario)|
|6||restricted Hutchinson and Mosier (rH/M)||use H/M-calculated fluxes when data did not meet the rejection criteria (∼17,500 fluxes for each Monte Carlo sampling scenario)|
Measurements of N2O, CO2, and CH4 concentrations of replicate air samples enabled calculation of means, SDs, and sample histograms (Fig. 2). The three gases were quantified with the same gas chromatograph, but all had different sampling precisions. The mean measured N2O concentration was 322 nL L−1 with an associated SD of 14.2 nL L−1, resulting in a sampling and analytical precision (CV) of N2O measurement of 0.044 at ambient concentration (Fig. 2A). The analytical variability associated with CO2 measurement (Fig. 2B) was lower than N2O or CH4 (Fig. 2C). The CO2 analysis precision was 0.014, whereas that of CH4 was 0.071. A normality test of the sample distributions indicated that none of the distributions was significantly different from a normal probability distribution at the probability levels (P) indicated in each figure panel.
If soil gas flux in the field were zero, then gas sampling of a chamber headspace over time would in essence be sampling ambient concentrations of the gas of interest over time. For example, in Fig. 3 we randomly selected three or four data points from the N2O data set of 35 air samples (presented in Fig. 2A). Plotting these selections as hypothetical time points collected from a soil chamber over a 1-h deployment time illustrates the possible result when a flux is calculated. Estimation of the “flux” was done using the different flux calculation methods. In the case of three time points (Fig. 3A), the calculated fluxes ranged from 77.5 nL L−1 h−1 for the H/M method to 14.0 nL L−1 h−1 for LR. Fluxes are expressed as nL L−1 h−1 to preserve generality and could be converted to units of moles (or mass) per unit area per unit time for any given chamber volume-to-area ratio and air temperature. Fluxes calculated from four random points selected from the distribution of measured N2O concentrations in air ranged from −76.0 (Quad method) to 15.1 nL L−1 h−1 (LR). The H/M procedure failed to produce a flux estimate for the data set of Fig. 4B because the quantity ln [(CA1,2 − C0)/(C3 − CA1,2)] = ln (−5.07) is not defined. Similarly, the HMR procedure produced a “no flux” recommendation, so a HMR-model flux was not reported. Because the data used to compute these fluxes were drawn from the same population of ambient N2O concentrations, considering any of these fluxes to be significantly different from zero would be committing a Type I error. This is the principle used in conducting the Monte Carlo simulations.
The Monte Carlo samplings of populations of N2O concentration enabled the generation of populations of N2O fluxes for each calculation method. The means and SDs of the populations of fluxes generated at each sampling precision for each computation method and chamber deployment time were calculated. Table 2 illustrates the results obtained when four time points were used with a chamber deployment time of 0.75 h. In theory, the means of the populations of fluxes should be zero over the entire range of analytical precisions, and the means associated with LR, Quad, H/M, rQuad, and rH/M range from −1.43 to 1.055. In contrast, the means associated with the HMR procedure range from −275,639 to 115,816; however, because of their large associated standard deviations, they are not significantly different from zero. The nonlinear regression used in the HMR procedure appears to be sensitive to small deviations in the data and occasionally produces extremely large or small flux estimates. For example, a data set of N2O concentrations of 297.1, 323.7, 305.7, and 314.1 nL L−1 (at times of 0, 0.333, 0.667, and 1.0 h, respectively) yields a HMR-recommended flux of 1.19e+08 nL L−1 h−1. However, if the first value in this data series is changed from 297.1 to 297.0 nL L−1, the HMR software provides a recommendation of “no flux.” In evaluating of the HMR method, we used all the HMR model–derived fluxes—even the extreme values—when they were recommended by the software. However, we observed that these extreme values occurred approximately 0.7% of the time (∼0.35% negative and 0.35% positive extreme fluxes). Thus, although these extreme fluxes affect the means and SDs of the populations of fluxes, their effect on the flux values associated with the fifth and 95th percentiles is <2%.
|nL L−1 h−1|
|0.02||−0.010 (11.5) ‡||0.011 (40.0)||0.294 (37.5)*||−0.230 (47.7)*||0.456 (49.4)*||−40,499 (1,874,000)*|
|0.04||−0.024 (22.9)||0.460 (80.3)||−0.219 (74.3)*||0.102 (95.6)*||−1.53 (96.2)*||−120,910 (4,176,000)*|
|0.06||−0.002 (34.4)||0.550 (121)||−0.294 (112)*||−1.34 (143.4)*||−0.150 (145.7)*||−275,639 (9,846,000)*|
|0.08||0.122 (45.9)||−0.401 (161)||0.325 (148)*||0.782 (191.2)*||−0.229 (198.8)*||−19,540 (8,979,000)*|
|0.10||−0.074 (57.2)||0.670 (201)||−1.46 (188)*||0.644 (239.2)*||0.728 (246.5)*||115,816 (11,120,000)*|
|0.12||0.231 (68.9)||1.055 (241)||0.957 (221)*||0.657 (286.9)*||0.127 (294.4)*||−238,600 (14,090,000)*|
The SDs of the populations of fluxes are influenced by analytical precision. With increasing CV (decreasing analytical precision), the SDs associated with the populations of Monte-Carlo derived fluxes increase. Chamber deployment time also influences the variability of the fluxes. For a given analytical precision, as chamber deployment time increases, SD decreases. For example, at a CV of 0.06, the SDs of the N2O flux populations derived from LR are 51.5, 34.4, and 25.7 nL L−1 h−1 for chamber deployment times of 0.5, 0.75, and 1.0 h, respectively (data for 0.5 h and 1.0 h times not shown). This pattern of decreasing standard deviation with increasing deployment time exists for all six flux calculation methods. In all cases, the populations of fluxes derived from LR and the Quad method are not significantly different from a normal distribution (P > 0.2); however, the distributions of fluxes generated by H/M, rQuad, rH/M, and the HMR procedure are significantly different from normality (P ≤ 0.001).
Because normality could not be assumed for each flux population, cumulative probability density curves were derived empirically. Figure 4 shows an example of the cumulative probability density curves for the six different flux estimation procedures evaluated using four-time-point data sets and a chamber deployment time of 0.75 h for a range of analytical precisions. For all the methods, the populations of fluxes are symmetric around zero. This is to be expected because (i) at each analytical precision, the concentration data points used to compute the fluxes were drawn from the same population; thus, on average, the concentration change vs. time should equal zero; and (ii) all the methods are capable of yielding positive and negative flux estimates. The cumulative probability density curves were generated for each analytical precision and deployment time and were used to compute the detection limit range for each flux calculation method. The example presented in Fig. 5 is the empirical cumulative probability density curves for the six flux calculation methods evaluated at an analytical precision of 0.04 and a chamber deployment time of 0.75 h for the different flux calculation methods.
Because the cumulative distributions are symmetric around zero, the fluxes corresponding to the 5% probabilities at the upper and lower tails were used to determine the upper and lower detection limit ranges at an α level of 0.05. This procedure is illustrated in Fig. 6, where the boxes at the upper and lower tails of the cumulative probability curves have been expanded. At the lower tail (Fig. 6A), the N2O fluxes corresponding to the 0.05 probability level, which were determined from the intersection of each curve and the horizontal line at P = 0.05, are −37.7, −114.7, −131.4, −149.1, −153.2, and −232.5 nL L−1 h−1 for LR, H/M, Quad, rQuad, rH/M, and the HMR methods, respectively. These fluxes represent the negative flux detection limits of the respective flux calculation procedures (α = 0.05). At the upper tail (Fig. 6B), the positive N2O flux detection limits are 228.2, 151.2, 149.1, 131.7, 111.3 and 37.6 nL L−1 h−1 for the HMR, rQuad, rH/M, Quad, H/M, and LR methods, respectively.
The positive and negative detection limits for each flux calculation method were determined for each level of analytical precision, deployment time, and sampling intensity as described above. For each method and for each deployment time, the magnitudes of the positive and negative detection limits were linearly related to analytical sensitivity. Examples of the precision vs. flux detection limit relationships for each flux calculation procedure are shown in Fig. 7. These relationships are presented for a chamber deployment time of 0.75 h and sampling intensities of three and four time points. In all cases, at a given sampling precision, the linear regression method for flux calculation has the narrowest limit of detection window. At a sampling intensity of three points, the rH/M method has the widest detection limit window, and when four time points are used in flux calculation, the HMR method has the widest detection limit window.
The slopes of the regression lines for the precision vs. detection limit relationships for all the deployment times tested are shown in Table 3 (three time points) and Table 4 (four time points). In all cases, the regressions yielded y intercepts that are not significantly different from zero (P > 0.1), and the regression coefficients are >0.999. The positive and negative slope factors presented in Tables 3 and 4 can be used to calculate positive and negative N2O flux detection limits for any given level of precision (for deployment times of 0.1, 0.25, 0.5, 0.75, 1.0, and 2.0 h). For example, for a sampling intensity of four time points and a deployment time of 1.0 h, the slope factor corresponding to the positive detection limit of the Quad method is 2479 (Table 4). If the sampling precision (CV) for N2O is 0.065, the positive detection limit of the Quad method would be 2479 × 0.065 = 161.1 nL L−1 h−1.
|Deployment time||Negative limit of detection slope
||Positive limit of detection slope
|h||nL L−1 h−1 CV−1|
|Deployment time||Negative limit of detection slope
||Positive limit of detection slope
|h||nL L−1 h−1 CV−1|
The conversion factors for computing flux detection limits shown in Tables 3 and 4 are only valid for N2O at ambient levels of 320 nL L−1. However, these results were scaled, so they are applicable to any gas at any ambient concentration. Scaling was done by dividing the slope factors by the 320 (the ambient N2O concentration used in the Monte Carlo simulations). Tables 5 and 6 show the scaled slope factors associated with each flux calculation method over the range of deployment times tested. The scaled slope factors can be used to compute flux detection limits for any gas if the ambient concentration of the gas of interest is known. For example, consider the scenario where three time points are collected over a chamber deployment time of 0.75 h and an ambient CH4 concentration of 1.75 μL L−1. From Table 5, it is determined that the scaled slope factors defining the precision vs. flux detection limits relationships for LR are −3.10 (negative flux detection limit) and 3.100 (positive flux detection limit). Conversion of these scaled slopes to factors to CH4 is done by multiplying by the ambient CH4 concentration. Thus, the negative and positive slope factors for CH4 flux detection at a chamber deployment time of 0.75 h are −5.425 and 5.425, respectively. If the analytical precision for CH4 detection at ambient levels is 0.096, then the positive and negative detection limits for CH4 fluxes calculated using LR would be −5.425 × 0.096 = −0.521 μL L−1 CH4 h−1 and 5.425 × 0.096 = 0.521 μL L−1 CH4 h−1, respectively. The scaled slope factors shown in Tables 5 and 6 are only valid for the chamber deployment times listed.
|Deployment time||Negative limit of detection factor
||Positive limit of detection factor
|Deployment time||Negative limit of detection scaled factor
||Positive limit of detection scaled factor
To extend the applicability of these results to other chamber deployment times, we developed regression relationships between deployment time and scaled slope factor (Fig. 8). Because the positive and negative slope factors for a given flux calculation procedure and deployment time are similar in magnitude and differ only in sign, the absolute values of the scaled slope factors were regressed against deployment time using the exponential model θ = a × DT−b, where θ is the scaled slope, DT is chamber deployment time, and a and b are regression coefficients. For all the methods, there is a significant exponential relationship between chamber deployment time and scaled slope factor. Table 7 presents the results of these regression analyses.
|Flux calculation procedure‡||Regression coefficients
|Three sampling points|
|Four sampling points|
The perturbations to soil gas flux resulting from placement of a chamber on the soil surface are well documented (Hutchinson and Livingston, 1993; Livingston and Hutchinson, 1995). Increases in chamber headspace gas concentration can slow the transport of gas from the soil, resulting in nonlinear concentration vs. time data. As a result, application of linear regression to diffusion-affected chamber data results in underestimates of trace gas flux (Hutchinson and Livingston, 1993; Livingston and Hutchinson, 1995). There have been numerous mathematical procedures developed to compute unbiased fluxes from curvilinear data (Hutchinson and Mosier, 1981; Wagner et al., 1997; Pedersen et al., 2001; Livingston et al., 2006; Pedersen et al., 2010). Other than bias, relatively little attention has been given to other statistical properties of flux estimation techniques. A notable exception is a recent study of Venterea et al. (2009). These investigators examined the influence of measurement variability on the bias and variance of flux estimates obtained from linear regression and from two nonlinear methods (Hutchinson and Mosier, 1981; Wagner et al., 1997). The variance of an estimator is a reflection of the estimator's sensitivity to the variability inherent in the underlying data. The degree of sensitivity will, in turn, affect the minimum detectable flux.
Past efforts to assess the minimum detection limits of soil gas emissions have focused on determining goodness-of-fit of regression procedures. For fluxes determined by linear regression, a t test of the slope of the regression line can be used to assess if the flux is significantly different from zero (Livingston and Hutchinson, 1995; Rochette et al., 2004). Because standard errors of the model parameters obtained in the Quad and HMR methods can also be calculated, a t test of significance can be applied to determine the significance of fluxes derived by these methods. The H/M flux procedure does not allow for calculation of an associated standard error directly. However, the stochastic application of the H/M procedure developed by Pedersen et al. (2001) does provide flux estimates with associated confidence limits, enabling the determination of regression significance. Although the above-mentioned statistical tests indicate whether a given flux is significantly different from zero, they do not provide an indication of the magnitude of the minimum detectable flux.
Hutchinson and Livingston (1993) describe a procedure for computing the minimum detection limit of individual chambers based on the standard error for each chamber (Eq. ):where t is the t statistic for n − 2 degrees of freedom at the user-desired probability level, and SE is the standard error of gas exchange rate for each chamber.
Although this approach attempts to include all the errors associated with a given chamber flux measurement (not just analytical errors), the practical utility of this technique may be limited. For each chamber and for each individual flux measurement, a different SE is derived. Thus, the minimum detectable flux will be different for each chamber each time a flux is measured. We applied this technique in our Monte Carlo study, and, for linear regression with three time points, we observed that Eq.  yielded minimum detectable flux values that ranged from 1945 nL L−1 h−1 to 0.0048 nL L−1 h−1 for an analytical precision of 0.05 and a chamber deployment time of 0.75 h.
Verchot et al. (1999) recognized the limitations associated with the minimum flux detection calculation identified in Eq.  and implemented a modification of the regression significance procedure whereby all their measured fluxes were ranked. The minimum detection limit was then determined as the flux threshold where >67% of the measured fluxes were significant. These authors reported a minimum detectable flux for N2O emission of 0.6 ng N cm−2 h−1. To compare this value with our results, we used the Verchot et al. (1999) chamber dimensions, deployment time, and sampling intensity to calculate that a flux of 0.6 ng N cm−2 h−1 corresponds to a headspace concentration change of detection limit of 29.9 nL L−1 h−1 (at 25°C and 1 atmosphere pressure). Results of our study indicate that this detection limit corresponds to a measurement precision of approximately 0.02 (at four sampling points, a 0.5-h deployment, and an ambient N2O concentration of 320 nL L−1). The approach adopted by Verchot et al. (1999) is more defensible than simply assessing regression goodness of fit; however, the assumption that 67% of the significant fluxes in any given data set exceed the actual minimum detectable flux is arbitrary.
Assessment of the minimum detectable flux could be determined experimentally by placing a chamber on a non–trace gas producing surface (e.g., a linoleum floor), sampling the headspace gas concentration of the chamber over time, and computing an apparent flux. However, this process would have to be repeated extensively to obtain reliable probability estimates associated with each flux determination. Monte Carlo sampling allows for the same analysis to be conducted mathematically and facilitates assessment of the effects of measurement precision on detection limit. In our evaluations we minimized the Type II error rate associated with the regression analyses by imposing the condition that all the regression-derived fluxes (LR, Quad, rQuad, and HMR) are significantly different from zero. Thus, the flux detection limits calculated by our procedure can be considered conservative estimates.
We observed linear relationships between analytical precision and minimum detectable flux. In addition to analytical and sampling precision, chamber deployment time, sampling intensity, and flux calculation method affect the detection limits associated with chamber headspace gas concentration changes over time. The regression relationships shown in Table 7 summarize these effects and can be used to compute a minimum detection limit for any gas.
As an illustration, we used the mean ambient concentrations and the associated sampling and analytical precision estimates shown in Fig. 2 to calculate the minimum detectable fluxes for N2O, CH4, and CO2 determined by the Quad method (Table 8). The first step in implementing this procedure is to calculate θ (the scaled slope factor using the coefficients of the exponential equation associated with the rQuad method with three time points). The slope factor for each gas is then calculated by multiplying θ by the ambient concentration of the gas of interest. The resulting slope factors for each gas are multiplied by the sampling and analytical precision of each gas to obtain the positive flux detection limit. Multiplication of the positive flux detection limits by −1 yields the negative flux detection limits. These flux detection limits have units of ηL−1 h−1. Conversion of these values to flux units of mass per unit area can easily be done if the temperature, pressure, and chamber dimensions are known. When this is done, chamber size (specifically height) affects how the ηL L−1 h−1 detection limit is translated into a gas flux value expressed on an areal basis. A detailed description of flux detection limit calculations is provided in the supplementary materials.
|Mean ambient concentration, ηL−1 (N2O); mL L−1 (CH4 and CO2)||323||1.79||385.5|
|Deployment time, h||0.667||0.667||0.667|
|θ (calculated from Table 7, rQuad method, three time points, 0.667 h deployment time)||10.61||10.61||10.61|
|Slope factor (θ × mean concentration)||3428||18.99||4091|
|Positive flux detection limit, ppb h−1 or ppm h−1 (slope factor × CV)||150.8||1.349||5.728|
|Negative flux detection limit, ppb h−1 or ppm h−1 (−1 × slope factor × CV)||−150.8||−1.349||−5.728|
Despite the method used to compute the minimum detection limit (MDL), consideration must be given to how fluxes that fall below the limit are treated. There are several options available to handle data that fall below the MDL (or within the detection limit band). These options are summarized by Gilbert (1987) and include (i) report the value as “below the detection limit,” (ii) report the value as zero, (iii) report some values between zero and the MDL (such as one half the MDL), or (iv) report the actual measured value even if it falls below the MDL. Of these options, Gilbert (1987) recommends the latter as the least biased course of action (i.e., report the measured value along with the stated MDL).
Due to variability in concentration vs. time data, no single flux calculation scheme will be applicable (or optimal) over the entire range of concentration vs. time data series encountered. For example, in cases where (C1 − C0)/(C2 − C1) < 0, the Quad, H/M, or HMR methods may not be applicable. In these cases, LR could be used. Therefore, rather than recommending a single flux calculation method to the exclusion of others, we suggest that “hybrid” flux calculation schemes be adopted whereby a linear regression is used when a given data set does not conform to a nonlinear model. Such a hybrid scheme is resident in the HMR software, which provides a linear flux when the HMR model fails. The coefficients provided in Table 7 can be used in the manner illustrated in Table 8 to compute the MDL for each of the individual flux calculation methods used within such a hybrid scheme.
Our study only considered sampling and analytical precision associated with trace gas concentration measurement. Other sources of variability (e.g., chamber leakage, changes in biological activity during the chamber deployment period) may also reduce measurement precision. However, our assessment of analytical precision mirrors the methodology used to collect chamber gas samples in the field. Our results on the minimum detectable fluxes associated with each flux calculation method should not be the only consideration in selecting a flux calculation method. For any given level of precision and deployment time, linear regression had the narrowest detection window, yet this method is known to yield biased estimates in some situations. However, methods that attempt to correct for diffusion effects on chamber headspace gas concentration data (i.e., Quad, H/M, HMR) may not be appropriate in all data sets. We are currently expanding our investigations of the bias and variance associated with different flux calculation methods to identify specific criteria that will enable recommendation of one technique over another. Finally, because of the increasing importance being placed on estimates of soil GHG emissions, we recommend that reports of soil trace gas flux include information about sampling and analytical precision and associated estimates of minimum detectable fluxes.