A comparison of three methods to assess natural source zone depletion at paved fuel retail sites

Natural source zone depletion (NSZD) encompasses all processes that result in petroleum hydrocarbon light non-aqueous phase liquid (LNAPL) mass loss. Vertical gas transport between the subsurface and atmosphere is a key component of NSZD. Gas exchange with the atmosphere may be restricted at sites with ground cover, which is typical for European fuel retail sites. This raises questions of whether, and to what extent, the generic NSZD conceptual model applies at these sites. Here, we present a study that evaluated how concrete and asphalt pavement affected NSZD processes and data interpretation for three NSZD assessment methods: soil gas concentration gradient, biogenic heat and carbon dioxide traps. All methods demonstrated that NSZD was occurring and NSZD rates were generally within the low end of values reported in the literature for unpaved sites. However, there was considerable variability in the rates, which highlights the need for careful examination of the conceptual site model and potential interferences for each method. The results demonstrate the viability of soil gas and temperature data collected from existing monitoring wells screened into the unsaturated zone without the need for additional, intrusive subsurface installations. The results also provide useful guidance for developing optimal long-term NSZD monitoring approaches, where necessary.

Light non-aqueous phase liquids (LNAPL) are immiscible organic liquids that are less dense than water. LNAPL-forming substances include petrol (gasoline), diesel, heating oils and jet fuel (kerosene). The occurrence of LNAPL in the subsurface can be the result of various kinds of releases at locations where these products are manufactured, stored or sold. When a release occurs, LNAPL will percolate downward under the influence of gravity and may spread laterally owing to geological heterogeneity or the presence of other preferential migration pathways. If a sufficient volume is released, LNAPL will continue to migrate downward into the saturated zone where it can spread laterally, often forming LNAPL bodies that are partially above and below the water table, similar to an iceberg . As LNAPL spreads, an increasing fraction of the LNAPL is trapped as a discontinuous non-wetting phase by capillary forces (i.e. residual LNAPL). Thus, an increasing volume of the released LNAPL is present as an immobile, residual phase, resulting in an overall decreasing volume of mobile LNAPL distributed over a larger volume of the subsurface.
Following a release, the vertical and lateral extent of LNAPL typically reach a stable condition on a timescale of weeks to years, depending on a number of parameters, including the release history, aquifer matrix characteristics, LNAPL physical properties and the rate at which LNAPL are depleted through natural and/or engineered processes (ASTM International 2014;CL:AIRE 2014;ITRC 2018). The distribution of LNAPL in the subsurface following a release is referred to as the 'source zone'. The source zone comprises both residual and potentially mobile LNAPL that can act as a source of contamination for groundwater or soil gas (ITRC 2018).
Natural source zone depletion (NSZD) encompasses all attenuation processes that result in LNAPL mass loss in the subsurface (Garg et al. 2017;ITRC 2018). These processes include physical mass transfer by dissolution and vaporization of chemical constituents to the aqueous (groundwater) and gaseous (soil gas) phases and biodegradation of LNAPL constituents through microbial-facilitated reactions.
The efficacy of natural attenuation of petroleum hydrocarbons in groundwater has been well established since the early 1990s (NRC 1993;Rice et al. 1995). Although there has long been evidence that microbiological degradation processes responsible for natural attenuation in dissolved phase plumes were also occurring within LNAPL source zones to 'weather' or change the composition of LNAPL (e.g. Christensen and Larsen 1993), there was a common historical perception that biodegradation of the source material itself was limited (Lyman et al. 1992;Newell et al. 1995). More recent research on NSZD at petroleum LNAPL sites (e.g. Johnson et al. 2006;Garg et al. 2017;CRC CARE 2020a) has demonstrated that the rate of natural LNAPL depletion is often of the order of thousands to tens of thousands of litres of LNAPL per hectare per year (l ha −1 a −1 ). The observation of natural depletion rates of this magnitude has highlighted the significance of NSZD in LNAPL conceptual site model (LCSM) development (e.g. Mahler et al. 2012;Lundy 2014) and site management decision making (CL: AIRE 2014). NSZD measurements are more frequently collected to better understand the relative benefit of active LNAPL remediation alternatives at LNAPL-contaminated sites, and can provide an alternative or supplement to approaches such as skimming or bailing of LNAPL (ITRC 2018;CL:AIRE 2019;CRC CARE 2020b).
Several institutions have published guidance and information regarding NSZD in different parts of the world (ITRC 2009;API 2017;CRC CARE 2018;CL:AIRE 2019); however, evaluation of NSZD is still limited outside the USA and Australia where the bulk of research on NSZD has been performed. The limited uptake may be due to a combination of different factors and gaps in our understanding of NSZD processes (see, e.g. CL:AIRE 2019). The aforementioned studies reported in the literature have largely focused on sites without any surface sealing. Effective gas exchange with the atmosphere may be limited at paved sites, which is typical for most service stations in Europe. Several recent guidance documents have highlighted the potential for measurement errors and bias using carbon dioxide (CO 2 ) efflux methods when lowpermeability surface cover materials are present (API 2017;CRC CARE 2018;ITRC 2018). When deployed at sites with lowpermeability artificial surfaces (e.g. concrete or asphalt), these methods typically involve measurements of CO 2 efflux from a hole drilled through the concrete or asphalt. Penetration through the concrete or asphalt creates a preferential pathway for soil gas flux, resulting in high-biased efflux results.
Additionally, the presence of low-permeability pavement materials can inhibit vertical gas migration and affect soil gas concentration profiles that are typically relied on for NSZD calculations using the soil gas concentration gradient method. For example, Coffin et al. (2008) found evidence of restricted soil gas exchange with the atmosphere beneath an asphalt surface at a fuelcontaminated site, resulting in build-up of methane and CO 2 in soil gas. However, Roggemans et al. (2001) compared soil gas, oxygen (O 2 ) and hydrocarbon profiles at 15 sites where paved conditions were present, and found that pavement did not necessarily inhibit vertical gas exchange between the atmosphere and the subsurface.
The influence of cracks and other preferential flow paths in pavement materials and the degree to which lateral transport of O 2 affects NSZD processes and data interpretation represent key unknowns at paved sites. Although the current body of literature and guidance on NSZD provides valuable insight into the effects of pavement on specific NSZD measurement techniques, there are very few field-based research papers that have directly addressed the issue using multiple data collection methods.
This work aims to help close the gap in understanding the effects of paved surfaces on NSZD measurements and data interpretation. The following three methods were applied to estimate NSZD rates at a former petroleum retail site to assess whether paved conditions inhibit NSZD, and to compare results of different NSZD measurement methods under similar conditions: • the soil gas concentration gradient method, which uses measurements of subsurface O 2 and CO 2 concentration profiles (these are the key chemical species involved in aerobic degradation of hydrocarbon and methane oxidation) to estimate NSZD rates (Johnson et al. 2006;Davis et al. 2009;Ririe and Sweeney 2018;Kulkarni et al. 2020); • the biogenic heat method, which uses measurements of subsurface temperature across the source zone to estimate heat flow gradients that arise from exothermic hydrocarbon biodegradation reactions (Mohr and Merz 1995;Subramanian et al. 2011;Ririe et al. 2013;Sweeney and Ririe 2014;Warren and Bekins 2015;Askarani et al. 2018;Kulkarni et al. 2020); • the CO 2 efflux method (CO 2 trap method), which measures CO 2 efflux across the ground-atmosphere interface (McCoy et al. 2014;Tracy 2015;Sihota et al. 2018;Kulkarni et al. 2020). The CO 2 traps used for this study were designed specifically for application at sites with impervious ground cover conditions (E-Flux 2017).
Multiple measurements were made using each of these NSZD methods over 1 year to evaluate seasonal influences on NSZD rates. The assessment included measurements in background (unaffected) areas as well as in gasoline-and diesel-affected areas. All of the measurements, including measurements at the background location, were made in areas covered by concrete and asphalt pavement. The use of multiple measurement methods provided insights into the conceptual site model and allowed for identification of site-specific interferences for some of the measurements.

