Published online 26 April 2005
Published in Vadose Zone J 4:225-239 (2005)
DOI: 10.2136/vzj2004.0029
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Measurement and Modeling of Temporal Variations in Hydrocarbon Vapor Behavior in a Layered Soil Profile
G. B. Davis*,
J. L. Rayner,
M. G. Trefry,
S. J. Fisher and
B. M. Patterson
CSIRO Land and Water, Private Bag 5, Wembley, W.A., 6913, Australia
* Corresponding author (Greg.Davis{at}csiro.au)
Received 29 January 2004.
 |
ABSTRACT
|
|---|
Fine-scale measurement of gasoline vapors, major gases (O2, CO2, N2, and CH4), residual nonaqueous phase liquid (NAPL) gasoline, and soil physical properties has allowed detailed assessment of the role of soil layering and seasonal variability on hydrocarbon vapor fate and biodegradation. In this study we conducted coring and static depth profile monitoring at the end of summer and end of winter for a layered sandy vadose zone in Perth, Western Australia. Transient on-line monitoring of vapors and O2 was also performed with in situ multilevel volatile organic compound (VOC) and O2 probes. For high soil moisture contents at the end of winter, vapors were shown to accumulate beneath a compacted, cemented layer approximately 0.3 m below the ground surface, and O2 penetrated only to depths of 0.4 m below ground. At the end of summer, when soil moisture was lower, O2 penetrated to depths of up to 1.5 m, and hydrocarbon vapors remained at or below this depth. Regardless of seasonal changes, sharp separations were seen between the depth of O2 penetration from the ground surface and the depth of penetration of the vapors upward from the hydrocarbon-contaminated zone. Modeling of steady-state O2 profiles indicated that a number of simple O2 consumption models might apply, including point-sink, distributed zero-order, or distributed first-order models, each leading to different biodegradation rates. Combining independent data with modeling helped determine the most appropriate model, and hence estimates of O2 consumption and hydrocarbon biodegradation. Also, reliable estimates of the biodegradation rate could only be calculated after consideration of the layered features.
Abbreviations: BTEX, benzene, toluene, ethylbenzene and the isomers of xylene • NAPL, nonaqueous phase liquid • NOM, natural organic matter • TOC, total organic C • TPH, total petroleum hydrocarbons • VOC, volatile organic compound
|
INTRODUCTION
|
|---|
HYDROCARBON VAPOR volatilizing to the ground surface from subsurface NAPL, such as residual gasoline, is a significant driver of health risk and the level of remediation required (Sanders and Stern, 1994; API, 1998). Ostendorf and Kampbell (1991) performed an initial field study of the behavior and biodegradation of petroleum vapors in soil impacted with petroleum NAPL. There were earlier studies by Karimi et al. (1987), who performed laboratory experiments, and Barber et al. (1990), who looked at a range of volatile contaminants in soil profiles near landfill leachate plumes. There have been several studies since looking specifically at gasoline-range vapors above petroleum impacted soil, and groundwater (Fischer et al., 1996; Laubacher et al., 1997; Davis et al., 1998b, 2000, 2001; Franzmann et al., 1999; Hers et al., 2000; Wang et al., 2003). Data from a number of studies were compiled by Roggemans et al. (2002). A range of modeling studies has also been performed (Jury et al., 1983; Baehr, 1987; Falta et al., 1989; Sleep and Sykes, 1989; Johnson and Ettinger, 1991; Lahvis and Baehr, 1996; Anderssen et al., 1997; Anderssen and Markey, 1997; Bekins et al., 1998; Johnson et al., 1999; Trefry et al., 2000, 2001; Turczynowicz and Robinson, 2001), some in combination with field studies (Barber and Davis, 1991a, 1991b; Öhman, 1999; Hers et al., 2000).
Despite the modeling efforts and related work, there are still only limited field data sets with sufficient detail for evaluating vapor processes in impacted soil profiles and for model validation. Additional well-documented studies are required (Johnson et al., 1999). Also, changes in soil moisture distribution and soil layering have been reported to impact vapor behavior and lead to complications when estimating biodegradation rates (Johnson and Perrott, 1991; Fischer et al., 1996; Davis et al., 2000).
We present the results of field research and modeling to quantify the role of a fine-scale moisture-retentive layer in a soil profile in changing the subsurface distribution of gasoline vapors and the major gases due to seasonal changes in moisture contents. Simple analytical and numerical modeling was performed to assess the impact of moisture variability on estimates of the biodegradation rate based on depth profiles. Coring, depth profiling of gases and vapors, and on-line monitoring of VOCs and O2 were used to determine changes from summer to winter. The moisture-retaining layer was found to have a dramatic effect on the distribution of vapors and major gases in the soil profile, and, if misinterpreted, calculated estimates of biodegradation rates.
|
MATERIALS AND METHODS
|
|---|
Field Site
The field site is located 50 km south of Perth, Western Australia. Leakage of gasoline has led to accumulation of light NAPLs in the zone of water table fluctuation and the volatilization of vapors into the overlying sandy vadose zone. The surficial soil profile and shallow aquifer described as the Safety Bay Sands (Playford and Low, 1972) comprise aeolian calcareous dune sands. The soil is predominantly leached sand with occasional limestone nodules and variably cemented layers, but with a pronounced shallow "cemented" limestone layer (described in more detail below). The climate is Mediterranean with an average annual rainfall of 771 mm and with typically hot dry summers (162 mm average rainfall from October to April) and mild wet winters (609 mm average rainfall from May to September). The water table is typically 2.6 to 3.3 m below ground. Soil coring and initial soil gas sampling indicated that, of the BTEX (benzene, toluene, ethylbenzene and the isomers of xylene) compounds, benzene was the least abundant with concentrations in the soil gas up to 200 µg L–1, and up to 25 mg kg–1 in gasoline NAPL in core material recovered from the site. Other gasoline compounds had concentrations that were typically at least an order of magnitude greater (see below). General characteristics of the site can be found in Johnston et al. (1998).
Sampling and Monitoring Installations
Multiport boreholes VMP1, VMP2, and VMP3 were installed in a west to east direction, 0.4 and 0.6 m apart, respectively, to allow sampling and measurement of soil gas depth profiles (see Fig. 1
for layout). The multiport installations consisted of bundles of 2-mm-i.d. and 3-mm-o.d. nylon tubes, with the ends covered with nylon mesh to avoid clogging by soil media. The nylon material has a low diffusion coefficient, which, together with adequate air purging reduces the potential for undue sorption of vapor compounds into the tubing. The depths of the monitoring ports on VMP1 were 0, 0.05, 0.15, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.25, 2.5, 2.75, and 3 m below ground surface. On VMP2 the ports were 0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, and 1.0 m, and on VMP3 the ports were 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, and 3 m. Samples were recovered primarily from selected ports of VMP2 and VMP3 on 14 Apr. and 1 Dec. 1999. Gas samples for organic analyses were collected in 20-mL glass syringes, and samples for major gases were collected in 50-mL plastic syringes after purging 100 mL of gas. Assuming radial symmetry, the extracted volume of 100 mL would be drawn from a radius in the soil of approximately 5 cm after taking porosity into account. The data from the uppermost sampling ports, at the 0- and 5-cm depth, has the potential to be affected by this sampling regime, but lower ports are at least 10 cm apart and so are unlikely to be affected. Analysis for the major gases was performed on site using gas chromatography with a thermal conductivity detector. Volatile organic compounds, including BTEX and total petroleum hydrocarbons (TPH), were analyzed using gas chromatography–mass spectrometry, after addition of a gaseous solution of fluorobenzene as an internal standard.
In December 1999 at three locations (VMP12, VMP15, and VMP16; Fig. 1), VOC and O2 probes were placed at the 0-, 0.5-, 1.0-, and 2.25-m depths below ground to provide near-continuous monitoring of total VOC vapor and O2 concentrations. See Patterson et al. (1999)( 2000) for details of construction, testing, and operation of the probes. The VMP12 and VMP15 locations were 7.5 and 2 m east of the midpoint between VMP1 and VMP2, respectively, and VMP16 was 5.3 m west of the midpoint between VMP1 and VMP2 (Fig. 1).
Coring and Physicochemical Parameters
Two soil cores were taken from within a 2-m radius of the VMP1-3 boreholes on 13 Apr. and 30 Nov. 1999 (i.e., the day before the principal dates of gas sampling). Their locations are noted on Fig. 1. Each core was recovered in two continuous sections over the soil profile by pushing an aluminium tube into the soil to a total depth of 3.4 to 3.5 m below ground surface. Soil bulk density, air-filled porosity, soil moisture, and NAPL contents were determined for 5-cm segments of the core. Core samples were analyzed for TPH, BTEX, 1,3,5-trimethylbenzene, naphthalene, and methylnaphthalenes. The soil was extracted using a diethyl ether/acetone solvent mixture. Analysis was performed using gas chromatography–mass spectrometry, with deuterated internal standards. Total carbon and total organic C (TOC) were also measured on core material. A subsample of the soil (typically 5–10 g) was dried until a constant weight was obtained, then crushed and homogenized using a ceramic mortar and pestle. The crushed sediment was divided into two duplicates. One of the duplicates was extracted using three successive portions (3 x 10 mL) of a hot solution of 3 M hydrochloric acid. The residue was well rinsed with deionized water then dried to a constant weight. The remaining duplicate sediment sample was untreated. Portions of a known mass from both subsamples were analyzed in an identical manner for C concentration using a LECO CS-444 C and S analyzer (LECO Corp., St. Joseph, MI). The C concentration of the extracted sediment was taken as the TOC, and the C concentration of the untreated sample was taken as total C, also comprising carbonate species.
Additional characterization of the site soils involved measurement of the hydraulic properties of the fine-scale (10–20 cm thick), shallow cemented limestone layer and the underlying sands. In the vicinity of the site, a small trench was dug to 0.25 to 0.30 m below the ground surface, exposing the top of the cemented layer. A portion of the trench was dug to 0.4 m to expose the underlying sands for testing. A ponded disc permeameter (Perroux and White, 1988) was used. The steel infiltration ring was placed on the cemented layer and sealed around the edge with bentonite slurry. A constant head of 0.014 m was established in the ring and the cumulative infiltration was measured through time. The unit was then transferred to the other side of the excavation onto the underlying sand where the cemented layer had been removed, and the procedure was repeated. The cumulative infiltration data were fitted to type curves to estimate the hydraulic conductivity (Perroux and White, 1988). Before and after the experiment, core material from the cemented limestone layer and sand portions of the trench was recovered to determine soil moisture.
|
RESULTS
|
|---|
Core Data and Physical Parameters
A cemented limestone layer of 0.10- to 0.20-m thickness appeared to cover much of the field site, with the upper surface at a variable depth below ground surface, typically 0.1 to 0.30 m below the ground surface. The topmost layer above the cemented zone comprised overlying sandy fill and assorted building rubble. Material below the limestone layer to the water table was typically fine to medium sands, with occasional thin coarse layers (Johnston et al. 1998).
Infiltration–permeameter tests showed that there was a consistent 20-fold contrast between the hydraulic conductivities of the cemented limestone and the sand beneath. At the latter time of measurement, the sand had a significantly larger moisture content (from 0.12 to 0.41), while the saturation state of the cemented limestone was essentially unchanged (going from 0.32 to 0.35). The permeameter testing indicated hydraulic conductivities of the order of 2 m d–1 (cemented limestone) and 35 m d–1 (sand), regardless of the measurement date.
For the two soil cores taken near VMP1, VMP2, and VMP3, volumetric fractions of soil (
s = "solid"), gas- or air-filled porosity (g = "air"), and liquid content (l = "liquid") and NAPL depth profiles are shown in Fig. 2 and 3
. Note that l consists of the volumetric soil moisture content (w) and the volumetric liquid NAPL fraction. In the soil profile above 2 m and below the base of the NAPL distribution, l = w; that is, there is no evidence of NAPL, although the soil moisture fraction may contain dissolved phase hydrocarbons. Note that uncertainties in the volumetric fractions are of the order of the second significant digit.
View larger version (33K):
[in this window]
[in a new window]
|
Fig. 2. Depth profiles of (left) soil volumes and (middle) NAPL contents determined from soil cores, and (right) soil gas measurements determined from VMP2 and VMP3 in April 1999.
|
|
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 3. Depth profiles of (left) soil volumes and (middle) NAPL contents determined from soil cores taken on 30 Nov. 1999, and (right) soil gas measurements determined from VMP2 and VMP3 on 1 Dec. 1999.
|
|
For the whole depth profile the soil core from April 1999 had a bulk density, 
b = 1.14 to 1.65 g cm–3 (with mean 1.38 g cm–3), and for the core recovered in November 1999, b = 0.95 to 1.61 g cm–3 (with mean 1.35 g cm–3). The TOC concentrations determined from soil subcores were low, typically below 0.1% by weight (0.04–0.07%), with one sample at the ground surface with a TOC of 0.13% (Table 1). Total C concentrations were 8.15 to 9.45% by weight. Temperatures were variable with depth and season, but typically were 20 to 25°C in summer and 15 to 20°C in winter across the bulk of the soil profile.
The depth profiles of s for April and November were similar. s was relatively uniform at approximately 52% by volume (or 0.52 m3 m–3) from 0.4 to between 1.5 and 2.0 m below ground, about 50% from between 2.25 and 2.5 to 3.5 m, and up to 60% around 2 to 2.25 m below ground. For the shallow depth interval (0–0.25 m below ground), s decreased from between 65 and 70% to about 45%. In April at the end of summer, w was typically 5 to 6% (or 0.05–0.06 m3 m–3) for the 0.3- to 1.0-m depth interval, whereas in November for this same interval w was 11 to 12%. For the shallow 0.1- to 0.3-m depth interval, w profiles showed a marked increase; in April the peak w was 20% while in November it had increased further to a peak of 25%. Concurrently, wherever w (or l) increased, g decreased from April to November in the same depth intervals. For example, in the shallow 0.1- to 0.3-m depth interval, g was about 30 to 35% in April but decreased to 25 to 32% in November, and in the bulk of the soil profile (0.3–0.4 to 2.25–2.75 m) g decreased from 33 to 43% in April to between 25 and 40% in November.
In April and November 1999, high NAPL concentrations were observed at similar depth intervals (2.2–3.25 m below ground) despite the water table being 3.3 m below ground in April (end of summer) and 2.6 m below ground in November (approximately the highest water level after winter rains). Note that in November 1999, a significant fraction of the NAPL mass, maybe up to 50% by weight, was effectively below the water table or in the capillary fringe, reducing the air–NAPL interfacial area and presumably volatilization. Benzene concentrations in the NAPL were low, relative to the xylenes and other compounds. Typically, the sum of the BTEX compounds made up <10% by weight of the NAPL TPH.
Soil Gas Data
We observed TPH vapor concentrations (Fig. 2, 3, and 4)
as high as 50000 to 70000 µg L–1 in the zone of NAPL contamination. Of the BTEX range of compounds, benzene and ethylbenzene (and 1,3,5 trimethylbenzene) had similar depth distributions, peaking at about 50 µg L–1, while toluene and m&p-xylene distributions were similar, peaking at about 1000 µg L–1. Interestingly, the TPH vapor concentration depth profile mimicked almost exactly the distribution of toluene adjusted by a factor of 50 (see Fig. 4). It is not possible to assess what vapor concentrations should be expected in equilibrium with the NAPL phase, since the mole fractions of each of the components in the NAPL are unknown. Summation of the concentrations of the BTEX range of compounds accounted for <5% by weight of the vapor phase TPH, based on the December 1999 data, compared with a weight ratio of nearly double that for the NAPL phase. This would imply that the NAPL consists of a larger mass fraction of compounds that are as volatile or more volatile than the BTEX compounds.
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 4. BTEX, 1,3,5 trimethylbenzene, and TPH hydrocarbon vapor concentration depth profiles for December 1999. Note that for these scales, the toluene data effectively coincide with the TPH data.
|
|
Distinct changes in the depth distributions of TPH vapor and major gases were observed between April and December 1999. In April (Fig. 2) with a low water table and a reduced w, O2 concentrations were non-zero from the ground surface to 1.5 m below ground, and TPH vapor concentrations were non-zero from the region of the NAPL to 1.25 m below ground. Figure 3 shows that, on 1 December, with a high water table and an increased w, O2 concentrations were non-zero to 0.4 m from the ground surface, while vapors penetrated from deep in the soil profile to within 0.4 m of the ground surface. The depth of overlap of the O2 concentrations and vapors was effectively the bottom of the depth interval of the zone of high soil moisture. It appears, therefore, that during high soil moisture periods the moisture-retentive layer at 0.1 to 0.3 m below ground was restricting vapor and O2 exchange. Note that CO2 concentrations for both dates were approximately 20 to 21% at depth, which is evidence of significant bioactivity, and most likely biodegradation of hydrocarbon vapors (e.g., McFarland and Sims, 1991; Franzmann et al., 1999). The December data also indicate that CO2 movement to the atmosphere may have been impeded by the higher moisture level in the shallow layer.
The in situ probe data at VMP12 for December 1999, as shown in Fig. 5
, were largely consistent with the 1 Dec. 1999 profile data in Fig. 3; that is, O2 concentrations were low at 0.5 m and zero at greater depths, and VOC concentrations were effectively zero at 0.5 m and non-zero and higher at greater depths. At later times and presumably as the soil profile dried, O2 concentrations increased even at the 1-m depth, and VOC concentrations at 1 m decreased to zero. In contrast, Fig. 6
shows that the probes at VMP16 reported non-zero O2 concentrations (0.01–0.125 atm or 1.0–12.5%) at the 0.5- and 1.0-m depths. Hydrocarbon vapors were absent at these depths. Hydrocarbon vapors were present at 2.25 m, but O2 was below detection levels. Hence, probe data from location VMP16 (for late December) are more consistent with depth profile data for April 1999 (Fig. 2) rather than for 1 December (Fig. 3). This seems to indicate that the shallow zone at the site is nonuniform and at VMP16 may retain less moisture or retain it for a shorter time than elsewhere across the site. This is shown by the intermittent effect of increased soil moisture on gas distributions in Fig. 6. Following 17 mm of rainfall on 14 and 15 Jan. 2000, O2 concentrations decreased sharply at the 0.5- and 1.0-m depths, but recovered to prerainfall conditions by late March 2000. It is surmised that the rainfall temporarily increased soil moisture contents in the shallow soil zone, decreasing the O2 flux from the atmosphere. Later in February and March (summer dry), drying of the soil profile decreased soil moisture contents to again allow increased O2 fluxes from the ground surface and atmosphere.
View larger version (27K):
[in this window]
[in a new window]
|
Fig. 5. Online monitoring of O2 and total VOCs at VMP12 from 24 Dec. 1999 to March 2000. A rainfall event occurred on 14 to 15 Jan. 2000.
|
|
View larger version (26K):
[in this window]
[in a new window]
|
Fig. 6. Online monitoring of O2 and total VOCs at VMP16 from 24 Dec. 1999 to March 2000. A rainfall event occurred on 14 to 15 Jan. 2000.
|
|
Note that during monitoring, O2 concentrations at VMP16 did not decrease to zero at the 0.5- or 1.0-m depths, and no change in the VOC concentrations was observed at these depths. This is consistent with the seasonal variability of the O2 and VOC vapor profiles, since from Fig. 2 and 3, wherever O2 was present vapors were apparently readily biodegraded; that is, there remained a relatively sharp separation between the presence of O2 and the presence of hydrocarbon vapors regardless of season.
GAS TRANSPORT MODELING
Modeling of such systems can be complex (e.g., Sleep and Sykes, 1989) or simple (e.g., Johnson and Ettinger, 1991). We chose to avoid complicated models of multiphase dynamics and multispecies reaction to focus solely on basic descriptions of the fate of a single gas species (O2) in the vadose zone in the presence of two potential lumped consumption mechanisms, biodegradation and mineralization. This is a common simplifying approach for assessing soil degradative capacity from field data. Even so, such simple transport models still engender considerable difficulties in interpreting consumption mechanisms and assigning degradation rates. Here simple analytical and numerical modeling is applied to the O2 data to aid interpretation in the presence of a shallow moisture retention layer and to estimate biodegradation rates. Steady-state depth profiles were assumed, so explicit seasonal, barometric, or temperature effects on the depth profiles are ignored. Jury et al. (1991) suggested that the effects of temperature, barometric pressure changes, and wind on vapor transport are minor (see also Massman and Farrier, 1992). Clearly, the depth profiles change seasonally and so are not strictly steady-state quantities. However, assuming an "active" soil profile depth of 2.75 m and an average diffusion coefficient for that depth of 2 x 10–6 m2 s–1 (see below), the characteristic time for vapor profiles to reach steady state after a step change disturbance would be <2 mo, which is less than the approximate 6-mo timeframe associated with seasonal transitions.
Assumptions made about the coupling between the O2 and TPH vapors and their degradation rate function drive the shape of the concentration depth profiles and the location of the zone or interface where biodegradation is occurring. Assuming a zero-order removal rate provides a quadratic (parabolic) form for the concentration depth profile, and a first-order rate function provides an exponential form, depending on details of the boundary conditions and assumed variability of the moisture content and diffusion coefficient. Ostendorf and Kampbell (1991) modeled coupled O2 and TPH removal assuming Monod kinetics (Monod, 1949) as the removal function. This provided a means of reproducing curvilinear O2 and hydrocarbon depth profiles and provided an overlapping zone of O2 and vapor concentrations. However, they did not model the individual fate of hydrocarbon compounds, such as benzene.
Ostendorf and Kampbell (1991) pointed out that if D is assumed to be uniform and constant, and if O2 and the hydrocarbon vapors are assumed to have the same (or similar) diffusion coefficients, then one can write a single equation for the net movement of a stoichiometrically summed surrogate X(z) given by
 | [1] |
