The noble gases, which are chemically inert under normal terrestrial conditions but vary systematically across a wide range of atomic mass and diffusivity, offer a multicomponent approach to investigating gas dynamics in unsaturated soil horizons, including transfer of gas between saturated zones, unsaturated zones, and the atmosphere. To evaluate the degree to which fractionation of noble gases in the presence of an advective–diffusive flux agrees with existing theory, a simple laboratory sand column experiment was conducted. Pure CO2 was injected at the base of the column, providing a series of constant CO2 fluxes through the column. At five fixed sampling depths within the system, samples were collected for CO2 and noble gas analyses, and ambient pressures were measured. Both the advection–diffusion and dusty gas models were used to simulate the behavior of CO2 and noble gases under the experimental conditions, and the simulations were compared with the measured depth-dependent concentration profiles of the gases. Given the relatively high permeability of the sand column (5 × 10−11 m2), Knudsen diffusion terms were small, and both the dusty gas model and the advection–diffusion model accurately predicted the concentration profiles of the CO2 and atmospheric noble gases across a range of CO2 flux from ∼700 to 10,000 g m−2 d−1. The agreement between predicted and measured gas concentrations demonstrated that, when applied to natural systems, the multi-component capability provided by the noble gases can be exploited to constrain component and total gas fluxes of non-conserved (CO2) and conserved (noble gas) species or attributes of the soil column relevant to gas transport, such as porosity, tortuosity, and gas saturation.
Thermal gradients, gravity, and the production and depletion of gases from biologic, magmatic, or other sources will alter the distribution and isotopic compositions of the indigenous atmospheric gases in unsaturated soils, regolith, firn, or other surface material (e.g., Severinghaus et al., 1996, 2003; Evans et al., 2001; Camarda et al., 2007; Freundt et al., 2013; Jones et al., 2014). Atmospheric noble gases, comprising a mass range from 3 to 136 atomic mass units, have the potential to probe and quantify the various processes and soil properties impacting gas transport (Severinghaus et al., 1996; Jones et al., 2014). Below a few decimeters, shallow soil horizons are quasi-isothermal (thermally buffered). Under these conditions, gas transport is driven by external pressure gradients (e.g., a CO2 or CH4 source, O2 depletion, etc.) resulting in laminar advective transport (Darcy’s law) and/or a gradient in the concentrations of soil gas components that drives molecular diffusion (Fick’s law). To simulate the mechanism of transport, two models are generally applied: the advection–diffusion (AD) and the dusty gas (DG) models. The models differ in formulation in that the DG model explicitly accounts for Knudsen diffusion (Thorstenson and Pollock, 1989a, 1989b). The apparent success of these two models in simulating transport and matching field-measured soil gas data has often been cited as validation of the models (e.g., Jones et al., 2014). However, field studies involve complex natural systems with many unconstrained variables; thus, a controlled, laboratory-based evaluation of the models would more directly test their validity.
Federico et al. (2010) analyzed the fractionation of 13CO2 in a volcanic field using both the AD and DG models. While several studies (e.g., Webb, 1998) have compared the DG and AD models for binary mixtures using laboratory data, laboratory multicomponent comparisons are rare. None, to our knowledge, have used noble gases.
Noble gases have recently been used to investigate processes of gas transport in ice caps (Huber et al., 2006; Orsi et al., 2015). The inert nature of noble gases and their presence at detectable levels in the atmosphere and in soil gas make them ideally suited for laboratory tests of transport models. We conducted a simple laboratory column experiment representing a shallow unsaturated soil horizon and determined the concentrations and relative compositions of noble gases and CO2 as a function of depth in the column for a range of steady CO2 fluxes (∼700–10,000 g m−2 d−1). The measured composition profiles were compared with the profiles predicted by the AD and DG models in the presence of a constant CO2 flux. The data were then used in conjunction with the AD and DG models to explore the utility of noble gases as process tracers of gas transport in soils. Due to the wide range in molecular weights of the noble gases, the present work provides a fairly rigorous test of the multicomponent capabilities of each model.
Experimental Methods and Numerical Simulations
Experimental Apparatus and Methods
A schematic of the laboratory test system is shown in Fig. 1. The laboratory apparatus was the same as the test system described by Evans et al. (2001) and consisted of a standard 208-L (55-gal) polyethylene drum, 56 ± 1 cm in diameter, cut to a height of 40 cm, with a perforated plastic disk supported 7 cm off the bottom. Above the perforated disk, the drum was filled to a height of 30 cm with kiln-dried, 0.595-mm (30-mesh) beach sand with a porosity of 0.40, a bulk density of 1.59 Mg m−3, and a volumetric water content of 0.0017. Carbon dioxide was injected at the base of the column, below the perforated disk, and maintained at a constant flow rate set by an external pressure regulator and capillary restriction. The regulator was calibrated and the stability tested using an in-line flow meter and/or by momentarily diverting the flow to a calibrated water displacement flask. Agreement between the techniques was generally within 2%, but at flows >3000 g m−2 d−1, only the water-displacement flask could be used. The diurnal drift in the flow rate was as much as ±4%, assumed to represent the stability of the regulator (Evans et al., 2001). Laboratory air temperature was controlled at 20 ± 1°C. Laboratory air pressure was maintained about 5 Pa below ambient building pressure as required by state and federal safety regulations. While not directly monitored, short-term fluctuations in laboratory pressure were evident as high-frequency noise in differential pressures as automated air-handling equipment switched on and off. Differential pressures were measured every 5 s at the highest flow rate and every second at lower flow rates, with averages recorded every minute with the same transducer that was cycled through the sampling ports to average out fluctuations in the reference (ambient laboratory air) pressure.
The sand column experiment was run at four different injected CO2 fluxes: 707, 2995, 5079, and 10,135 g m−2 d−1. The injected CO2 was allowed to flow through the column for at least 48 h (lowest flux) and 6 h (highest flux) before measurements, times shown to be sufficient to allow the system to approach steady state (Evans et al., 2001). For example, during a test run at an imposed flux of 507 g m−2 d−1 and 24 h after the initialization of flow, flux at the sand surface was repeatedly measured by a conventional infrared gas analyzer (IRGA)-based accumulation chamber both at the center and around the edges of the sand surface. All 18 flux measurements gave results within 10% of the imposed flux, the average (n = 6) for measurements made at the center was 495 g m−2 d−1, and the average for the edge measurements (n = 12) was 496 g m−2 d−1. This demonstrates a lack of edge effects impacting the flow regime (and see Evans et al., 2001, Fig. 5) and that a near steady-state condition was achieved in times much shorter than the 48-h wait time used for the low-flux experiments.
Sampling ports within the sand were set at depths of 6, 12, 18, and 24 cm below the sand surface and from an outlet tube at a depth of 30 cm from the sand surface and below the perforated plastic disk (Fig. 1). At each sampling depth, following the initial flow period, pressures were measured, gas samples for CO2 and noble gas analyses extracted, and after sample extraction, pressures were measured again. Samples for noble gas analysis were collected in 9.8-cm3 copper tubes sealed with refrigeration-style clamps at each end. The CO2 concentrations and pressure profiles were measured at the same depths as the noble gas samples. The CO2 concentrations (±3%) were determined with a portable IRGA manufactured by Drager and Bacharach, calibrated against seven standards across the 0 to 100% concentration range. Pressure profiles were measured and recorded using the most sensitive commercially available closed-system differential pressure transducer (Setra 265; full-scale range of ±62 Pa) together with a Campbell Scientific datalogger. All pressure readings were taken with the transducer diaphragm positioned at the same elevation as that of the sand surface and the reference port open to ambient laboratory air. The gradient measured in this way represents the pneumatic pressure-head gradient driving the viscous component of flow (Evans et al., 2001), thus accounting for gravitational effects (Thorstenson and Pollock, 1989a, 1989b).
The He, Ne, Ar, Kr, and Xe analyses were performed at the Center for Isotope Geochemistry, Lawrence Berkeley National Laboratory. The Cu-tube samplers were attached directly to a gas purification line, and reactive gases (e.g., CO2, N2, O2) were removed using a Ti-sublimation pump followed by exposure to a Ti–Zr alloy getter held at 300°C. Small fractions of the purified noble gas sample were introduced into a quadrupole mass spectrometer (Pfeiffer QMS 422) for analysis of the elemental abundances and Ar isotopic compositions. The noble gas abundances were determined from peak height measurements of 4He, 22Ne, 36Ar, 84Kr, and 132Xe. Peaks at masses 44 (CO2), 35 and 37 (Cl isotopes), and 18 (H2O) were monitored to correct for doubly charged ion interferences on 22Ne+, Ne isotopic compositions from H2O, and 36,37Ar isotope contamination from H35Cl and H37Cl. The necessary correction factors (production efficiency of doubly charged CO2 relative to singly charged ions and the H35,37Cl/35,37Cl ratios) were determined from mass spectrometer background analyses, and in all cases corrections were less than ∼1%. The analytical procedure and mass spectrometer were calibrated using aliquots of air, the conventional standard for noble gas analyses. During the measurement, 36Ar is the dominant peak, with a peak height ranging from ∼50 to 1000 times the peak heights for the other noble gases. For this reason, 36Ar is the reference isotope of choice and data are reported (Table 1) as F(i) values, defined as the (i/36Ar) ratio measured for a sample normalized to the same ratio measured in the air standard. The latter is determined from an average of multiple analyses (n ≥ 10) of the air standard.
During the runs at 707 and 10,135 g m−2 d−1, we used commercially available CO2 with a purity of 99.5%. However, it became clear that the CO2 was contributing a small amount of noble gases, particularly Kr and Xe. Whenever the CO2 concentration in the sand column exceeded 95%, the large relative amounts of CO2–tank-derived Kr and Xe and their putative fractionation prevented confident correction of those samples for these components. For the other two runs, we used a second source of higher purity CO2 (99.9995%). However, analyses of each CO2 tank found little difference in the Kr and Xe levels of the two different purity tanks (Table 1). Thus samples with CO2 concentration >95% were considered unreliable with respect to Kr and Xe and were treated as such during simulations of soil-gas profiles.
Both the DG and AD models were used to simulate one-dimensional steady-state CO2 and noble gas transport and corresponding concentration profiles within the sand column. Flow geometry in the simulation (Fig. 1b) was defined by a series of volume elements and boundaries between them. Mass balance was attained by requiring that mass accumulation in each volume element to equal the net mass flux across the element boundaries. Initially, with no CO2 flow, the pore space was filled with air at ambient atmospheric concentrations. During the experiment, air at the top of the sand column and pure CO2 at the bottom set the boundary conditions.
The DG model evaluates multicomponent transport simultaneously. By solving a set of equations, the molar diffusive flux of each component can be calculated. The MIN3P (Molins et al., 2008, 2010) software package, which implements the DG model, was used for simulations.
Using the AD model as implemented in TMVOC and the DG model in MIN3P, simulations were run using theoretical binary diffusion coefficients (noble gas–CO2) calculated after Bird et al. (2002). Molecule–molecule interaction parameters, such as collision integrals and characteristics length scales, were estimated using Lennard–Jones potential intermolecular force laws, which approximate potential energy as a function of the distance between atoms to calculate the binary diffusion coefficients in Eq.  (Bird et al., 2002). We found almost no difference between the Bird et al. (2002) binary diffusion coefficients and those calculated from other published sources (e.g., Reid et al. ; see Table 2 for a comparison). The impact of the porous medium attributes is included through a porosity–tortuosity–gas saturation factor applied to the free-air binary diffusion coefficients, as described above.
The pressure measurements and analytical results are provided in Table 1. Pertinent characteristics of the sand column, calculated binary diffusion coefficients, and gas viscosities used in the simulations are provided in Table 2.
As anticipated, the pressure gradient across the sand column was very small, ∼6 Pa at the highest CO2 flux (Table 1), compared with ambient air pressure (∼105 Pa), and at the extreme low end of our measurement capability even with the most highly sensitive transducer we could acquire. For the three experiments with lower applied fluxes, the total range in pressure across the sand column dropped below 3% of the full-scale range of the pressure transducer, and pressure gradients could not be accurately resolved. As a result, only the pressure differences measured under the high-flux condition are presented in Fig. 2 and Table 1. The measured pressure difference between the 30-cm depth and the surface (∼6.3 Pa) is near the values predicted from MIN3P and TMVOC (Fig. 2a). The deviation of measured and theoretical profiles reflects the difficulty of determining exceedingly small differential pressures as well as possible impedance from the perforated plate. Despite these difficulties, the results do indicate a pressure gradient through the sand of about the anticipated magnitude.
Carbon Dioxide Profiles
The vertical profiles of CO2 concentrations as a function of depth and flux are shown in Fig. 3. The depth of each sampling point is well constrained (±2 mm), but additional uncertainty could arise due to the relatively large volume of gas sampled (∼60 cm3). These large sample sizes were needed to stabilize the CO2 concentration measurements and to accommodate the 9.8-cm3 sampler used for noble gas abundances. Assuming that the 60 cm3 was extracted from a spherical volume with a porosity of 0.40 (Table 2), the radius of the sphere corresponds to a depth uncertainty of approximately ±3.3 cm. The vertical error bars in Fig. 3 have been drawn to reflect this value, but if the gas sampled derives from a spherical volume in a concentration gradient, its composition should approximate that at the center of the sphere, i.e., reducing the impact of any depth uncertainty.
The flux-dependent concentration profiles are consistent with those found by Evans et al. (2001) under similar conditions. At the lowest injected CO2 flux (707 g m−2 d−1), the concentration of CO2 increased linearly with depth and was consistent with a total CO2 flux dominated by molecular diffusion. At the highest injected CO2 flux (10,135 g m−2 d−1), the CO2 concentration profile clearly showed the transfer from a total flux dominated by molecular diffusion (linear concentration trend at shallow depths) to an advection-dominated regime (rapidly increasing CO2 concentration gradients at greater depths in the column; Table 1). The experimental data were compared with concentration profiles simulated by the AD and DG models. Although in constructing these curves we were unable to input a measured pressure gradient for each flux (as discussed above), both the MIN3P and TMVOC packages can determine the anticipated pressure gradients from the CO2 fluxes and the properties of the sand column. A good agreement was observed between both the AD and DG models and the experimental data. Chi-squared tests showed that model misfit was statistically insignificant relative to gas-concentration measurement error, with no clear winner between the two models. A detailed statistical analysis comparing the measured and predicted values was not conducted. However, calculated residual values (measured − predicted) suggest that the measured CO2 concentrations were consistently ∼1 to 5 mol% higher than the values predicted by both the AD and DG models, with no correlation with flux. Measured and modeled concentrations are presented in the supplemental material.
Noble Gas Profiles
Figure 4 shows the concentration profiles of the noble gases plotted as fractionation factors [F(i) values], defined as the concentration of the ith gas normalized to 36Ar and the i/36Ar ratio in air: F(i) = [(i/36Ar)sample/(i/36Ar)air]. The F(i) plots are arranged from the lowest to the highest injected CO2 flux. At a given CO2 flux, the noble gas fractionation factors depend on the molecular weight of the noble gas with respect to 36Ar and depth. It is also clear that the degree of fractionation (with respect to the i/36Ar ratio in air) depends on the magnitude of the CO2 flux (note the difference in vertical scales in Fig. 4). With increasing CO2 flux, the noble-gas fractionation profiles become increasingly curved and the magnitude of fractionation increases. To compare the ability of the AD and DG models to predict the noble-gas profiles, the F(4He), F(22Ne), F(84Kr), and F(132Xe) percent residuals, defined as [(measured/predicted) − 1] × 100, are shown in Fig. 5. The F(i) residuals for He and Ne appear to be randomly distributed at any applied flux, suggesting that the AD and DG models fit the data about equally well. The same holds true for Kr except at the high flux. The non-random structure in the two heavier gases at high flux occurred when the CO2 concentrations were >95 mol% (Table 1), reflecting the contribution of CO2 tank-derived Kr and Xe (as discussed above).
The isotopic composition of individual noble gases can provide additional constraints on gas transport. Illustrated in Fig. 6 are variations in Ar isotopic compositions (ratios) as a function of depth in the sand column and the imposed CO2 flux, clearly showing the expected decrease in the heavier 40Ar isotope with depth and CO2 flux. Even at the lowest CO2 flux, the gradient in Ar isotopic ratios (∼10‰ at 24 cm) is easily resolved given the precision of modern Ar isotopic measurements (e.g., Severinghaus et al., 2003). As the injected CO2 flux increases, changes in 40Ar/36Ar ratios with depth become systematically larger. At the highest flux, the 40Ar/36Ar ratio is ∼100‰ lower than the atmospheric value at the 24-cm depth. Comparable shifts in N2 and CO2 isotopic compositions would also be expected (e.g., Severinghaus et al., 1996). Furthermore, it is evident from Fig. 6 that both the AD and DG models can be used to predict the depth- and flux-dependent isotopic variations in the same manner as gas abundances.
Advection–Diffusion and Dusty Gas Models Applied to High-Permeability Systems
The sand column experiment, albeit simple in nature, clearly demonstrates the viability of modeling soil gas transport in binary and multicomponent systems. The model simulations using either the DG or AD models adequately predicted concentration soil gas profiles for multicomponent systems in the presence of varying CO2 fluxes. Although binary systems have been studied previously (Webb, 1998; Oldenburg et al., 2004; Federico et al., 2010), to our knowledge, the multicomponent capabilities of the MIN3P (DG) and TMVOC (AD) software packages have never before been verified experimentally. It is important to note that this demonstration is only valid within the permeability regime of the experiment (∼10−10–10−11 m2). Therefore a full test of the AD and DG models across the spectrum of soil permeability has not been achieved. Nonetheless, within the permeability regime of the experiment, both models appeared to predict CO2 and noble gas compositions within the sand column with comparable skill.
Advection–Diffusion and Dusty Gas Models Applied to Low-Permeability Systems
Given that the Knudsen diffusion coefficient (DK) is a function of permeability, Eq.  can be used to evaluate the relative contributions of the two terms on the left to the total flux as a function of permeability for a given total flux (NT). In Fig. 7, the percentage contribution of the second term on the left (the Knudsen regime) is plotted as a function of permeability for Kr. The curves reflect different sampling depths within the column. At the permeability of the sand column experiment (dashed line), the Knudsen regime contributed <1% of the total flux. The permeability at which the Knudsen regime becomes important (accounts for more than ∼5% of the transport) is on the order of ∼10−13 m2. In low-permeability lithologies, such as dolomites and shales, the Knudsen regime dominates the diffusive flux. The permeability of most unconsolidated soils is likely to fall within a permeability regime where Knudsen diffusion is not a measureable factor. However, Knudsen diffusion is expected to be a factor in clay-rich soil columns or soil columns with interbedded clay layers.
The AD model accounts for Knudsen diffusion by adjusting the intrinsic permeability using the Klinkenberg factor (b, Eq. ). Webb and Pruess (2003) argued that at a permeability where the Knudsen regime begins to dominate (Ko < 10−13), the AD is clearly inadequate unless adjustments to the model are made. One approach has been to add an additional diffusive flux that affects non-equimolar species (e.g., He, Ar,…, Xe). The flux arises because lighter species have greater velocities than heavier ones, which results in an “internal” pressure gradient acting in the opposite direction to the ordinary diffusive flux (e.g., Scanlon et al., 2002). It is clear that to provide a more comprehensive comparison of the DG and proposed adjusted AD models, future column experiments need to be conducted at much lower permeability.
Utility of Noble Gases asConservative Tracers
The effective binary diffusion coefficients have been replaced by their free-air equivalents and the tortuosity (τ), porosity (ϕ), and degree of gas saturation (Sg) have been inserted using the relationship given in Eq. . In the column experiment, the sand was dry; therefore Sg has a value of unity and can be ignored. It is clear from Eq.  that if the total flux is known (e.g., measured at the surface), then the measured F(i) values at depth z provide an in situ measure of the soil attributes (τϕSg) that impact gas transport. Determining F(i) depth profiles can provide a vertical profile of the soil attributes. Furthermore, collecting samples at multiple depths provides additional constraints for determining in situ soil attributes even for the case of much lower permeability where the assumptions used to simplify Eq.  are no longer valid. Alternatively, if the soil attributes are known, the depth-dependent F(i) values provide a measure of the total flux.
The preceding discussion provides an illustration of how noble gases can constrain parameters influencing gas transport in soils. The utility of the noble gases lies in the fact that there are four stable noble gases (Ne, Ar, Kr, and Xe), covering a mass range from ∼20 to 132 atomic mass units that behave coherently in natural systems. As the one exception, 4He is excluded because of a potential source in soil and groundwater systems from radioactive decay of U and Th and therefore will not always meet the criteria defining stagnant gases in soil horizons.
The development above did not take gravitational settling into account due to the short height of the experimental column. For deep unsaturated zones, the effects of gravitational settling will be superimposed on those of any non-conservative flux (Severinghaus et al., 1996).
We also note that the F(132Xe) values in Fig. 5 do not conform well with the plotted curves and that the deviations occur when the CO2 concentrations are >95 mol%. This is an experimental artifact reflecting contamination of the ambient-air-derived noble gases with CO2–cylinder-derived Xe, as discussed above.
Within the high-permeability regime of the sand column CO2 flux experiment (5 × 10−11 m2), the DG model (MIN3P software package) and the AD model (TMVOC software package) accurately predicted the depth-dependent concentration profiles of the non-stagnant CO2 and the stagnant noble gases across a range of CO2 flux from ∼700 to 10,000 g m−2 d−1. To our knowledge, this is the first time the multicomponent capabilities and software packages of both transport models have been laboratory tested and confirmed. The primary difference between the two models is that the AD model does not explicitly include the impact of gas molecule–particle (sand grain) interaction on gas diffusivity, referred to as Knudsen diffusion. Molecule–particle interaction becomes important at very low permeability, and at the permeability of the sand column, Knudsen diffusion accounted for <1% of the total diffusive flux (Fig. 7). For a more rigorous evaluation of the DG model and modified versions of the AD model that take into account Knudsen type diffusivity, column experiments at much lower permeability need to be conducted. For instance, in column experiments with a permeability of ∼10−16 m2, Knudsen diffusion is expected to account for ∼70% of the total flux.
When applied to natural systems, the multicomponent capability built into the DG model and added to the AD model through the TMVOC software package can be exploited to constrain the total flux of non-conserved species, such as CO2, and attributes of the soil column relevant to gas transport, such as porosity, tortuosity, and gas saturation. In the presence of a gas flux driven by a source or sink, conserved species in the soil column, including the noble gases, will redistribute themselves in accordance with their respective mass-dependent advection and diffusion velocities. The resulting mass fractionation with respect to the initial compositions provides multiple equations that can be solved for either total flux or soil attributes. In the case of the noble gases, after excluding 4He (which may be locally radiogenic), an initial atmospheric composition in the soil column can be safely assumed for the four conserved species (Ne, Ar, Kr, and Xe) to yield three equations defining their relative fractionation with respect to the atmosphere. For a single non-conserved species like CO2, the system is well constrained.
Natural systems are more likely than not to be very complex, with multiple non-conserved species added and/or removed from the soil column. The present work provides experimental confirmation for the approach taken by Jones et al. (2014), who applied DG modeling of noble gases to deconvolve subsurface biogeochemical and transport processes associated with hydrocarbon degradation at an oil-spill site. The ability to constrain both the processes and rates of gas transport and generation in complex systems is an excellent example of how the multicomponent conserved species, such as the noble gases, can provide valuable information regarding gas dynamics in vadose and saturated zones.
Variables and Units
Special thanks to Stefan Finsterle and Sergi Molins for providing access to the MIN3P and TMVOC software packages and educating us in their use. The collective experience and insight of Stefan and Sergi was invaluable. This work was supported by the Director, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, of the U.S. Department of Energy under Contract no. DE-AC02-05CH11231. USGS reviewer Chris Green offered helpful suggestions for improving the paper, as did Associate Editor Peter Nico and three anonymous reviewers. We acknowledge support from the National Research Program and Toxic Substances Hydrology Program of the USGS. Mention of trade names is for identification purposes only and does not constitute endorsement by any entity mentioned herein.