Field site description
The study site, located in Mediterranean Europe, was formerly occupied by a retail fuelling station built in 1966 that stored fuel in nine 20 000 l underground storage tanks (USTs). Above-ground equipment and buildings were removed as part of decommissioning in 2009. The underground infrastructure was left in place, emptied of liquid fuels and filled with solid foam ( piping) and grout (tanks). A site layout showing the former retail station features is shown in Figure 1. Additional detail on site infrastructure has been given by Concawe (2020).
The entire surface of the study site was paved. Concrete covered most of the central portion of the site, and asphalt covered the surrounding areas (Fig. 1). Based on soil boring logs, the concrete slab in the central portion of the site is c. 20-30 cm thick, and the asphalt varies from c. 30 cm thick to the west (near the background locations SV5, S-21 and C5) to 10 cm thick to the east, where monitoring wells T1A-C and T2A-C are located (Fig. 2). The pavement was generally in a poor state of repair, and therefore may not reflect engineering conditions at a modern fuel retail site. Cracks were present in all of the pavement materials. Although no attempt was made to quantify specific characteristics of the cracks (e.g. aperture or density), it was noted that larger, fully penetrating cracks were generally more abundant in the concrete, suggesting that there may have been greater potential for gas exchange between the subsurface and atmosphere in these areas. Below the pavement and surface fill, native unconsolidated material of pebbles and gravels in a fine sand to clayey matrix extends down to the water table at 8-12 m below ground surface (m b.g.s.) with groundwater flow to the SE (Fig. 2). Additional detail on subsurface characterization data has been given by Concawe (2020).
Sampling of soils for analysis of total petroleum hydrocarbons (TPH) identified the presence of shallow, residual LNAPL impacts from grade to 3 m b.g.s. near former USTs and dispensers (Fig. 1). Additionally, residual gasoline LNAPL impacts were present on the western and central portion of the site, comprising a dissolved phase plume of lighter fraction petroleum hydrocarbons (C 6 -C 10 -range TPH), benzene, toluene, ethylbenzene and total xylenes (BTEX), and methyl tert-butyl ether/ethyl tert-butyl ether (MTBE/ETBE). Diesel LNAPL was present on the eastern portion of the site with gauged LNAPL thicknesses up to c. 50 cm, and heavier TPH fractions (C 10 -C 40 ) present in groundwater. Polycyclic aromatic hydrocarbons (PAH), primarily naphthalene, were found in both affected areas. The extents of the LNAPL and dissolved phase petroleum hydrocarbon plumes are delineated, and historical fluid level gauging and dissolved phase monitoring data indicate that LNAPL and associated dissolved phase plumes are stable. An area that was identified to have no hydrocarbon impacts was established on the west side of the site (identified as Background in Fig. 1), which is hydraulically up-gradient of the hydrocarbon-affected zones. Additional detail on subsurface geological conditions and distribution of LNAPL impacts at the study site have been given by Concawe (2020).