where X = CH – CO/f, with CH as the hydrocarbon vapor concentration, CO the O2 concentration, and f the stoichiometric coefficient, and where z is depth. Note that Eq. [1] is independent of the functional form of the degradation rate. The boundary conditions are simply assumed to be held at constant concentrations, either zero for O2 at some depth L or for hydrocarbons at the ground surface (z = 0), and effectively held at C0 = 20.9% for O2 in air at z = 0 and the hydrocarbon concentration held at C1 at z = L. These conditions can generally be expressed as
 | [2] |
 | [3] |
where CH is the hydrocarbon concentration, CO is the O2 concentration, and L is the depth of the vadose zone. Note that C1 for hydrocarbons is assumed to be the highest hydrocarbon concentration observed at the base of the depth profile. With these boundary conditions, Eq. [1] can be solved to yield
 | [4] |
However, when X is calculated as a function of z/L for the April (and November) data, linearity, as indicated by Eq. [4], is not apparent. This nonlinear nature for the measured data is shown in Fig. 7
and implies that other transport or reaction processes or variable diffusion coefficients may be important in the vadose zone. Incorporation of discrete zones, spatial detail, or other reactive processes into modeling seems justified.
Hers et al. (2000) largely assumed coupling between the O2 and hydrocarbon concentration and degradation via zero- and first-order rates (not Monod), but also extended the modeling to two dimensions to accommodate house slab geometry. In this spirit, we now consider the following transient gas transport equation:
 | [5] |