Relevant stoichiometry and conversion factors for NSZD rate estimation
The methods used for quantifying NSZD rates in this study rely on mass-and/or energy-balance approaches. NSZD rates are not determined by directly measuring changes to the source zone. Rather, rates are inferred based on stoichiometric relationships between products and reactants involved in NSZD processes, specifically, by measuring the flux of electron acceptors (e.g. O 2 ) into the source zone, or by measuring the flux of petroleum degradation products such as CO 2 or excess heat from biodegradation out of the source zone. Following Johnson et al. (2006) and Garg et al. (2017), the dominant hydrocarbon biodegradation reactions for this study were assumed to be methanogenesis (equation (1)) followed by aerobic oxidation of methane (equation (2)), or direct aerobic oxidation (equation (3)). In equations (1)-(3) a and b represent the number of carbon and hydrogen atoms in a given hydrocarbon compound, respectively: These equations provide a basis for estimating NSZD rates by measuring the flux of O 2 into a source zone (soil gas concentration gradient method) and/or measuring the flux of petroleum degradation products such as CO 2 (soil gas concentration gradient method and CO 2 trap method) out of the source zone. Additionally, the change in enthalpy or heat of reaction can be calculated from the internal energy of the products and reactants in equations (1)-(3). The heat of reaction is then used to convert measurements of subsurface heat flux associated with petroleum biodegradation into NSZD rates.
Whether degradation occurs through methanogenic degradation followed by methane oxidation (equation (1) followed by equation (2)) or through direct aerobic oxidation (equation (3)), the reactants and products are ultimately the same. The stoichiometric relationships can be expressed in terms of the mass of hydrocarbon degraded per unit mass of O 2 consumed, the mass of hydrocarbon degraded per unit mass of CO 2 produced, or the amount of energy released per unit mass of hydrocarbon degraded. For this study, octane (C 8 H 18 ) was used as a representative hydrocarbon for gasoline and hexadecane (C 16 H 34 ) was used as a representative hydrocarbon for diesel fuels (Bacha et al. 1998). The resulting stoichiometric coefficients for O 2 utilization (S O 2 ; 0.29 g C 8 H 18 per g O 2 and 0.29 g C 16 H 34 per g O 2 ), CO 2 production (S CO 2 ; 0.32 g C 8 H 18 per g CO 2 and 0.32 g C 16 H 34 per g CO 2 ), or heat released to the formation (DH rxn ; −47.9 kJ g −1 C 8 H 18 and −47.2 kJ g −1 C 16 H 34 ) were calculated using molecular weights and standard heat of formation for each of the compounds represented in equations (1)-(3) presented in Table 1 (Haynes 2012). The stoichiometric coefficients for O 2 utilization, CO 2 production and the heat released from biodegradation are relatively invariant for a broad range of hydrocarbons on a mass basis. Thus, although C 8 H 18 and C 16 H 34 were used to represent the gasoline and diesel LNAPL, respectively, for this study, the use of alternative representative hydrocarbons would not significantly affect the results for any of the three methods utilized.
Equivalent volumetric LNAPL depletion rates were calculated using LNAPL density (r n ) values of 0.77 g ml −1 for gasoline area locations, and 0.85 g ml −1 for diesel area locations.