where Ck(z,t) is the concentration of species k (kg m–3 of air-filled volume) at depth z (m), Deffk is the effective diffusion coefficient (m2 s–1) for species k, R(z) is a volumetric partitioning factor, and 
k is a lumped degradation term (kg m–3 soil s–1). At steady state we then have
 | [6] |
with the boundary conditions given by Eq. [2] and [3].
Diffusion Coefficient Formulation
It is clear from the governing Eq. [6] that if Deff(z) is constant with depth, then D and are linearly related; that is, if D changes, then changes similarly, given similar curvature of the concentration depth profiles. Refinement of the diffusion coefficient model should improve estimates of the degradation rate . Following an earlier soil characterization study (Öhman, 1999), we employ the following mixed soil tortuosity model:
 | [7] |
Equation [7] employs a Penman–Moldrup (Penman, 1940; Moldrup et al., 1996) tortuosity model for the gas phase and a Millington–Quirk tortuosity model (Millington and Quirk, 1961) for the fluid phase. Here the Henry coefficient, Kk, describes the partitioning of species k into the gas and fluid phases, and the Moldrup exponent was estimated from soil water potential characteristics of the local Safety Bay Sand via scaling relationships (Burdine, 1953; Campbell, 1974; Öhman, 1999). For O2 partitioning between air and soil water phases, the appropriate Henry constant is KO2 = 34. We neglect codiffusive (multi-species) corrections to the diffusion coefficient, since we will primarily be concerned with modeling the O2 dynamics in the zone where hydrocarbon vapors are absent.
Gross Features of the Effective Diffusion Coefficient
Equation [7] is a reasonable representation to assess the relative changes in Dsoil,effk resulting from moisture content changes between April and November 1999, and the likely impact on O2 fluxes and estimates of biodegradation. For the bulk of the soil profile (0.4–2 m below ground) for April (Fig. 2), average g = 41% and n = 47%; then from Eq. [7], Dsoil,effk/Dg,freek = 0.168. For November (Fig. 3), g = 35% and n = 48%; then from Eq. [7], Dsoil,effk/Dg,freek = 0.076. So for the bulk of the soil profile, Dsoil,effk declined by approximately 55% from April to November.
For the shallow zone of the soil profile (0.1–0.3 m below ground), for April, g = 33% and n = 53%; then from Eq. [7], Dsoil,effk/Dg,freek = 0.041. For November and December, g = 26% and n = 54%; then from Eq. [7], Dsoil,effk/Dg,freek = 0.013. For this narrow zone the diffusion coefficient for April is more than three times that for November.
Values of Dsoil,effk/Dg,freek for composite stratified media can be determined from the discrete harmonic mean:
 | [8] |