Soil gas concentration gradient method
As indicated in equations (1)-(3), biological depletion of LNAPL consumes O 2 and produces CO 2 , creating a chemical gradient that drives diffusive gas flux in the direction of lower concentrations. Assuming gas transport in the subsurface can be adequately modelled as a steady-state, 1D (vertical) diffusion process, the mass flux of O 2 into the subsurface or transport of CO 2 from an LNAPL source at depth toward ground surface can be estimated using Fick's first law (Johnson et al. 2006;Davis et al. 2009;API 2017): where J i is the mass flux of O 2 or CO 2 (g m −2 s −1 ), DC i =Dz is the soil gas concentration gradient (g m −4 ), and D eff i is the effective gas diffusion coefficient (m 2 s −1 ). Soil gas composition profile data collected from nested soil gas probes (Johnson et al. 2006) and monitoring wells screened into the unsaturated zone (Sookhak Lari et al. 2017; Sweeney and Ririe 2017) provide a direct measurement of vertical concentration gradients, and effective gas diffusion coefficients (gas diffusivity values) were estimated using the Millington-Quirk expression (Millington and Quirk 1961) as a function of soil total porosity (θ T ) and soil gas saturations (S g ): Average values for soil porosity (0.18 cm 3 void per cm 3 soil) and gas saturation (0.58 cm 3 gas per cm 3 void) were determined from laboratory analysis of site soils (Concawe 2020) and D air i , the molecular diffusion coefficient for O 2 or CO 2 in air, was defined as 2 × 10 −5 or 1.6 × 10 −5 m 2 s −1 , respectively (Hillel 1998).
O 2 and CO 2 mass flux values from the soil gas concentration gradient method were converted into equivalent hydrocarbon depletion rates (R GM NSZD ) in units of litres of LNAPL per hectare per year (l ha −1 a −1 ) using the stoichiometric coefficients and density values described above, using equation (6): For this study, soil gas composition profiles were measured during five quarterly monitoring events between June 2017 and July 2018. During each event, soil gas composition was measured at five multilevel soil vapour probes (SVPs) installed in the source zones (SV1 and SV2 in the diesel area, SV3 and SV4 in the gasoline area) and background area (SV5) (Fig. 1). Multi-level SVPs were constructed using polytetrafluoroethylene (PTFE) tubing with 6.4 mmouter diameter and 4.8 mm inner diameter, and stainless steel AMS probe tips packed in sand and separated from one another by layers of bentonite, in accordance with ITRC (2014) guidance. Each SVP included three vertical sampling intervals for assessing changes in soil gas composition with depth, including one probe installed in the upper portion of the unsaturated zone, beneath the pavement (1.0-1.5 m b.g.s.), one probe installed at an intermediate depth (4.8-5.2 m b.g.s.), and one probe installed near the base of the unsaturated zone, c. 1-2 m above the water table (6.8-10.0 m b.g.s.). A 15-20 cm thick sand filter pack was installed to allow unrestricted gas flow to the probes, and a bentonite seal was installed in the annular space between each of the soil gas probe intervals. Additional details on SVP construction have been given by Concawe (2020). Soil gas samples were analysed for O 2 , CO 2 , methane (CH 4 ) and hydrogen sulphide (H 2 S) with an MSA Safety Services Altair 5X multigas meter and volatile organic compounds (VOC) with a handheld miniRAE 3000 photoionization detector (PID). All instruments used for field soil gas measurements were calibrated in accordance with manufacturer specifications prior to use. Additionally, soil gas profiling for O 2 and CO 2 was completed within several monitoring wells across multiple depths starting at the top of the well screen and continuing down to the base of the unsaturated zone in 1 m intervals. Soil gas profiling in wells was completed using low-volume purge methods described by Sweeney and Ririe (2017). The foundation for using measurements of increased temperature in the subsurface attributable to biodegradation to determine hydrocarbon depletion rates was established decades ago (Mohr and Merz 1995), and has more recently become a tool for monitoring remediation performance (Subramanian et al. 2011) and NSZD (Ririe et al. 2013;Sweeney and Ririe 2014;Warren and Bekins 2015;Askarani et al. 2018). The biologically mediated NSZD processes that destroy hydrocarbons and alter the composition of soil gas also release heat and create temperature gradients away from zones where microbial degradation is occurring. Depending on soil thermal properties and climate, seasonal atmospheric temperature changes also modify soil temperatures to depths of 8-25 m (Mount and Paetzold 2002;Busby 2015). At sites with LNAPL impacts at depths that are affected by seasonal temperature changes or other background thermal influences, subsurface temperatures can be conceptualized as a combination of heat from LNAPL depletion and background heat transport processes. Calculation of NSZD rates from subsurface temperature measurements requires an understanding of the background temperature distribution to identify temperature increases attributable to NSZD. The increase in temperature attributable to NSZD at a given depth and time of year, DT NSZD (z, t), is calculated using equation (7) (Sweeney and Ririe 2014;Askarani et al. 2018): where T SZ (z, t) and T BKGD (z, t) are the temperatures at depth z and time t within the LNAPL source zone and at a non-affected location, respectively. For this study, temperature data from monitoring well S-21 (identified as Background in Fig. 1) located in an area where LNAPL impacts were absent, was collected to measure background subsurface temperature fluctuations to be subtracted from temperatures recorded in gasoline-and diesel-affected areas.
Using the background-corrected temperatures calculated from equation (7), the sum of the upward and downward conductive heat flux (q T ) away from the depth at which the maximum backgroundcorrected temperature difference is observed can be calculated using Fourier's law of heat conduction (equation (8)): where K u and K d are the effective thermal conductivity values (W m −1 K −1 ) for subsurface materials above and below the depth of the maximum background-corrected temperature, respectively, and (DT =Dz) u and (DT =Dz) d are the upward and downward temperature gradients. Vertical temperature profile data were collected every 3 months from October 2017 to July 2018 in 14 wells located in LNAPLaffected areas and one background well location (Fig. 1). Snapshot temperature profiles were manually measured using a HOBO S-TMB-M017 temperature sensor with accuracy of ±0.3°C and HOBO U14-002 LCD data logger. The temperature sensors were lowered in the site monitoring wells, and temperature measurements were recorded in 1 m increments from the surface after allowing the temperature of the thermistor to equilibrate for up to 3 min at each depth interval in the unsaturated zone, and 1 min for each depth interval in the saturated zone, using methods provided by Subramanian et al. (2011), Ririe et al. (2013) and Sweeney and Ririe (2014). Continuous subsurface temperature data were collected using HOBO TidbiT® MX2203 and Thermochron iButton® DS1922L temperature data loggers with an accuracy of ±0.5°C set to 1 h intervals installed in background well S21, gasoline area well S5 and diesel area well S7 (Figs 1 and 2). A total of 36 temperature data loggers (including 12 in the background well S21, 11 in the gasoline impacted well S5, and 13 in the diesel-affected well S7) were installed at 1 m intervals from the surface to the bottom of the wells. An additional data logger was installed above ground with no direct exposure to sunlight or wind to record the ambient temperature at the site.
Conductive heat flux from biodegradation was estimated for source zone locations using measured, background-corrected temperature gradients and thermal conductivity values derived through analysis of continuous temperature profile data collected from the background location. The method for determining representative thermal conductivity values involved estimates of thermal diffusivity (the ratio of thermal conductivity to volumetric heat capacity) determined using a combination of the phase lag and amplitude ratio methods described by Carson (1963) and de Jong van Lier and Durigon (2013). An average thermal conductivity of 1.6 W m −1 K −1 was determined from site-specific estimates of thermal diffusivity and volumetric heat capacity (Concawe 2020).
For this study, the maximum background-corrected temperatures were observed at a depth of 10 m b.g.s., which corresponded to the total depth of the background well (S-21). Thus, only upward heat flux could be estimated (the downward heat flux terms, K d (DT =Dz) d , in equation (8) were not included), providing a minimum estimate of the total heat flux attributable to biodegradation.
Heat flux estimates were then used to estimate equivalent LNAPL depletion rates based on stoichiometric relationships (equations (1)-(3)) and the heat of reaction (Table 1) using equation (9): CO 2 trap method CO 2 flux from the subsurface to atmosphere was measured using CO 2 traps installed in autumn (October 2017) and spring (April 2018) to estimate the biodegradation rates in multiple seasons. Four traps were placed in the source areas (C1 and C2 in diesel area and C3 and C4 in gasoline area) and one trap (C5) was placed in the background area (Fig. 1). CO 2 traps were constructed by E-Flux, LLC (Fort Collins, CO, USA). They comprised a 10 cm diameter polyvinyl chloride (PVC) receiver pipe that is inserted into the soil (in this case, through the asphalt or concrete pavement materials) and a main trap body that contains two layers of CO 2 sorbent consisting of a mixture of calcium and sodium hydroxides. The receiver pipe provides a defined measurement area through which all upward soil gas flow was directed to the trap. The upper layer of sorbent material in the trap captures ambient CO 2 entering from the atmosphere, and the lower layer captures CO 2 entering the trap from the subsurface. Additional details on CO 2 trap design and deployment for measurements of biodegradation have been documented by McCoy et al. (2014).
The conventional CO 2 trap design, as described by McCoy et al. (2014), does not account for the presence of low-permeability pavement materials at the ground surface. Pavement materials (e.g. asphalt and concrete) have much lower permeability and gas diffusivity than the CO 2 trap sorbent material and most noncohesive soils. When conventional CO 2 trap receivers are installed through pavement, this can create a preferential pathway for gas flow (i.e. a chimney or stack effect) that magnifies the CO 2 flux from the subsurface into the trap, leading to an overestimation of the biodegradation rates (API 2017). To address this potential bias at the study site, CO 2 traps were installed inside a sealed vault and the tops of traps were connected to the soil under the pavement with a vapour pin to equalize pressure conditions (Fig. 3). This specific method developed by the CO 2 trap supplier (E-FLUX) for paved sites was developed with the goal of creating a closed system to ensure that the traps were exposed to natural gas flux and did not become a preferential pathway for flow of gas into or out of the subsurface. The method was designed with the assumption that there is a sufficient number of cracks in the pavement material that gas transport in the unsaturated zone is primarily 1D (vertical), except at shallow depths directly beneath the pavement where gas flux is concentrated within cracks (E-Flux 2017). The method also does not account for the presence of LNAPL impacts directly beneath the pavement.
Sorbent recovered from the traps was dried, homogenized and analysed for carbonate and fossil fuel-derived carbon content by E-Flux. The 14 C ratio in the captured carbon dioxide ( 14 CO 2 ) from the traps was analysed by the laboratory to discriminate between CO 2 resulting from degradation of petroleum hydrocarbons and CO 2 resulting from degradation of naturally occurring organic matter. The 14 CO 2 content evolved from biological degradation of petroleum hydrocarbons will be negligible (non-detect). By contrast, 14 CO 2 resulting from degradation of naturally occurring organic matter in near-surface environments where organic matter has more recently been in equilibrium with the atmosphere will reflect the mean age of the organic material being degraded.
The 14 CO 2 data were used to correct the NSZD rates using a two end-member mixing model, assuming that all CO 2 captured by the traps represents a combination of petroleum hydrocarbon degradation and degradation of organic matter that has recently been in equilibrium with the atmosphere. The assumption that all carbon sources other than petroleum hydrocarbons contain 14 C concentrations near modern atmospheric levels will over-predict the percentage of petroleum-derived CO 2 and NSZD rate for sites with buried organic materials of sufficient age to have a lost a significant fraction of 14 C isotopes (Lundegard et al. 2000;Trumbore 2000). The 14 CO 2 -corrected flux is used to estimate the NSZD rate using the same relationship as presented for measurements of CO 2 flux using the soil gas concentration gradient method (equation (6)).
Additionally, 14 C and 13 C isotopes in soil gas CO 2 , soil and LNAPL were analysed to differentiate between CO 2 generated through biological degradation of naturally occurring soil organic matter and degradation of LNAPL constituents.
Samples were submitted to the Center for Isotope Research at Groningen University (The Netherlands) for isotope analysis. The results of the 13 C analyses are reported using delta notation relative to the international PDB (PeeDee Belemnite) standard, expressed in per mil (‰) units. The results of the 14 C analyses are reported as a fraction of modern carbon (F 14 C), where modern carbon is defined relative to a 1950 baseline.