where L = L1 + L2 + ... + Lm. Most of the gas transport of interest occurs in the top 2 m of the profile, so in the current context Eq. [8] may be applied to the two-stratum (m = 2) system defined by the moist top 0.3-m zone and the underlying bulk zone to the 2-m depth. In this way, the effective two-stratum diffusion coefficient for April is calculated to be Dsoil,effk/Dg,freek = 0.127 and for November and December is Dsoil,effk/Dg,freek = 0.050. In essence then, the change in moisture content from April to November and December more than halves the effective diffusion coefficient for the upper 2 m of the soil profile.
Simplistically then, the degradation–consumption rate is directly proportional to changes in D, if steady state is assumed. Therefore degradation rates may decrease by two- to threefold from April to November and December based on moisture content changes and changes in D. However, this is only true if the curvatures of the O2 depth profiles are similar in April and December, which they are not, and that D is constant with depth. Regardless, it is clear that soil moisture has a strong influence on potential estimates of O2 consumption and biodegradation. This is explored more fully below, by modeling of moisture and curvature changes for the two dates.
Smoothed Soil Property Functions
Descriptions of the effective diffusion coefficient may be improved by explicitly incorporating information of the soil volumetric fractions into Eq. [7]. This was achieved by fitting each of the volumetric data sets of Fig. 2 (n, w, g 
n – w) with simple sums of Gaussian functions:
 | [9] |