Field data acquisition
NSZD field data acquisition comprised the following: • quarterly measurements of subsurface O 2 and CO 2 concentration profiles from soil vapour probes (SVPs) installed at five locations (see Fig. 1) and existing monitoring wells screened within the unsaturated zone to estimate NSZD rates using the soil gas concentration gradient method (Johnson et al. 2006;Davis et al. 2009); • quarterly subsurface temperature measurements (15 wells) and continuous subsurface temperature measurements in three wells to support analysis using the biogenic heat method ( Fig. 1), designed for application at sites with impervious ground cover conditions.
SVPs and CO 2 traps were located within 1-2 m of monitoring wells in the background area (SV5 and C5 located adjacent to well S-21), gasoline area (SV3 and C3 located adjacent to well S-5) and diesel area (SV1 and C1 located adjacent to well S-7).

CO 2 provenance
The concentrations of 14 C and 13 C in soil and LNAPL samples, and the concentration of 14 C in soil gas were evaluated to differentiate between CO 2 generated through biological degradation of naturally occurring soil organic matter and degradation of LNAPL constituents (Table 2). Among the soil vapour samples collected at the site, the highest concentration of 14 C is reported for the background SV5 (0.696 at 1.2 m b.g.s.). Although this result indicates that contributions from petroleum degradation in the background location are lower than in source areas, the background concentration is much lower than modern atmospheric levels, which are slightly above unity (ASTM International 2016; Turnbull et al. 2017). The background value is very close to the F 14 C value for soil gas collected from SV3 in the gasoline source area (0.621 at 6.8 m b.g.s.). The similarity between F 14 C measured in soil gas from the background and gasoline areas suggests that the low F 14 C value observed at the background location may be due to lateral migration of soil gas beneath the asphalt pavement. Similar results could also be obtained for degradation of organic matter c. 3000 years old at the background location (SV5) and very limited biodegradation of hydrocarbon at SV3. Field soil gas measurements and CO 2 trap results discussed below provide additional lines of evidence pointing to lateral soil gas transport.
The F 14 C result for the CO 2 sample collected from SV4 (0.174 at 1.5 m b.g.s.) is indicative of shallow contamination (fossil carbon) within the UST secondary containment area (Fig. 2). This surface contamination is consistent with the O 2 /CO 2 profile measured in the area, as discussed in the 'Soil gas concentration gradient method' results below.
F 14 C values of 0.350 (SV1 at 5.2 m b.g.s.) and 0.252 (SV2 at 10 m b.g.s.) for the soil gas CO 2 samples collected from the diesel source area are consistent with the presence of a fossil source related impact at depth in the unsaturated zone. The low F 14 C value of 0.057 for carbon in the LNAPL sample (monitoring well S-7) is consistent with a fossil carbon source, possibly slightly altered by additives with higher 14 C content, such as biodiesel components derived from modern plant biomass (Thomas et al. 2017), and the 13 C value of the LNAPL is typical for a fuel. Carbonates are widely present in the soil matrix of the site, as reflected by the δ 13 C values close to zero (Aelion et al. 2010). Obtaining conclusive information about recent organic carbon fractions based on 14 C results of the soil samples is hindered by the high carbon content documented, which is mainly related to inorganic carbon.

Soil gas concentration gradient method
Under conditions of vertical, 1D gas transport, NSZD results in decreasing concentrations of O 2 and increasing concentrations of CO 2 with depth. However, soil gas profile data collected from several SVPs installed in the gasoline and diesel areas exhibited higher concentrations of O 2 and lower concentrations of CO 2 at intermediate depth intervals. This observation suggests that the standard assumption of vertical (1D) soil gas transport (Fig. 4) at the study site is not a good simplification. Multiple lines of evidence support the conclusion that lateral soil gas transport was occurring in the unsaturated zone at the study site, as follows.
• O 2 content in soil gas increased and CO 2 content decreased with depth at the SVPs and monitoring wells from the shallow intervals to depths of c. 5 m b.g.s. • CO 2 was consistently present at concentrations between c. 1.3 and 2.2 vol% and O 2 between 18.3 and 20 vol% over the full vertical profile at the background locations (S-21 and SV5). The lack of O 2 depletion and/or CO 2 increase with depth is inconsistent with O 2 being supplied by ingress at the ground surface. • Whereas CO 2 concentrations at the background location were very low, 14 C concentration in soil gas collected from the background location (SV5 at 1.2 m b.g.s.) was similar to 14 C concentration in soil gas collected from the deep probe at SV3 (Table 2), consistent with a portion of the CO 2 from each location being derived from petroleum biodegradation.
Soil gas data collected from monitoring wells screened across the water table indicated that there was a transition to a vertical soil gas transport regime near the base of the unsaturated zone, where consistent decreasing O 2 and increasing CO 2 concentrations were observed with depth (Concawe 2020). Thus, biodegradation rates were estimated using vertical concentration gradients for O 2 and CO 2 with data collected from the lower two sampling points (c. 5 and 8 m b.g.s.) of SVPs and soil gas composition data points collected from the lower 3 m at the base of the unsaturated zone in monitoring wells. Soil gas results for SVP SV4 and monitoring well S-17, located c. 4 m apart in the gasoline area (Figs 1 and 2), are compared in Figure 5.
Vertical sampling points at each SVP location did not provide adequate resolution within the vertical transport regime at the base of the unsaturated zone, and thus NSZD rate estimates obtained from SVPs are probably underestimates of the true values. Moreover, the tank vault structure (Fig. 1) had a potential influence on SV4 and S-17 soil gas profiles.
Average NSZD rates were estimated from O 2 and CO 2 profiles of the lower two probes of the SVPs and samples collected at the depth of the unsaturated zone in monitoring wells where O 2 decreased and CO 2 increased with depth. The estimated biodegradation rates are presented in Tables 3 and 4, respectively. The average NSZD rates estimated from O 2 gradients were higher than rates estimated from CO 2 gradients, particularly in the diesel area. The observation that NSZD rates calculated from CO 2 gradients are uniformly lower than rates calculated from O 2 gradients indicates that more O 2 is being utilized than would be predicted from the measured increase in CO 2 using the stoichiometric relationships summarized in equations (1) and (3). The difference between O 2 depletion and CO 2 production cannot be explained by lateral gas transport. Rather, these data suggest that additional chemical reactions may be occurring that consume CO 2 produced during hydrocarbon biodegradation (e.g. carbonate buffering), or that consume O 2 without producing CO 2 (e.g. oxidation of reduced soil minerals). Given the documented presence of carbonate-rich soil and increased alkalinity of groundwater within and downgradient of the source areas at the site (Concawe 2020), it is likely that CO 2 produced from biodegradation of petroleum hydrocarbons may be reacting and going into solution. Similar discrepancies between O 2 utilization and CO 2 production rates have been documented at a number of sites in the USA where in situ respiration tests were completed in the  . Soil gas profiles for soil vapour monitoring point SV4 and monitoring well S-17, located adjacent to one another in the gasoline area. Poor resolution of O 2 and CO 2 concentration gradients was observed using the fixed-depth SVP data, whereas profiling in wells provided higher resolution of gradients near the base of the unsaturated zone. 1980s and 1990s, attributed to the formation of carbonates from the evolution of CO 2 produced by biodegradation (Leeson and Hinchee 1996). The NSZD rates estimated from soil gas data collected from monitoring wells were generally higher than rates estimated from SVPs. Using monitoring wells screened over the interval of interest in the unsaturated zone to measure soil gas concentrations provides a better understanding of vertical gas diffusion gradients and can be used to select an appropriate depth interval where the assumption of vertical gas transport is valid (i.e. where lateral transport is absent) and calculate a more representative NSZD rate. Whereas SVPs are suitable for sites where uniform, vertical soil gas transport dominates, data from this site demonstrate that SVPs installed at pre-determined depth intervals may be insufficient for resolving concentration gradients near the base of the unsaturated zone at sites where lateral gas transport is important.

Biogenic heat method
Plots of temperature data versus time for wells S-5, S-7 and S-21 in gasoline, diesel, and background areas, respectively, instrumented with data loggers are presented in Figure 6. Temperatures measured near the base of the unsaturated zone or near the water table (>8 m b.g.s.) in source zone wells (S-5 and S-7) were generally higher than temperatures measured at the same depth from the background location. This observation is consistent with the soil gas data interpretation, indicating that biodegradation is occurring in this depth interval. On average, temperatures measured in the upper few metres of the subsurface (above the zone of more intense biodegradation) were similar in the background and source zone locations S-5 and S-7. The closer agreement between temperatures at S-21, S-5 and S-7 in the upper few metres of the subsurface is consistent with the lack of shallow impacts observed at S-5 and S-7, and suggests that data collected from S-21 are suitable for understanding background temperature distribution resulting from seasonal variability.
Manual temperature measurements are compared with the data logger readings for S-5, S-7 and S-21 in Figure 6. Whereas the manual measurements recorded near and within the saturated zone generally agree with the data logger results, the manual readings from the unsaturated zone differ from the data logger measurements. Further, the manual readings are consistently warmer than the data logger readings during periods of higher air temperature, and cooler during periods of lower air temperature. This observation suggests that the thermistor may not have reached thermal equilibrium with the subsurface in the unsaturated zone. The annual mean temperature profiles from quarterly manual measurements are compared with manual measurements at the background location in Figure 7. Comparison of the annual average background temperature profile and the annual average temperature profiles observed in the monitoring wells located in the source areas shows an anomalous behaviour, as higher temperatures were measured in the shallow soil column (<6 m b.g.s.) in the background well compared with source area wells.
There is no known hydrocarbon source in the background location that could affect the temperature profile between the surface and the depth of 6 m b.g.s. The differences of solar exposure between the background and the source areas are not considered significant, as only those monitoring wells closer to buildings are affected by a diminished solar exposure.
The observation of elevated temperatures in the source area provides a strong line of evidence that NSZD is occurring. In all cases, the maximum background-corrected (elevated) temperature occurred at 10 m b.g.s., which was the deepest temperature measurement at the background well (S-21). Therefore, the dataset did not permit direct estimation of downward thermal transport. Thus, estimates of heat flux include only the heat from biodegradation that is transmitted upward by conduction to the ground surface. This approach ignores downward and lateral heat flux away from zones where heat is being generated, which probably results in an underestimate of the biodegradation rates.
Given the potential bias in manual temperature measurements recorded in the unsaturated zone, it was assumed that the corrected background temperature at the ground surface (DT (z ¼ 0)) was zero, and upward thermal gradients were estimated by dividing the maximum temperature difference by the depth where the maximum annual average temperatures were observed (10 m), as described by Warren and Bekins (2015). NSZD rates estimated from discrete, quarterly measurements ranged from 2000 to 2700 l ha −1 a −1 , with an average value of 2500 l ha −1 a −1 in the gasoline area and from 1900 to 3200 l ha −1 a −1 with an average value of 2700 l ha −1 a −1 for the diesel area.
For the continuously recorded temperature data, the annual average temperatures recorded at each depth interval were plotted against the background temperature profile to calculate average upward thermal gradients. Temperature profiles along with fluid levels, well construction details and results of soil screening during well installation using a PID are presented in Figure 8.
Average temperatures measured using data loggers are similar to those obtained from the quarterly snapshot measurements, with a maximum temperature increase of c. 2.1-2.6C°relative to background observed at 10 m b.g.s. The upward thermal gradients shown in Figure 8 over the depth range of 5-10 m b.g.s. are probably representative of biodegradation occurring at the base of the unsaturated zone. However, temperature data from S-5 exhibit increased temperatures relative to background in the shallow unsaturated zone that may reflect the presence of shallow impacts, consistent with soil gas concentration gradient and CO 2 trap results from the gasoline area.