to capture the essential features of the volumetric measures, and thereby yielding continuous (with depth z) effective functions for the volumetric quantities presented in Fig. 8
. In this way, Dsoil,effk was expressed as a rational function of Gaussian terms, yielding the curves in Fig. 9
. Gaussian parameters for these functions plotted in Fig. 9 are listed in Table 2. The presence of the moist soil layer around z = 0.2 m is clearly visible in the April and November curves as steep reductions in Dsoil,effk, which rebound to bulk values by z = 0.5 m. Further down the profile, Dsoil,effk slowly declines to liquid phase values near the water table at z = 3.5 m. Refined estimates of the throttling effect of the moist layer can be gained via mean diffusion coefficients 
defined over depth intervals by the continuous harmonic form
 | [10] |
For April 1999 and for the interval z 
[0, 0.3], 
/Dg,free=0.068, while for z[0.3, 1.75], /Dg,free=0.146, which is more than twice the diffusion coefficient for the upper moist layer. Thus, based on the measured soil properties and the above tortuosity model, and at steady state, we might expect the upper O2 profile to have twice the slope of the profile in the lower zone.
Modeling the Steady Oxygen Profiles
In the remainder of the paper, we restrict our attention to modeling and interpreting the O2 profiles to elucidate processes and understanding of temporal variations as found in our summer vs. winter experiment. This is attractive for assessment of vapor biodegradation at impacted sites since O2 is relatively inexpensive to analyze and simple to model compared with volatile hydrocarbons. Often, too, hydrocarbons vapors may be measured as TPH rather than speciated to give the behavior of individual compounds (e.g., Ostendorf and Kampbell, 1991). This makes linkage of O2 consumption to hydrocarbon degradation somewhat difficult. Even so, it is recognized that coupled modeling is desirable (Öhman, 1999) to provide constraints on biodegradation rates; indeed O2 consumption via natural organic matter (NOM) and reduced species in the vadose zone such as sulphides may occur (Johnston et al., 1998; Hers et al., 2000). Here we assume that NOM and other reduced species that may consume O2 over the depth profile did not change significantly in the year of monitoring (1999). This allows estimation of relative summer vs. winter O2 consumption rates and inference of hydrocarbon degradation rates.
In the following sections we apply several models to the steady O2 profiles to assess and understand the controlling mechanisms on the vapor transport. First, simple analytical models are compared with the steady O2 data to obtain gross estimates of lumped flux rates. Then a more complex model incorporating descriptions of spatial variability in the soil parameters is employed, refining the flux rate estimates and further elucidating other reactive processes. Finally, these models are extended to include distributed consumption terms such as may arise from intrinsic soil O2 demand.
Diffusive Descriptions of the Steady-State Oxygen Profiles
For a constant diffusion coefficient and assuming instantaneous O2 consumption occurs at a single depth in the soil profile (i.e., the O2 concentration is held at zero at that depth), steady-state solutions to Eq. [5] are linear with depth. Such a profile is clearly not representative of the O2 data measured in April 1999. However, in the light of the measured soil properties discussed above, a zonal or composite structure to the diffusion coefficient may be appropriate.
Steady Solutions for Piecewise-Constant D
Consider a soil column with a piecewise-constant (composite) effective diffusion coefficient given by
 | [11] |
where zc is the location in the column of the change in diffusion coefficient. The steady concentration profile is given by solving the homogeneous form of Eq. [6] (i.e., with zero consumption term, = 0). This may be done by solving Eq. [6] in each of the two zones where the diffusion coefficient is constant. Each zone solution has two constants of integration, one of which is fixed by application of the boundary conditions C1(0) = C0 for Zone 1, and C2(L) = 0 for Zone 2 for some depth L (based on the profile data). The remaining constants are fixed by requiring that the concentration solution is continuous at zc [i.e., C1(zc) = C2(zc)] and that the Fickian flux across zc is conserved; that is, D1dC1/dz = D2dC2/dz evaluated at zc. Solving these four conditions simultaneously yields the following steady composite solution:
 | [12] |
where 
= D1/D2. Note that the magnitudes of the diffusion coefficient components D1 and D2 do not appear in Eq. [12]—only their ratio is important in determining the steady concentration profile. Equation [12] may be fitted to the April measured O2 profile (C0 = 20%, v/v; L = 1.75 m) by seeking least-squares optimal values for and zc. The results are given in Table 3. This solution is plotted in Fig. 10a
and is an excellent fit to the April data. Taking this result in isolation, we may reasonably suspect that the April profile is diffusion controlled, with the top 30 cm of the soil profile having an effective O2 diffusion coefficient approximately one-fifth of the effective diffusion coefficient at depth. Recalling the earlier discussion on diffusion coefficient values from the tortuosity formulation section, we see that the approximate ratio of extremal diffusion coefficients for the upper moist and lower bulk layers is 0.041/0.168 = 0.24, in good agreement with the fitted value of 0.19. However, the earlier refined harmonic mean estimates of the diffusion coefficient, for the same vertical intervals specified by zc and L, yield a ratio of 0.068/0.146 = 0.47.
View this table:
[in this window]
[in a new window]
|
Table 3. Regression results for the different transport models. Asymptotic standard errors are supplied in parentheses.
|
|
View larger version (36K):
[in this window]
[in a new window]
|
Fig. 10. Oxygen steady-state profiles for April and December 1999, together with the fitted profiles from (a, b) the piecewise-constant D models of Eq. [12] and (c, d) variable D models of Eq. [13] and [14].
|
|
The same composite model may be fitted to the November data; see Fig. 10b and Table 3. In this case, the results are more equivocal, since by z = 0.4 m the measured O2 concentration declines to levels comparable to the measurement error. There is weak evidence for a change in slope of the profile between z = 0.2 and 0.3 m, based on a single data point. The location of the change in slope is consistent with the April results.
Steady Solutions for Variable D
Steady O2 profiles were calculated using the detailed diffusion coefficient defined by Eq. [7] and the April property functions of Table 2. The calculations were performed by noting that, in the absence of explicit consumption terms, steady solutions to Eq. [6] must satisfy Fick's Law:
 | [13] |
where the integration constant q is identified with the vertical O2 flux in the column. In the absence of an O2 consumption term, the flux q is constant throughout the column. In fact, q can be thought of as the volumetric O2 consumption or demand at the interface between the upper O2-rich zone and the lower hydrocarbon vapor zone. At a conceptual level, this implies that all the biodegradation would be occurring in a narrow reactive zone where O2 and hydrocarbons coexist. Integration of Eq. [13] generates a family of curves parameterized by q and by a boundary condition. For convenience, the atmospheric boundary condition is applied [i.e., C(0) = 20%, v/v]. The flux parameter q may be chosen to provide best least-squares fit to the measured profile.
Practically, the calculated C(z) profiles may match the measured data near the ground surface (q = 3.8 x 10–5% v/v m s–1), or the zero point near the 1.75-m depth (q = 2.7 x 10–5% v/v m s–1), but not both (Fig. 10c, dashed curves). In the assumed absence of a distributed consumption term, this may be rationalized in terms of a local O2 consumption process present at depth, such as might be expected from microbial degradation of hydrocarbon vapors. From Fig. 2 it is seen that the TPH concentration declines from high levels near the water table to zero at approximately z = 1.25 m. Accordingly, we may model the effect of microbial metabolism of TPH in this lower zone by an O2 consumption term that is zero higher in the column but non-zero at depth (i.e., at the overlap between the O2 and TPH profiles). In this case, Eq. [13] becomes
 | [14] |
The effect of this consumption-at-depth term would be to reduce the local O2 flux, and hence the absolute local O2 profile gradient, at depth, lifting the (fixed q) O2 profile to non-zero values. At the same time as these adjustments are made to the deeper part of the O2 profile, the upper profile (where = 0) remains essentially unchanged. Figure 10c (solid curve) shows such a result where an arbitrary Lorentzian functional form {i.e., 0/[1 + (z – z0)2/
2]} of amplitude 0, half-width = 5 cm, and centered at z0 = 1.25 m was assumed for . Least-squares estimated values for q and 0 are listed in Table 3.
In this way, we may ascribe the shape of the April O2 profile to two different effects operating at two different elevations in the column. First, in the column above z = 1 m, the O2 profile is diffusion-controlled, with the detailed shape determined by the spatial variation of the soil gas and moisture contents and porosity. Of particular note here is the role played by the upper moist soil layer. Second, at depths below 1 m, the O2 profile displays some characteristics of local consumption. This is inferred from the inability of constant-flux curves to match the O2 profiles over the full column. Matches are improved by the addition of localized consumption terms deep in the column, corresponding to the measured region of overlap between the O2 and TPH distributions.
We can try the same approach with the November O2 profile. Fitting Eq. [13] to the first four measured points yields q = 14.7 x 10–5% v/v m s–1, with the corresponding profile rapidly decreasing beyond z = 0.2 m (Fig. 10d). Adding a Lorentzian consumption term centered at z = 0.227 m allows the deeper part of the profile to be corrected. Unfortunately, the uncertainty in the O2 measurements obscures the nature of the profile beyond z = 0.25 m, so few quantitative conclusions can be drawn for the November profile.
Effects of Intrinsic Oxygen Consumption
Here we consider models incorporating variable diffusion coefficients, with either distributed zero-order or first-order O2 consumption throughout the column. These models apply to situations where the intrinsic soil biogeochemical consumption of O2 is significant, even in the local absence of petroleum hydrocarbon contaminants.
Variable D and Distributed Consumption
We incorporate into the general governing Eq. [5] detailed descriptions of the variation with depth of the effective diffusion coefficient and partitioning factor. The aim is to assess the influence of the upper moist soil zone on the gas transport properties of the vadose zone. In this case, we conform to the approach of Öhman (1999), who regarded the consumption of O2 as being related to the local gas fraction; thus, Eq. [5] becomes
 | [15] |