CO 2 trap method
NSZD rates were calculated from CO 2 flux data collected from CO 2 traps deployed in autumn and spring at the study site, as shown in Table 5.
The highest NSZD rates were measured in the gasoline source area, whereas the lowest rates were observed within the diesel source area. The standard assumption for interpreting CO 2 flux results is that uniform, vertical gas exchange between the unsaturated zone and atmosphere occurs, and the CO 2 flux measured at the ground surface represents the depth-integrated biodegradation rate. As noted previously, the entire site is paved with concrete in the central portion of the site, where the majority of gasoline and diesel impacts are present, and asphalt in the surrounding areas (Fig. 1). Moreover, the tank vault structure (Fig. 2) had a potential influence on C4 (near SV4 and S-17). Estimated gas diffusivity coefficients for competent asphalt and concrete pavement are of similar magnitude (Comité Euro-International du Béton 1990; Peng et al. 2012;Yoon 2018) and are generally more than an order of magnitude lower than the average diffusivity estimated for site soils. However, cracks and joints appear to be more abundant in the concrete pavement areas, and thus there is probably greater exchange of gases between the subsurface and atmosphere in these areas. This may explain why there is evidence of significant lateral diffusion of petroleumderived CO 2 at the background trap location (where asphalt pavement is present), but not at the CO 2 traps located in the diesel area.
As discussed with respect to interpretation of O 2 and CO 2 concentration gradients measured in the unsaturated zone, the presence of inorganic carbon in the form of carbonates in soil may complicate interpretation of the CO 2 data, as CO 2 produced from biodegradation of petroleum hydrocarbons at depth may react and go into solution within a bicarbonate-buffered system. In these circumstances, CO 2 traps provide results that need to be evaluated based on site-specific conditions, where • the higher biodegradation rates in the gasoline area may represent degradation of shallow contamination immediately below the paved surface that would be less likely to be influenced by CO 2 consumption by carbonate minerals; • the low biodegradation rate in the diesel area may be due to lateral offset of CO 2 diffusing upward from the source, or reactions between CO 2 produced at depth with carbonate minerals in soil and porewater; • the higher biodegradation rates observed at the background location are probably the result of lateral vapour transport from areas of the site where hydrocarbon impacts are present.
It appears that the CO 2 trap deployment method using the vapour pins to equilibrate pressure and avoid the stack or chimney effect that often results from the contrast in gas diffusivity between soil and pavement was effective. Other than the high CO 2 flux measured at location C4 in the gasoline area in October 2017, the measured NSZD rate results do not appear to be biased high relative to NSZD rates reported in literature (e.g. Garg et al. 2017). Although the approach appears to have mitigated the stack effect, the CO 2 trap method includes a simplifying assumption of 1D gas transport (i.e. it is assumed that the CO 2 measured at the surface is derived from petroleum present in the subsurface directly below the trap). However, in this case there is probably lateral transport of CO 2 to fissures in the pavement. Additionally, CO 2 produced through biodegradation of petroleum hydrocarbons near the base of the unsaturated zone may be consumed by reactions with carbonate minerals before reaching the traps. These conditions undermine the reliability of CO 2 traps for quantifying NSZD at the site. However, the observation that petroleum-derived CO 2 is being produced is a direct line of evidence that NSZD is occurring.
Comparison of NSZD rates from soil gas concentration gradient, biogenic heat, CO 2 trap methods Average NSZD rates by the three methods applied at the site are presented in Table 6, which shows a higher than expected variability in NSZD rates determined per method. Variability in NSZD measurements are not unique to the study site. Guidance such as API (2017) and CRC CARE (2018) have identified a variety of Fig. 6. Summary of temperature data recorded using data loggers and manual thermistor readings in background (S-21), gasoline (S-5) and diesel (S-7) area wells. Gaps in temperature data were modelled using a best-fit sinusoidal function, shown as dotted lines.
factors that can affect NSZD measurements and suggested that NSZD rate measurements should be considered to have one order of magnitude precision. The variability between methods in this study was above an order of magnitude. CRC CARE (2018) went as far as providing summary tables describing the potential effects of certain site conditions on NSZD measurements (Table 2; CRC CARE 2018) and comparing attributes of different measurement methods to support method selection (Table 6; CRC CARE 2018). Here, we compare the methods deployed for this study to give a sense of the accuracy and precision of the results. For this evaluation, accuracy is evaluated based on whether the primary assumptions of the analytical models  for each method were met at the site. Where assumptions for a given method appear consistent with site data, precision is assessed on the basis of temporal and spatial repeatability of the results.
Soil gas composition data showed that the fundamental underlying assumption of vertical soil gas transport that underpins the soil gas concentration gradient method was not valid in the upper 5 m of the unsaturated zone. The CO 2 trap results were also probably affected by lateral soil gas transport, as evidenced by the depleted 14 CO 2 signal observed in the trap placed in the background area of the site. The elevated NSZD rates measured from CO 2 traps in the gasoline area (17 000-130 000 l ha −1 a −1 ) were probably the result of CO 2 contributions from biodegradation of shallow, residual LNAPL. Results for CO 2 traps deployed in the diesel area, where shallow LNAPL was not present, were much lower (0-1000 l ha −1 a −1 ). The limited CO 2 flux measured in the diesel area was also suspected to be, at least in part, the result of CO 2 buffering from bicarbonate in porewater. The lower rates estimated from the soil gas concentration gradient method for CO 2 as compared with O 2 support this hypothesis.
Although there is evidence that lateral gas transport in the upper 5 m of the unsaturated zone probably affected the accuracy of soil gas concentration gradient measurements made using SVPs, soil gas composition profile data from monitoring wells indicated that there was a transition to a vertical soil gas transport regime at the base of the unsaturated zone (>5 m). The observed monotonic decrease in O 2 and increase in CO 2 concentrations at these depths are consistent with the assumption of 1D gas diffusion. Microbial generation of heat is affected by gas transport mechanisms; aerobic oxidation of hydrocarbons and/or methane occurs where hydrocarbon impacts (and/or methane) and O 2 are both present in the unsaturated zone. However, uniform vertical soil gas flux is not a core underlying assumption in the calculation method. The primary assumption for the biogenic heat method is that all heat transfer processes other than those associated with NSZD are approximately equivalent for background and source zone locations, allowing for clear identification of the NSZD heat signal. For this study background-corrected temperatures measured close to the water table, where more intense biodegradation was expected, ranged from c. 2 to 3°C, whereas background-corrected temperatures within the upper 2 m of the subsurface were close to zero in wells instrumented with data loggers (S-5 and S-7). The relative agreement between temperatures measured at the background and source zone locations within the upper few metres of the subsurface suggest that the background data are adequate for identifying heat transfer associated with NSZD at the study site. Greater variability in shallow, background-corrected temperatures was observed for the manually collected temperature profiling data, many of which indicated negative background-corrected temperatures in the upper 4-7 m of the subsurface for source zone locations, but had similar background-corrected temperatures to the wells with data loggers from 8 to 10 m b.g.s.. More sophisticated approaches for modelling background temperature distribution were considered to refine biogenic heat estimates. However, all published backgroundcorrection approaches known to the authors (e.g. de Vries 1963; Sweeney and Ririe 2014;Askarani et al. 2018) assume a zero net annual heat transfer to the subsurface, and/or fixed annual mean background temperatures at depths where seasonal temperature variations are negligible. The data collected in this study (summarized in Figs 6-8) indicated that the mean temperature with depth was cooler than the mean temperature over the study period as well as the longer-term climatic mean. This indicates that there is a net downward heat flux in this area. Use of the background should have accounted for this but highlights that this method is relatively new and as such the theoretical basis does not consistently match site conditions. Multiple sites have experienced issues with background surface temperatures not corresponding to ideal behaviour (Askarani et al. 2018;CRC Care 2020a).
Of the measurement approaches applied, methods that relied on subsurface measurements using existing monitoring wells (soil gas concentration gradient method and biogenic heat method) provided results that were consistent with method assumptions. The spatial and temporal variability in results for these methods provides a general indication of method precision. Results presented in Table 6 show the range in annual mean values for both methods collected from 13 locations. These data indicate a range of values observed for the soil gas concentration gradient method that vary by a factor of about five (440-2000 l ha −1 a −1 ) and a slightly lower range of rates for the biogenic heat method that vary by a factor of less than two (1900-3200 l ha −1 a −1 ). The slightly lower spatial variability observed for the biogenic heat method may be explained by soil heterogeneity. Soil gas diffusion rates are sensitive to subtle changes in soil texture and fluid saturations, whereas conductive heat transport is relatively consistent for all soil types as the thermal conductivity of the soil is a function of the minerals that make up soil grains and the water content in the soil matrix. Temporal variability, as indicated by differences in rate estimates over time at the same location, varies by a factor of three or less for both methods. These  (12000) ranges for both methods are relatively narrow compared with the order of magnitude precision suggested by CRC CARE (2018), which suggests that both methods probably provide reasonable estimates of NSZD rates at the site. Although accuracy concerns have been identified for the carbon traps and data from the shallow SVPs, it is worth noting that results from all methods provided valuable lines of evidence that LNAPL degradation is occurring at the site. In these authors' experience, the magnitude of NSZD rates has not been used to make site decisions when LNAPL and associated dissolved and vapour phase bodies are stable or decreasing, and risk to receptors is appropriately addressed. Quantification of rates has contributed to the industry understanding the significance of these natural remedial mechanisms. However, where potential land use changes or other factors do not require an accelerated pace of remediation, the uncertainty associated with the observed variability in this study would not necessarily change site decisions. Even in scenarios where transition from active remediation to NSZD is being contemplated, the magnitude of rates may not be a key metric unless there are specific time constraints. Rather, the stability of LNAPL and associated dissolved and vapour phase plumes, and assurance that risks can be appropriately managed by NSZD, will more probably be criteria.