with Deff given by Eq. [7] and R(z) by g(z) + w(z)/K. At steady state, Eq. [15] becomes
 | [16] |
where
 | [17] |
In practice, was considered to consist of either zero- or first-order terms only. For zero-order consumption terms, Eq. [16] was integrated to the form:
 | [18] |
and solved using the same method as for the earlier local consumption model. For first-order terms, the transient form of Eq. [15] was integrated numerically to steady state. Then, for both consumption terms, the steady solutions were fitted to the O2 profiles by least-squares regression. Table 3 lists the regression results. Both zero- and first-order models gave good fits to the data, certainly within the bounds of analytical error, plotted in Fig. 11
.
View larger version (32K):
[in this window]
[in a new window]
|
Fig. 11. Observed steady O2 profiles for April and November 1999, together with the fitted profiles from variable D models incorporating distributed consumption.
|
|
In particular, the April data are well represented by both curves, but there were slight differences between the consumption models for the November data. At z = 0.4 m, where O2 supply is limited, the first-order model was better able to account for the curvature in concentration. However, this feature is small in comparison to analytical error, so it should not be used to favor one consumption term over the other. Also plotted in Fig. 11b and 11d are the November variable D profiles calculated assuming that the consumption rates are held fixed at the rates determined from the variable D regression for the April data (Fig. 11a and 11c, respectively). There are noticeable discrepancies between these hybrid curves and the regressed April profiles (Fig. 11a and 11c, respectively), indicating that the changes in soil properties, especially the moisture content, have greatly affected the O2 transport to the soil column. The hybrid curves themselves are reasonable approximations to the measured November O2 profiles; however, the agreement is improved, in a least-squares sense, by the November regressions on the consumption rate parameter 1. Unfortunately, due to the steepness of the November profile in the first 50 cm below the ground surface, adjusting the variable D profile by even a small amount as part of the optimization procedure corresponds to a large inferred change in consumption rate. This is the reason for the large standard errors in Table 3 and for the significant differences between the April and November consumption rates for the two consumption models. As a final comment on the regressions, we note that the diffusion coefficient formulations themselves are potential sources of uncertainty in the estimated rates. A gross sensitivity analysis was performed by varying the magnitudes of the diffusion coefficients used in the regression of Table 3. Relative changes of 50% in Deff led to changes in estimated fluxes and consumption rates of similar magnitudes (40–60%); that is, the response to scale perturbations in Deff is approximately linear. While this does not mean that the mixed diffusion model of Eq. [7] is necessarily appropriate for the soil type, it does suggest that other comparable diffusion models would yield similar flux and consumption rate estimates, at least to within an order of magnitude.
|
DISCUSSION
|
|---|
Literature Comparisons of Estimated Rates
A direct comparison of the zero-order rates can be made with laboratory-derived estimates from Franzmann et al. (1999), who took four soil samples from 0 to 2.8 m below ground at the site and performed spiked 14C-benzene mineralization experiments. They obtained zero-order estimates of 27 to 83 µmol kg–1 soil d–1. This converts to an O2 consumption rate of 2 to 6 x 10–5% s–1, which is comparable to the distributed consumption rate of 6.1 x 10–5% s–1 determined for April 1999 in our study (Table 3). The 6.1 x 10–5% s–1 O2 consumption rate can be converted to a hydrocarbon degradation rate of 6.6 mg hexane kg–1 soil d–1 or 6.0 mg hexane L–1 aqueous phase h–1, assuming a soil bulk density of 1.4 kg m–3 and a moisture filled porosity of 0.06 m3 m–3. Hers et al. (2000) estimated BTX zero-order mineralization rates of 0.6 to 1.4 mg BTX L–1 aqueous phase h–1. Fischer et al. (1996) estimated rates from field data of 0.5 to 40 mg kg–1 soil d–1, which is comparable to the other estimates here.
The first-order rates are more variable. The first-order benzene half-lives determined by Franzmann et al. (1999) were 11 to 72 d (1–7 x 10–7 s–1), which convert to an O2 consumption rate of 0.3 to 2.2 x 10–6 s–1 (accounting for stoichiometry). From data in Fischer et al. (1996), first-order rates of 1.4 to 3.6 x 10–6 s–1 (or for O2 = 0.45–1.1 x 10–5 s–1) can be estimated. Hers et al. (2000) gave BTX first-order rates of 0.5 to 1.2 h–1 (1.4–3.3 x 10–4 s–1), which convert to an O2 consumption rate of 0.42 to 1.1 x 10–3 s–1. Here first-order model fits gave estimates of O2 consumption of 0.09 to 2.6 x 10–4 s–1 (Table 3).
Diffusion vs. Consumption
Taken together, Fig. 10 and 11 permit direct comparison of the two transport approaches: diffusion + local consumption vs. diffusion + distributed consumption. Unfortunately, a priori it is not possible to quantify the precise processes that gave rise to the observed profiles, but we are in a position to compare the contributions to the profiles from the different models, and perhaps assess the relative likelihood of the applicability of the models in light of the independent data.
Absolute Consumption Rate Estimates
Care must be taken in assessing and comparing the parameter values listed in Table 3. The local consumption model yields estimates of q, the net O2 demand of the soil column. The distributed consumption models report values of the parameters, which represent the rates of O2 consumption in units of the partitioning coefficient (see Eq. [23]). Total O2 consumption rates integrated over the soil column, 
0, can be inferred from the 0 values by the integral
 | [19] |
where L is the column length of the O2 profile, taken to be 1.75 m for April and 0.4 m for November. 0 must have units of percent (v/v) per second, and the integral over the length of the soil column in Eq. [19] adds an additional length scale, ensuring that 0 is expressed in units of percent (v/v) meters per second. The April and November 0 values are calculated to be 3.7 x 10–5% v/v m s–1 and 5.3 x 10–5% v/v m s–1, respectively (Table 3). These values can be compared directly with the relevant q values from Table 3 of 4.5 x 10–5% v/v m s–1 and 21 x 10–5% v/v m s–1, respectively. We see that the soil O2 demand is somewhat less than the fitted influx, so there is potential for some O2 supply to deeper zones in the profile (e.g., for hydrocarbon degradation).
Similarly, the net O2 consumption rates for the soil column from the first-order model, 1, can be estimated from
 | [20] |
with the same column lengths as for the zero-order calculations. The April and November 1 values are calculated to be 4.1 x 10–5% v/v m s–1 and 32 x 10–5% v/v m s–1, respectively (Table 3). The entry flux rate (i.e., at z = 0 m) can be calculated by finite differences using Eq. [13], yielding q = 3.8 x 10–5% v/v m s–1 and 27 x 10–5% v/v m s–1, respectively. Note that no standard errors are reported for these finite difference estimates. In this first-order case, the influx rates are of the same order as the consumption rates, albeit slightly lower, meaning that the possibility of O2 supply to the deeper hydrocarbon-contaminated zones is excluded.
Interpretation of Consumption Rates
The above analyses show that using local and distributed consumption models leads to different conclusions. The local consumption results emphasize that the shape of the O2 profile is largely determined by the soil saturation state, and that inclusion of a localized consumption term at depth in the column, as may be inferred from consideration of the hydrocarbon data, can lead to acceptable matches with measured O2 profiles. This approach does not need to invoke any distributed, intrinsic O2 consumption process in the soil to match the data to within the analytical error. This is not to say, however, that distributed consumption is not occurring within the column, rather that after diffusion is taken into account, the contributions from any distributed consumption would likely be less than the error bounds on the data (this is discussed further below).
On the other hand, cosmetic matches with the data can be improved, sometimes substantially, by employing a distributed consumption model, especially a first-order model that is most appropriate for cases of limited O2 supply. For these models, the fitted consumption rates are such that the soil O2 demand essentially balances the O2 influx, providing little potential for O2 transport to the hydrocarbon-contaminated zone at lower depths in the column. This is counter to the intuitive model of hydrocarbon vapors being attenuated by microbial activity at the interface between the O2 and hydrocarbon distributions. Significantly, calculations showed that much of the change in O2 profile from April to November could be accounted for by simple changes in soil properties. This reinforces the idea that the O2 profile is essentially diffusion-controlled. In this sense, the distributed consumption terms may be regarded as parameters to be varied to fit the measured data points, with the variations being of the order of the analytical measurement error.
This de-emphasis on distributed O2 consumption is supported by the low site soil organic fraction and sulphide content, in contrast to the higher values found by Hers et al. (2000). Additional support for limited "background" O2 demand and consumption in the soil profile is provided by Davis et al. (1998a), who reported negligible O2 consumption following aeration of an uncontaminated vadose zone location nearby to the current site (within 800 m). Davis et al. (1998a) monitored consumption over only a 24-h period, so zero-order consumption on the basis of measurement errors was <1% v/v d–1 or 1.2 x 10–5% v/v s–1. This is <20% of the total O2 consumption inferred from modeling the O2 profiles assuming variable D and distributed zero-order consumption (Table 3), which gives greater confidence that secondary consumption of O2 by NOM or other reduced species, such as sulphides, is limited in this instance.
Another assessment to support this de-emphasis on distributed O2 consumption is to infer O2 consumption rates from hydrocarbon disappearance (or biodegradation) estimated from the vertical depth profiles of vapors in Fig. 2 and 3 and compare these with the O2 consumption estimates generated directly from the O2 depth profiles, as above. The net flux of hydrocarbon vapor from the residual NAPL zone at depth toward the ground surface, qh, can be written as
 | [21] |
For the April data, the hydrocarbon concentration decrease was largely linear with a gradient of approximately 65000 µg L–1 for a 1-m vertical depth interval of 1.25 to 2.25 m. The diffusion coefficients for the BTEX and naphthalene range of hydrocarbons in air have been estimated as 7 to 9 x 10–6 m2 s–1 (Lugg, 1968; Karimi et al., 1987; Grathwohl 1998). From Eq. [7] and the likely range of Dsoil,effk/Dg,freek 
, the effective hydrocarbon vapor diffusion coefficient estimates in the sandy vadose zone at the site lie in the range 0.7 to 2.3 x 10–6 m2 s–1. Substituting into Eq. [21] gives qh = 0.05 to 0.15 µg vapor L–1 air m s–1. Assuming 3.2 g of O2 reacts with each 1 g of hydrocarbon vapor (usually hexane is taken as a model hydrocarbon), and converting molar volumes from Standard Temperature and Pressure (STP) to 25°C, allows us to convert qh to qo = 1.3 to 3.8 x 10–5% O2 m s–1, which is the apparent O2 flux or consumption rate at depth 1.25 m.
We can take the same simple approach for the O2 profile for April 1999 over the depth interval from 0.3 to 1.5 m. Approximating Dsoil,effk/Dg,freek , assuming the diffusion coefficient for O2 in air is approximately 2 x 10–5 m2 s–1, and multiplying by the O2 gradient (10% to the 1.2-m depth), yields qo = 1.7 to 4.2 x 10–5% O2 m s–1. The O2 flux (and implied consumption at 1.25–1.5 m depth) calculated from the hydrocarbon distribution and the O2 distribution are a close match. The consumption rates are also comparable to those generated from simulation of variable D with local consumption (Table 3). Although simplistic, this analysis adds considerable weight to the argument that most if not all of the O2 is being consumed in the soil profile by aerobic biodegradation of hydrocarbon vapors within the narrow zone of contact between the O2 and hydrocarbon vapors, and only limited O2 consumption is likely to be occurring in the bulk of the soil profile due to processes other than hydrocarbon vapor biodegradation.
Of course, this conclusion is not easily confirmed by November depth profile data. As discussed earlier, hybrid calculations using November soil property functions and April consumption rates were able to account for most of the changes in soil O2 profile for the April to November period. The large magnifications (with high standard errors) in consumption rates obtained for November resulted mainly from attempting to fit the steep November data. Other factors may also add to the uncertainty. For example, soil heterogeneity may play a role, since the soil core and VMP2 and VMP3 data are from different locations, albeit within a 2-m radius. Soil temperature changes, especially in the shallow soil profile for November, may also influence biodegradation rate estimates.
|
CONCLUSIONS
|
|---|
In a sandy vadose zone contaminated by gasoline NAPL, a shallow moisture retentive layer has seemingly had a dramatic effect on the depth distribution of hydrocarbon v
TOP
ABSTRACT
INTRODUCTION