Conclusion
Where NSZD needs to be demonstrated as occurring, any of these methods would have provided sufficient justification. Where more accurate and precise estimates are desired, the variability in NSZD rate estimates using the different methods applied in this study demonstrates the need for careful examination of the conceptual site model and consideration of potential interferences prior to selecting a single method for NSZD measurements at paved sites. Although the study results demonstrate that the presence of pavement was not inhibiting O 2 ingress or limiting biodegradation of LNAPL constituents, several site-specific conditions that contributed to the observed variability between measurement methods were identified, including the following: • evidence of lateral soil gas transport in the upper 5-6 m of the unsaturated zone; • utilization of O 2 and production of CO 2 within multiple depth horizons owing to the presence of localized shallow, residual LNAPL in the top few metres of the subsurface in addition to LNAPL present at the water table; • indicated reactions between CO 2 produced from biodegradation of petroleum hydrocarbons with carbonate minerals in soil and porewater.
The results from this case study highlight the utility of measurements in existing monitoring wells screened across the base of the unsaturated zone. At sites with existing, appropriately screened monitoring wells located within the source zone, an initial screening approach that includes temperature profiling and measurement of soil gas composition (O 2 , CO 2 and CH 4 ) at the base of the unsaturated zone (Ririe and Sweeney 2018) provides an opportunity to obtain useful NSZD data without the need for additional, intrusive subsurface work. Depending on site-specific needs, these data may provide adequate information to estimate NSZD rates. The results may also provide valuable information for selecting method(s) that are best suited for future or long-term NSZD measurements, if necessary.