Natural Source Zone Depletion (NSZD) Quantification Techniques: Innovations and Future Directions

Abstract

Natural source zone depletion (NSZD) is an emerging technique for sustainable and cost-effective bioremediation of light non-aqueous phase liquid (LNAPL) in oil spill sites. Depending on regulatory objectives, NSZD has the potential to be used as either the primary or sole LNAPL management technique. To achieve this goal, NSZD rate (i.e., rate of bulk LNAPL mass depletion) should be quantified accurately and precisely. NSZD has certain characteristic features that have been used as surrogates to quantify the NSZD rates. This review highlights the most recent trends in technology development for NSZD data collection and rate estimation, with a focus on the operational and technical advantages and limitations of the associated techniques. So far, four principal techniques are developed, including concentration gradient (CG), dynamic closed chamber (DCC), CO₂ trap and thermal monitoring. Discussions revolving around two techniques, “CO₂ trap” and “thermal monitoring”, are expanded due to the particular attention to them in the current industry. The gaps of knowledge relevant to the NSZD monitoring techniques are identified and the issues which merit further research are outlined. It is hoped that this review can provide researchers and practitioners with sufficient information to opt the best practice for the research and application of NSZD for the management of LNAPL impacted sites.

1. Introduction

Natural source zone depletion (NSZD) includes natural loss of light non-aqueous phase liquid (LNAPL) hydrocarbon contamination in oil spill sites. With the current focus of the global market on fossil fuel as the main source of energy [1], oil spills are inevitable. Oil spills have been happening from the 1900s [2,3,4], when a gradual shift from coal to oil consumption occurred [5], to this date [6,7]. For example, South Korea suffers from over 300 cases of marine oil spill incidents annually [8]. On 21 July 2016, a damaged Husky Energy pipeline (16 TAN line) resulted in leakage of 225 m³ heavy crude oil, of which about 40% ended up in the North Saskatchewan River [7].

Management of oil spill sites cannot be ceased even after a major shift to renewable sources has happened. Evidences show that petroleum hydrocarbon (PHC) contamination lasts for decades to over a century [3,4,9,10,11,12]. A zero-order NSZD rate model (The Glide Path Model developed by [13]) suggested that on average, LNAPL longevity can range 35 to 105 years, with plus or minus factors of 2 and ½ respectively [13]. Stable isotope analysis estimated a minimum life expectancy of 110 years for an oil plume near Bemidji, Minnesota [14].

Compared to the conventional “active” remedial techniques, NSZD can be an efficient and sustainable soil and groundwater bioremediation alternative. Some examples of active systems include: excavation (dig and haul, D&H), skimming, multi-phase extraction (MPE), water flooding, surfactant-enhanced subsurface remediation (SESR), air sparging/soil vapor extraction (AS/SVE), physical or hydraulic containment, in situ chemical oxidation, in situ soil mixing (stabilization) and steam injection [15,16]. These remediation strategies are typically costly, energy intensive and environmentally disruptive [17] and often leave behind significant LNAPL residues [18]. More importantly, they are generally associated with large greenhouse gas (GHG) emissions. Clean-up of 1 kilogram (kg) of soil and groundwater contaminants using these methods can emit up to 5 and 130 tonnes of CO₂ to the atmosphere, corresponding to geometric means of 0.015 and 1.3 tonnes, respectively [19]. To put it into context, in 2019, about one (0.975) billon kg of toxic chemicals was released into land due to industrial activities in the United States (U.S.) only [20]. The issue is compounded given the consistency in release (e.g., substance release of 1.045 billion kg happened in the U.S. about ten years earlier, in 2010 [20]) as well as the persistence of PHC contamination.

The NSZD is a passive LNAPL bioremediation technique; hence, it does not contribute more to GHG emissions than what LNAPL contaminated sites are already emitting to the atmosphere. The NSZD is mainly a method of monitoring the current CO₂/heat emissions and correlating these emissions stoichiometrically with LNAPL decomposition in the subsurface, rather than actively interfering with the LNAPL removal. As such, developing the knowledge and utilizing NSZD technique can closely tie into Sustainable Development Goal (SDG) #13-Climate Change. More specifically, study and application of NSZD bioremediation technique can meet one of the objectives of SDG #13, “improvement of education and institutional capacity on climate change mitigation, adaptation, impact reduction” [21].

Since NSZD seems to be a less expensive remediation technique compared to active methods, it can help to reinforce the capacity of institutions, including governments, businesses and non-governmental groups involved in remediation, decommissioning and restoration research and industry, in ways that they can improve and use the technique, as well as plan and manage LNAPL sites cost effectively.

The NSZD consists of a series of natural processes (refer to Section 2 for further information), with minimal environmental footprints and diverse desirable impacts on the LNAPL body, including decreasing LNAPL mass, saturation and mobility, shortening LNAPL longevity and altering LNAPL composition [22]. Depending on spill magnitude, site conditions and regulatory requirements, the NSZD can either help to reduce the deployment time of active remedial techniques or can be completely substituted for these methods [23]. Although NSZD is an ongoing process in almost all LNAPL release sites, caution must be taken when considering it as the sole or primary remediation strategy. Careful assessment is required to determine if the NSZD occurring at the site is effective enough to meet the regulatory objectives or if it requires the support of auxiliary techniques [24].

This natural phenomenon has certain characteristic features (see Section 2) that can serve as NSZD fingerprints and be surrogates for the quantification of NSZD rate (i.e., rate of bulk LNAPL mass depletion; see Section 2.2 for more information) [22]. The question is how to collect NSZD data with the highest accuracy (close to true values) and precision (reproducible data) in the most cost-effective way, to ensure that continuous monitoring of NSZD is feasible, efficient and economical. Yearlong data collection is critical given the temporal (seasonal and diurnal) and spatial (vertical and horizontal) variability of NSZD rates [25,26] for numerous, inherent reasons (which are discussed in this paper). When NSZD is used as the primary or sole remedial technique, analyzing the longevity and the composition of LNAPL body that both can be estimated by NSZD rates is an integral part of the remediation strategy [22].

This review highlights the most recent trends in technology development for NSZD data collection and rate estimation. The objectives of this literature review and analysis were to provide a state of the knowledge revolving around the available NSZD rate quantification techniques, including an overview of the underlying NSZD concept, definition of NSZD rate, influencing factors and operational and technical advantages and limitations of the techniques. Discussions revolving around two techniques, “CO₂ trap” and “thermal monitoring”, are expanded (see Section 4 and Section 5, respectively), due to particular attention to them in the current industry. The gaps of knowledge relevant to the NSZD monitoring techniques are identified and the issues which merit further research are outlined. A few items were out of the scope of this article including: (i) critical review of the underlying factors and providing thorough explanation of the underlying phenomena; (ii) analyzing techno-economic aspects of NSZD rate quantification techniques; (iii) bibliographic analysis of the topic. These are included in the concluding remarks (Section 6) as a few additional directions for future studies.

A review article with this scope sounds to be an important and a much-needed addition to the current literature. The NSZD is an emerging technology at embryonic stage and the literature currently lacks a document that can provide the audience with a snapshot of the available techniques developed for resolving NSZD rates and the associated gaps of knowledge. The audience of this paper are researchers from Research and Development (R&D) businesses and academia as well as practitioners from remediation and reclamation industry who are interested to learn more about and use the NSZD for regulatory advisory and services. The NSZD has the potential to be used as an effective and sustainable oil spill management technique. However, its rate should be quantified reliably, accurately and precisely, to avoid the risk of regulatory non-compliance resulted from inaccurate rate estimations. Identifying and discussing the gaps of knowledge pertaining NSZD rate quantification is prudent to give direction to academicians for further research and to increase the acceptance of this emerging technique in the industry. It is hoped that this review article can provide researchers and practitioners with enough information to opt the best practice for the research and application of this sustainable LNAPL management technique.

2. Overview of NSZD

The current literature carries a wealth of information about NSZD concept, microbiology and mechanisms [22,25,27,28,29,30,31,32,33,34,35,36,37,38], reviewed by a number of researchers [10,25,36,39,40]. As such, this section only intends to provide an overview of these NSZD aspects, to give context to the next discussions pertaining to NSZD rate quantification techniques. For further information, the readers may refer to the sources introduced above and cited in the following sub-sections.

2.1. NSZD Conceptual Model

Figure 1 illustrates NSZD conceptual model which has been developed based on observations and efforts of a number of researchers [25,27,33,35,40,41,42,43,44,45].

According to evidences gathered about microbial community structure and geochemical properties of a shallow hydrocarbon-contaminated aquifer at a former petroleum refinery, Irianni-Renno et al. [41] identified four distinct zones at an LNAPL impacted subsurface:

• An aerobic, low LNAPL mass zone at the upper part of unsaturated (vadose) zone, extending to ground surface. This is where only soil respiration and modern carbon (modern-C) generation occurs (Region I in Figure 1);
• A moderate to high LNAPL mass zone in the middle of the vadose zone where low-oxygen environment transitions to an anaerobic condition (Region II in Figure 1);
• An anaerobic, high LNAPL mass zone, encompassing the lower part of the vadose zone including capillary fringe and smear zone and the upper portion of saturated zone (Region III in Figure 1);
• An anaerobic, low LNAPL mass zone below LNAPL plume (Region IV in Figure 1).

The term “source zone” was coined by Amos et al. [27] referring to the subsurface area including saturated and unsaturated zones where free-phase oil is present (Regions II & III in Figure 1).

Natural loss of LNAPL occurs in these regions through diverse physicochemical and biogeochemical processes, including volatilization, dissolution and most importantly biodegradation. The extent to which each of these processes occurs greatly depends on geologic and microbial conditions of each region [11,41,46]. LNAPL composition is an important influencing factor as well [39]. Equilibrium simulations conducted on water–oil and air–water systems indicated that volatilization and biodegradation were dominant removal pathways for C₆–C₉ n-alkanes and cyclohexanes and dissolution and biodegradation were principal removal pathways for the other PHCs [9].

The NSZD begins with methanogenesis and anaerobic reduction of nitrate (NO₃⁻), iron (Fe³⁺) or sulfate (SO₄²⁻), if present, in saturated and anaerobic parts of vadose zones (Regions III & IV). A syntrophic consortium of bacteria and archaea (fermenters and methanogens) converts hydrocarbons to acetate (CH₃OOH) and then methane (CH₄) and carbon dioxide (CO₂) [11,41,46]. Fermenters convert LNAPL to acetate (reaction 3) which can be taken up by acetate-using methanogens (reaction 4, known as acetotrophic or acetoclastic methanogenesis). Syntrophic organisms oxidize acetate as well, converting it to hydrogen (reaction 5) which serves hydrogenotrophic methanogenesis (reaction 6) as the main substrate [36].

Biogenic gases generated through biodegradation (i.e., CH₄ and CO₂) are released to the gaseous phase of the subsurface through two fundamental phenomena: direct outgassing and degassing [25,44]. Direct outgassing pertains to less soluble yet biodegradable PHCs (e.g., tri- and tetra-methylbenzenes, naphthalenes and cyclohexanes [27]). Due to low solubility, these hydrocarbons are most likely physically isolated and distant from the groundwater body and are removed through a distinct pathway in soil pore spaces adjacent to the entrapped oil, which is referred to as “oil-contact biodegradation” and is defined as “biodegradation near the oil without significant migration”. The gaseous by-products of oil-contact biodegradation are emitted directly to the gaseous phase, without an intermediate dissolution step into the bulk aqueous body (i.e., the groundwater) [44]. The mechanisms of PHC uptake by microorganisms in the pore spaces are described by Hua and Wang [40] (in particular, refer to Figure 3 of Hua and Wang [40]). Briefly, to take up hydrophobic hydrocarbons, microorganisms “pseudo-solubilize” the hydrocarbons to pass them through their membranes. Pseudo-solubilization of hydrocarbons is performed by producing biological surfactant molecules (in short, biosurfactants), to emulsify hydrocarbons and produce micro-droplets suitable for microbial assimilation [40].

Degassing, in contrast, includes gas–liquid partitioning and bubble formation in the groundwater. Degassing is followed by ebullition which happens when gas saturation increases, bubbles enlarge and buoyancy outcompetes capillary forces [27]. While gas release in the pore spaces in the smear zone happens mainly through direct outgassing, gas release in the saturated zone can happen through both mechanisms [25,44]. Direct outgassing in groundwater is due to the hydrophobic nature of LNAPL which positions them atop the water table, exposing them to the gaseous phase of the aquifer [44].

Garg et al. [25] argues that direct outgassing is different from degassing. While this is true at macroscopic level, these two phenomena seem to share the same principles at microscopic levels. In smear zone, microorganisms reside either in the water droplets entrapped inside the oil patches or the water entrapped in the pore spaces which surrounds the oil droplets/patches. With the help of advanced analytical techniques (e.g., nuclear magnetic resonance and Fourier transform ion cyclotron resonance mass spectrometry), it is found that the microorganisms remain metabolically active in miniscule water droplets (1 to 3 μL) encased within oil [43]. In either case, biogenic gases must partition and bubble out of the liquid phase whether it is water or oil.

Dissolution of CH₄ and CO₂ causes exceeding local liquid-phase gas solubilities, which results in bubble formation, ebullition of gas bubbles and upward gas fluxes [25,46]. Thereafter, advection and diffusion transport the released gases toward the ground surface [47]. According to Irianni-Renno et al. [41], upward flux of CH₄ is most likely the cause of anaerobic zone formation in the vadose zone because it limits downward flux of O₂. As depicted in Figure 1, CH₄ and O₂ fronts converge at the plane Z_A/A in the subsurface. In the presence of oxygen diffused into the subsurface (i.e., O₂ influx), aerobic methanotrophic bacteria oxidize CH₄ as the source of carbon and energy and O₂ depletion occurs (reaction 1). Due to low solubility, CO₂ produced through methanogenesis and methane oxidation diffuses upward to escape to the atmosphere (i.e., CO₂ efflux) [41].

Aerobic processes (reaction 2) can partially contribute to PHC degradation and O₂ depletion. Oxygen can find its way into the subsurface through surface recharge of oxygen-bearing meteoric waters [48] and directly from the atmosphere in gaseous form. Highly soluble benzene, toluene, ethylbenzene and the xylenes (BTEX) and less soluble PHCs such as tri- and tetra-methylbenzenes, naphthalenes and cyclohexanes are susceptible to aerobic and anoxic biodegradation [49]. Nonetheless, in NSZD sites where a significant mass of LNAPL exists, the dominant long-term PHC degradation pathway is supposed to be methanogenesis [27]. Methanogenesis was estimated to contribute 16% to PHC degradation in 38 sites, indicating the general importance of anaerobic degradation in PHC contaminated sites [50]. It agreed with the laboratory study of Salminen et al. [51] which suggested that anaerobic processes played key roles in PHC degradation at a boreal, lightweight fuel and lubrication oil contaminated site. However, it seems that biodegradation can only degrade lighter hydrocarbons. Only 10% of n-alkanes that were depleted through aerobic and anaerobic processes belonged to a C₁₅–C₄₀ group, while ~90% were lighter n-alkanes with shorter carbon chains (C₁₀–C₁₅) [51]. The methanogenic process starts with the degradation of longer chain, biodegradable n-alkanes, working through to the degradation of shorter chain consecutively. Whereas, the aerobic degradation sequence is the opposite [28]. Similar to long chain PHCs, branched alkanes are refractory, tending to remain in source zone for a longer period than n-alkanes [52,53]. Garg et al. [25] suggest that aerobic degradation of PHC is likely an important process in the vadose zone of sites with recent spills. On the contrary, at sites with relatively old releases, an anaerobic core forms in the vadose zone [25]. Aerobic processes are likely predominant in shallow depths (e.g., up to 1.4 m below ground surface, bgs) on NSZD sites [41] where major oxygen sink in an unsaturated zone is linked to methane oxidation [27].

Aerobic and anaerobic degradations of LNAPL and the associated metabolites (e.g., methane oxidation) produce heat (▲ symbol in Figure 1). Biogenic heat can be used as a signature of these exothermic processes occurring at NSZD sites (read more in Section 5). Heat produced through aerobic methane oxidation and aerobic LNAPL degradation varies between 45–50 KJ/g, while anaerobic LNAPL degradation generates 300–800 times less energy (0.06–0.16 KJ/g), depending on the representative PHC chosen for conceptual model construction. In thermal monitoring technique (read more in Section 4.4), this heat is measured and stoichiometrically translated to NSZD rate [11,12,23,25,41,46].

It should be noted that what is described here only constitutes the components of a general model. Prior to customizing this model for any NSZD site, it is imperative to develop a site-specific LNAPL Conceptual Site Model (LCSM) through careful assessment of the site geology (soil type, soil particle size, subsurface fractures, stratigraphy, etc., which can limit or enhance LNAPL migration in the vadose zone), hydrostratigraphy (groundwater characteristics, such as flow pattern, capillary fringe and seasonal water table fluctuations), LNAPL spatial distribution (horizontal extent, vertical extent, smear zone, dissolved-phase plume and vapor plume) and residual LNAPL mass and composition [54,55]. Site-specific LNAPL fingerprinting and LCSM development is a critical step before embarking NSZD technique for LNAPL site management, as LNAPL compositional changes (i.e., transitioning between different dominant classes of chemicals) influence short-term NSZD rates by shifting between dominant processes (i.e., volatilization, dissolution and biodegradation) at different site stages, which consequently dictate long-term NSZD rates, overall LNAPL mass loss and the transition points [39]. Other information that should be considered for LCSM development include: (i) potential exposure pathways and receptors; (ii) potentials of explosivity; (iii) LNAPL recoverability [54,55], which define the longevity of risk factors. As an example, the development of the LCSM model is practiced in Sihota et al. [26].

2.2. Definitin of NSZD Rate

NSZD rate is defined as the rate of bulk LNAPL mass depletion [22], expressed as LNAPL volume degraded per unit area per year (in U.S. gallons/acre/year or liters/hectare/year) [56]. NSZD rate is affected by many factors, which are discussed in the next section. To understand how these factors affect NSZD rates, a distinction should first be made between actual and apparent NSZD rates.

Microbial NSZD degradation rate is defined as “actual NSZD rate” [13,46]. The sequence of anaerobic methanogenesis and aerobic methane oxidation is the major PHC biodegradation pathway at NSZD sites, which seems to contribute over 99% to overall LNAPL volume degraded per unit area per year. The remainder (<1%) is attributed to biodegradation capacity of electron acceptors, such as oxygen, nitrate, iron and sulfate. Methanogenesis is the rate-limiting process and thus it is supposed that the actual NSZD rate is roughly equal to the methanogenesis rate [25].

There is a time lag (5–7 months) between the occurrence of the actual NSZD rate and the manifestation of NSZD as surficial CO₂ efflux. This is attributed to CO₂ travel time in the subsurface [26]. This concept is extrapolated to thermal monitoring quantification technique as well, assuming that a gap also exists between the actual rate and thermal expression of NSZD due to heat transport mechanisms [46]. As such, an additional parameter, “apparent NSZD rate”, is defined to account for this delay [26,46]. Drawing from the information in the literature [25,26,46], apparent NSZD rate can be defined as the rate associated with a measurement performed at a certain time with a certain quantification technique.

2.3. Factors Affecting NSZD Rates

Factors regulating actual and apparent NSZD rates are shown in Figure 2 and are divided into three categories:

• Factors regulating actual NSZD rate;
• Factors regulating apparent NSZD rate;
• Factors regulating actual and apparent NSZD rates simultaneously.

Figure 2 depicts a brief description of these factors and how they regulate NSZD rates. Appendix A provides the state of the knowledge about these factors available in the current literature. The readers interested in additional information are referred to this Appendix A.

It seems that interactions exist between these factors. These interactions are shown in Figure 3 with dash lines. The anticipated principles of these interactions are underpinned in the discussions provided in Appendix A.

The entangled effects of factors can explain temporal and spatial heterogeneity in the subsurface, which results in variability of the NSZD rate from site to site, spot to spot (in a single site) or even time to time (at one spot). The degree of subsurface heterogeneity and NSZD rate variability depends on how many factors are interacting at a certain time in a single location. This emphasizes the challenges that one might face in driving site-specific LCSM and NSZD conceptual models. If the effects of factors affecting apparent NSZD rate are minimized, actual and apparent NSZD rates converge and the accuracy of the relevant quantification techniques improve.

2.4. Microbiology of NSZD Sites

Understanding microbial ecology of NSZD sites is critical for optimizing LNAPL bioremediation plans and developing NSZD conceptual and LNAPL longevity models [41]. Many factors affect microbial NSZD degradation rates (i.e., actual NSZD rate; see Section 2.2 for the difference between actual and apparent NSZD rates) which in turn regulates LNAPL longevity. Those factors are thoroughly discussed in Section 2.2.

2.4.1. Microbial Ecology in Recent and Aged Oil Spill Sites

During the early stages of petroleum release, indigenous microbial communities are most likely not equipped with the metabolic tools (e.g., enzymatic system, bio-surfactants producing abilities, etc.) or physical arrangements (e.g., biofilm structure, physical distance from LNAPL body) to decompose LNAPL. Over time, these abilities are developed through shifting the microbial community structure [13,62,63]. The rate of acclimation greatly depends on: (i) LNAPL composition; for example, the degradation of paraffins (n-alkanes) and non-n-alkanes (branched and aromatic constituents) with carbon numbers of C₁₀ and less are metabolically preferred to non-n-alkanes with C₁₁ and more [11]; (ii) inhibition by certain LNAPL components (e.g., gasoline [55]) and other sources (e.g., volatile hydrocarbons (VHCs) and acetate [25]); (iii) temperature; seasonal fluctuations in temperature equal to or greater than 4 °C can affect NSZD rates [64]. Due to subsurface temperature changes caused by LNAPL degradation, NSZD can be a self-stimulatory process [13,64]. Hence, the acclimation phase can last days, weeks, or even months based on site-specific characteristics [13], during which hydrocarbon oxidation is concomitant with microbial growth [52].

NSZD mathematical simulations are mostly based on a biofilm concept, with a few models considering motile or suspended microorganisms [39]. Microbial agglomeration in porous media (e.g., subsurface soil and rocks and riverine hyporheic zones) has been observed and studied, albeit mostly without the NAPL phase [65,66,67,68]. Mapping physiologic zones of an LNAPL plume suggested that only 15% of the total population was in suspended form [69]. Through scanning electron microscopy (SEM), extensive biofilm development was observed on the mineral surfaces of an LNAPL contaminated site adjacent to a former refinery in Carson City, Michigan [10,70]. In a laboratory experiment, biofilm development was observed in a bio-stimulated column which was packed with 20–30 mesh silica sand and fed with nutrient broth and diesel fuel. Bio-stimulation included seeding the column with a mixed bacterial culture that was retrieved from a refinery contaminated site and had hydrocarbon-degrading strains from Variovorax and Stenotrophomonas genera [65].

Biofilm formation reduces soil porosity and alters soil tortuosity due to biofilm growth and sloughing (pore clogging). This problem is compounded by the presence of LNAPL which have adverse effects on electrical conductivity and permittivity of soil. LNAPL have high electrical resistivity and act as insulators (0.001 mS/m) when compared to highly conductive, natural formation fluids, such as freshwater (0.1–1 mS/m). LNAPL also show low permittivity (e.g., 2–3) relative to those of clean water (~80) and water-saturated sands (20–30) [10,63]. Spatial distribution of attached microorganisms influences the structure of interfaces, the connectivity of pore spaces and the bulk geometric properties of the medium [67]. Interfacial properties of the sediment grains are altered by microbial colonization and biofilm development, production of extracellular polymeric substances (EPS) and formation and precipitation of mineral complexes (e.g., metal sulfides) due to microbial activity [71,72]. Biofilm distribution can be either continuous or patchy (in the latter form, they mainly accumulate in pore throats), depending on physical, chemical and biological properties and history of the medium [73,74,75]. In an LNAPL impacted subsurface, site history can be divided into two fundamental eras: before and after LNAPL release. The combined effects of the LNAPL presence and biofilm formation modify subsurface mass and momentum transport dynamics, potentially leading to biofilm isolation from LNAPL plume [13,63]. As a consequence, not only geologic features but also microbial characteristics of biofilm-dominated zones in an LNAPL contaminated subsurface can be distinct from those areas conquered by motile, free-floating microorganisms which can actively “chase” food regardless of flow directions or may precipitate to the bottom of the aquifer due to Stokes law [39,63].

In the mid and late stages of LNAPL release, thermodynamic principles and diauxie may govern microbial ecology of NSZD sites. Coined by Jacques Monod, the term “diauxie” refers to preferential microbial metabolism in a dual-substrate medium which promotes a relatively fast and less energy-intensive growth pathway [13,62]. When more than two substrates are available, microbial culture encounters a “multiauxic” situation [76]. In multiauxic media, a cascade of reactions happens from the most to the least thermodynamically favorable substrate [13,62,77,78]. The example provided above about the comparison between paraffins and non-n-alkenes with C₁₁ and more is relevant to this concept as well. Thermodynamics of NSZD sites are discussed in the next section.

2.4.2. NSZD Thermodynamics

Among possible metabolic pathways in an LNAPL contaminated subsurface (

Table 1. Thermodynamics of LNAPL biodegradation at NSZD sites through different possible pathways. Decane (C₁₀H₂₂) is chosen as model hydrocarbon (compiled based on Interstate Technology and Regulatory Council (ITRC) [22]; Karimi Askarani and Sale [81]; Stockwell [78]).

Rank a	Redox Reactions	Equation #	Free Energy b ΔGro (kJ/mol-C10H22) c	Heat Exchange b ΔHro (kJ/mol-C10H22) c
1	Aerobic methane oxidation: $7.75{ \mathrm{CH}}_4+15.6{ \mathrm{O}}_2\to 7.75{ \mathrm{CO}}_2+15.5{ \mathrm{H}}_2\mathrm{O}$	Equation (1)	−6696 d	−7444 d
2	Aerobic respiration: ${\mathrm{C}}_{10}{\mathrm{H}}_{22}+15.5{ \mathrm{O}}_2\to 10{ \mathrm{CO}}_2+11{ \mathrm{H}}_2\mathrm{O}$	Equation (2)	−6505	−6978
3	Denitrification: ${\mathrm{C}}_{10}{\mathrm{H}}_{22}+12.4{ \mathrm{NO}}_3^{-}+12.4{\mathrm{H}}^{+}\to 10{ \mathrm{CO}}_2+17.2{ \mathrm{H}}_2\mathrm{O}+12.4{ \mathrm{N}}_2$	Equation (3)	−6369	−6314
4	Manganese reduction: ${\mathrm{C}}_{10}{\mathrm{H}}_{22}+31{ \mathrm{MnO}}_2+62{\mathrm{H}}^{+}\to 10{ \mathrm{CO}}_2+42{ \mathrm{H}}_2\mathrm{O}+31{ \mathrm{Mn}}^{2+}$	Equation (4)	−6556	−6559
5	Iron reduction: ${\mathrm{C}}_{10}{\mathrm{H}}_{22}+62 \mathrm{Fe}{\left(\mathrm{OH}\right)}_3+124{\mathrm{H}}^{+}\to 10{ \mathrm{CO}}_2+166{ \mathrm{H}}_2\mathrm{O}+62{ \mathrm{Fe}}^{2+}$	Equation (5)	−4441	−5160
6	Sulfate reduction: ${\mathrm{C}}_{10}{\mathrm{H}}_{22}+7.75{ \mathrm{SO}}_4^{2-}+15.5{\mathrm{H}}^{+}\to 10{ \mathrm{CO}}_2+11{ \mathrm{H}}_2\mathrm{O}+7.75{ \mathrm{H}}_2\mathrm{S}$	Equation (6)	−951	−230
7	Methanogenesis: ${\mathrm{C}}_{10}{\mathrm{H}}_{22}+4.5{ \mathrm{H}}_2\mathrm{O}\to 2.25{ \mathrm{CO}}_2+7.75{ \mathrm{CH}}_4$	Equation (7)	−203	−23

^a The ranking is determined by thermodynamic principles and tends to show the sequence of these reactions. The reaction which yields the highest energy occurs first, followed by the less yielding energy reactions [77]. Based on this concept, methane oxidation is the most thermodynamically favorable chemical reaction and is the primary source of biogenic heat at NSZD sites (see Section 4.4). ^b Changes in free energy in fact show change in entropy and heat exchange is resulted from changes in enthalpy. The values are calculated by the half-reaction bioenergetic values assembled by Stockwell [78] (refer to Table 1 in that reference). Note that changes in input parameters (i.e., enthalpies of products and reactants) can alter the overall enthalpy of the reaction. ^c Convert kJ/mole to kJ/g using decane molar weight of 142.29 g/mole. ^d For methane oxidation, the free energy and heat exchange are expressed in decane equivalent, to provide unit consistency, comparability of the values and convenience in comparison. In calculations, the denominator was on decane basis (i.e., kJ/mole decane), instead of methane.

), methanogenesis is the least thermodynamically favorable because as demonstrated by ΔG_r^o values in Table 1, it yields least energy for cell growth. Due to limited concentrations of electron acceptors in most contaminated zones (see Table 2.3 in Wiedemeier et al. [79]), reductive pathways associated with these electron acceptors most likely leave a significant portion of LNAPL unprocessed for methanogenesis. It is anticipated that methanogenesis is responsible for depletion of up to two orders of magnitude greater PHC degradation than it is typically estimated by monitored natural attenuation (MNA) models (e.g., BIOSCREEN model [80]). The MNA is the predecessor of NSZD and only considers electron-acceptor-mediated degradation of dissolved BTEX [25].

If we base the discussion on thermodynamic principles, methanogenesis can theoretically begin when all electron acceptors (i.e., oxygen, nitrate, manganese, iron and sulfate) are depleted [78], which might make the impression that methanogenesis can only occur in mature oil spill sites. However, recent studies showed that syntrophic PHC degradation can proceed even in the presence of electron acceptors such as nitrate, Fe(III) and sulfate [36,82,83,84]. Using most probable number (MPN) technique, Bekins et al. [69] determined population distributions of six physiologic types (i.e., aerobes, denitrifiers, iron-reducers, heterotrophic fermenters, sulfate-reducers and methanogens) in the anerobic portion of an LNAPL plume. Culturable methanogenic and iron-reducing bacteria were found together, albeit with a strong inverse relation between their populations. These observations anticipated that the methanogenic condition evolved transitioning from iron-reduction to methanogenic condition through more than 20 years of contamination [69]. Finding iron-reducing and methanogenic organisms usually together [69,85,86] suggests that the associated processes can be active simultaneously [69,87].

2.4.3. Microbial Community Structure and Syntropy

The methanogenic process is not as straightforward as Equation (7) (see Section 2.4.2). In anaerobic groundwater aquifers, like any other anaerobic ecosystem, syntropy exists to ensure microbial consortia perform as a single catalytic unit for the community’s common benefit [36]. According to metagenomic and metaproteomic analyses, Pelotomaculum first degrades terephthalate via decarboxylation, producing intermediates such as acetate and butyrate. Secondary syntrophs (Candidatus Cloacamonas acidaminovorans affiliated with Thermotogae and OP5) oxidize butyrate to acetate and H₂. Butyrate oxidizers can also be fed by Caldiserica (WWE1) which fixes H₂/CO₂ to butyrate. Thereafter, methanogens such as Methanosaeta and Methanolinea produce CH₄ and water from processing H₂ and CO₂ [88,89]. Although terephthalate is not a hydrocarbon, its methanogenic degradation can be used as a model to shed light on possible PHC contamination pathways in a contaminated subsurface [36].

Studying soil core at 15-cm intervals, Irianni-Renno et al. [41] demonstrated that LNAPL degradation was mediated by the syntrophic relation of fermenters and methanogens in Region III and methanogens and sulfate-reducing bacteria in Region IV (refer to Figure 1 to see these Regions). This supports Bekins et al. [69] who observed that methanogenic and iron-reducing organisms coexist in anaerobic LNAPL plume (see Section 2.4.2). It posited that anaerobic benzene degradation completes through unusual, complex syntropies formed between the microorganisms present in iron-reducing enrichment cultures [86,90], such as Peptococcaceae and Desulfobulbaceae [86] and Cryptanaerobacter and Pelotomaculum phylotypes [91]. Benzene degradation can happen under denitrifying conditions in enrichment dominated by syntrophic Peptococacceae and Betaproteobacteria species [92]. Gram-positive Peptococcaceae species are closely related to clones globally recovered from contaminated aquifers [86].

Syntropy can exist at deeper levels, i.e., among methanogens. Syntrophaceae spp. (e.g., Syntrophus and Smithella) are thought of as key players in syntrophic hydrocarbon metabolism, as they are frequently found in PHC impacted environments [93,94]. Molecular analysis of PHC plume at a former refinery in Carson City, Michigan, indicated that 21% and 32% of the microbial communities which were respectively isolated from free-phase and dissolved-phase contaminated zones had over 97% similarity to Syntrophus sp. strain B2 from Deltaproteobacteria subdivision [95]. In a shallow PHC contaminated aquifer, Syntrophus- and Methanosaeta-related species were respectively dominant bacterial and archaeal organisms below the water table [41]. Members of the Syntrophus genus bridge hydrocarbon fermentation and hydrogenotrophic methanogenesis through converting fermentation by-products (propionate, acetate and butyrate) to CO₂, H₂ and formate [95]. Callaghan et al. [96] studied Desulfatibacillum alkenivorans AK-01 as a model organism for anaerobic alkane biodegradation due to its versatile metabolism (e.g., sulfate reduction and n-alkane biodegradation via fumarate). In the lack of sulfate, this strain could produce methane from n-hexadecane when cocultured with H₂/formate-using methanogens. Its genetic signatures suggested that interspecies formate transfer might play a critical role in the syntrophic pairing of this strain and other species [96]. Syntrophus buswellii and Syntrophus gentianae are characterized as bacteria catalyzing syntrophic oxidation of aromatic compounds (e.g., benzoate, crotonate, gentisate, hydroquinone) via interspecies hydrogen transfer [57]. In a binary mixed culture composed of Syntrophus gentianae and Methanospirillum hungatei, the former almost completely converted benzoate to acetate, methane and CO₂, while the latter scavenged on hydrogen [57].

Microbial evidences suggest that temporal and spatial fluctuations of NSZD rates (refer to Garg et al. [25] and Sihota et al. [26] to see variable NSZD rates) are not necessarily rooted in the limitations or “inaccuracy” of NSZD rate quantification techniques, but are either inherent or caused by environmental conditions (see Section 2.3 and Section 3). Multiple methanogenic pathways may occur in an LNAPL impacted subsurface, with temporospatial preference for one over another due to many factors which cause variable microbial NSZD degradation rates over time and space. For instance, it is proposed that acetate depletion in PHC contaminated environments is mostly due to acetate oxidation coupled with hydrogenotrophic methanogenesis (reactions 5 and 6 in Figure 1), rather than acetotrophic methanogenesis (reaction 4 in Figure 1) [89,93,97,98]. High CO₂ concentrations can disrupt this syntropy and shift the metabolic reactions towards acetotrophic methanogenesis [99]. Extrapolating from the dynamics of redox reactions in rumen fermentation dictated by the balance between hydrogen ions and H₂ [100], it is thought that high H₂ partial pressure can impose similar effects on fermentation pathways as high CO₂. A variable CO₂ concentration of 6–14 v/v% (the denominator is defined as overall vadose zone gases) has been observed in soil column over LNAPL body [26], which can be one of the many reasons for spatial heterogeneity of microbial community observed in LNAPL impacted areas. Temperature is another factor altering syntropic microbial relationships at NSZD sites. In a laboratory study, molecular analysis demonstrated that not only microbial composition but also microbial activity (i.e., biodegradation rates) shifted with temperature. One-year exposure to temperatures 4 and 9 °C ceased biogas generation and sulfate reduction took over. At temperatures 22 °C and 30 °C, dominant archaeal genera included Methanosaeta, Methanobacterium and Methanosarcina. Based on the presence of these species, it was inferred that acetoclastic methanogenesis was the main methane generation mechanism [64], with the possibility of hydrogenotrophic methanogenesis by Methanosarcina [101].

These implications are further supported by the findings of Allen et al. [95], where clone libraries of 16S rRNA gene sequences retrieved from four depths of a contaminated field site and two depths of an uncontaminated background site were examined. Correspondence analysis proposed vertical microbial gradients, composed of various aromatic hydrocarbon-degraders (e.g., Sphingomonas aromaticivorans and Brachymonas petroleovorans); syntrophic hydrocarbon-degrading, sulfate-reducing and dissimilatory iron-reducing species; large populations of methylotrophs and methanotrophs; and fermenting bacteria. This diverse microbial community contrasted with the background’s common soil microbiology (e.g., Acidobacteria and Actinobacteria). On a side note, this change in microbial structure was concomitant with a spectrum of geoelectrical properties at the examined area, suggesting that the evolution of microbial populations in a subsurface imposes strong effects on their geophysical surroundings. Multivariate statistical analysis revealed that a specific microbial culture composed of syntrophic δ- and ε-Proteobacteria caused elevated bulk electrical conductivity. This observation pinpointed that geoelectrical measurements may be used as a cost-effective monitoring tool for investigating the microbial ecology of remediation sites [95].

3. NSZD Rate Quantification Techniques

So far, four principal methods are developed to resolve NSZD rate, including concentration gradient (CG), dynamic closed chamber (DCC) and CO₂ trap and thermal monitoring. The first three methods use either at-grade reactant gas fluxes (e.g., O₂ [11]) or product gas fluxes (i.e., CO₂) associated with LNAPL degradation as a surrogate for NSZD rate measurement. Thermal monitoring adopts biogenic heat as an NSZD signature.

The CG method is based on Fick’s law of diffusion and quantifies LNAPL mass loss due to volatilization and biodegradation. Nested gas probes in a soil column or multi-level samplers are used to develop vertical soil gas concentration profiles in the subsurface [23,29,47]. The DCC measures surficial CO₂ effluxes and the associated device consists of a soil gas chamber which is placed at grade and seals over a soil collar. It is connected to an infrared gas analyzer (IRGA) which measures CO₂ concentration accumulated in the chamber. An increase in CO₂ in the chamber’s headspace is concomitant with an increase in absorption of infrared light and a decrease in light transmission which is correlated to CO₂ accumulation. Using the in situ readings, CO₂ efflux can be quantified [102]. The device and principles of CO₂ trap and thermal monitoring are thoroughly discussed in detail in Section 4 and Section 5, respectively. The discussion about these two techniques is expanded due to particular interest in them in the current remediation and reclamation industry.

Advantages, limitations and niche applications of these methods are presented in a matrix presented in

Table 2. Summary of advantages and limitations of NSZD rate quantification techniques.. Appendix B discusses each technique individually and provides a thorough literature summary of each technique’s advantages and limitations.

Technique	Advantages	Limitations	Niche Applications
Concentration gradient (CG)	Gives instantaneous measurement. No need for background correction (depending on the installation depth of the sampler). Can resolve NSZD rate via different physicochemical pathways. Can determine soil properties and texture in a soil column.	Intrusive. Time consuming and labor intensive. Only provides a snapshot of subsurface gas profiles. Quantifies diffusive gas transport only, overlooking the effect of advection. Accuracy can be affected by many environmental conditions. Accuracy can be affected by anthropogenic interferences.	Capped sites (e.g., covered landfills).
Dynamic closed chamber (DCC)	Gives instantaneous measurement. Short runtime (i.e., surficial CO2 efflux measurement within 90–120 s). Can take multiple measurements in a short period of time, so can also be used for long-term measurements. Minimally intrusive. Provides real-time field values. Easy to transport device. Accounts for both diffusive and advective gas transport mechanisms.	Time consuming and labor intensive (less than the CG though). Correction for natural soil respiration (NSR) is required, using either stable carbon and radiocarbon isotope analysis or background correction. Stable carbon and radiocarbon isotope analyses require sample collection, expensive analytical procedure and trained personnel for the analysis. Finding a background location which can be representative of the NSR of NSZD spot is a challenge. Accounts for NSZD rate via biodegradation only. Accuracy can be affected by many environmental conditions. Accuracy can be affected by anthropogenic interferences. Estimate NSZD rate via biodegradation only. Cannot be used throughout the entire year in cold regions like Canada. It is highly likely that it captures NSZD rates with a few months’ delay.	Site-wide surveys (e.g., screening LNAPL impacted areas).
CO2 trap	Provides time-integrated, long-term measurement of surficial CO2 effluxes. Simple. Cost-effective. Passive (minimal labor requirement). Long-term measurement can offset the adverse effect of barometric pumping. Minimally intrusive. Accounts for both diffusive and advective gas transport mechanisms.	Time consuming. Detection limit depends on the device’s opening. Accuracy can be affected by many environmental conditions. Accuracy can be affected by anthropogenic interferences. Estimate NSZD rate via biodegradation only. Cannot be used throughout the entire year in cold regions like Canada. It is highly likely that it captures NSZD rates with a few months’ delay.	Long-term monitoring of NSZD sites.
Thermal monitoring	Short-term and long-term measurements, depending on the type of biogenic heat sampler utilized. Cost-effective. Yearlong measurement. High-resolution data collection. Averaging out negative and positive biases. Converging apparent and actual NSZD rates. Higher accuracy (compared to other methods).	Intrusive. Labor intensive for field installations and desktop analysis (data reduction, programming, mathematical skills). Requires background correction, which is a challenge to spot the correct location. Accuracy of this method is highly sensitive to environmental factors that can either suppress or cease methanotrophic activity. Estimate NSZD rate via biodegradation only.	Long-term monitoring of NSZD sites, especially at regions with natural (e.g., ice cover) and artificial (e.g., urban developments with pavements) impervious surfaces.

and the supporting literature information is discussed in Appendix B.

4. Passive CO₂ Trap: A Cost-Effective Technique

Due to simplicity and cost-effectiveness, use of CO₂ trap for NSZD rate quantification is economically advantageous over other techniques presented in Section 3. Additionally, this method allows board-scale, replicated measurements for resolving not only NSZD rate but also natural ecosystem carbon exchange at a global scale [30]. As such, it is worth giving additional attention to details associated with its usage, such as device configuration, method accuracy and deployment conditions.

4.1. Components and Configurations

The main component of CO₂ trap is an alkali chemical that CO₂ is absorbed into. The predecessor of this method, which was mostly used for studying soil CO₂ evolution in agricultural ecosystems, was equipped with a caustic solution, such as NaOH/BaCl₂ and KOH (see [103] and the references therein). Porous alkaline solid (i.e., soda lime) has been used as a replacement for caustic solution [104,105], to eliminate handling, freezing and spillage of the liquid [105,106].

Soda lime is a variable mixture of sodium hydroxide (NaOH), calcium hydroxide (Ca(OH)₂), water (~20 wt%) and in some cases potassium hydroxide (KOH) which forms (sodium, potassium and calcium) carbonates when reacting with CO₂ [104]. Water is a crucial constituent [107] because, as shown in Equations (8)–(10), it initiates CO₂ absorption reactions to form the different forms of carbonate [30].

$${\mathrm{CO}}_{2\left(\mathrm{g}\right)}+{ \mathrm{H}}_2\mathrm{O} \to { \mathrm{H}}_2{\mathrm{CO}}_3$$

(8)

$$2{\mathrm{H}}_2{\mathrm{CO}}_3+2\mathrm{NaOH}+2\mathrm{KOH} \leftrightarrow { \mathrm{Na}}_2{\mathrm{CO}}_3+{ \mathrm{K}}_2{\mathrm{CO}}_3+4{\mathrm{H}}_2\mathrm{O}$$

(9)

$$2\mathrm{Ca}{\left(\mathrm{OH}\right)}_2+{ \mathrm{Na}}_2{\mathrm{CO}}_3+{ \mathrm{K}}_2{\mathrm{CO}}_3 \leftrightarrow 2{\mathrm{CaCO}}_{3\left(\mathrm{s}\right)}+2\mathrm{NaOH}+2\mathrm{KOH}$$

(10)

Soda lime should be placed in a confined space, like beneath a chamber, to capture gas emissions at grade. Two main configurations of CO₂ trap used for NSZD rate estimation are proposed by McCoy [56] and Keith and Wong [30] (Figure 4). The McCoy model resulted in a patent in 2014 [108], which has been commercialized by E-Flux [109] and its performance has been compared to other NSZD rate quantification techniques [11,46,47]. This model is a cylindrical device composed of a Schedule 40 polyvinyl chloride (PVC), with internal diameter (ID) of 10 cm, surface area of 0.008 m², height of 30 cm and rubber O-rings to create air-tight seals [56]. It is featured with two novel characteristics, incorporated into the design patented by Zimbron et al. [108]: (1) an unrestricted advective flow-through (in simpler terms, an open-top design); (2) top and bottom trap elements. It is supposed that the open-top design prevents interferences of concentration and/or pressure build up in the sealed chamber [56]. Two soda lime absorbent elements mounted at two cross sections of the pipe and each containing 50.0 g of granular soda lime convert gaseous CO₂ to solid carbonate. While the bottom element absorbs the soil efflux, the top element prevents atmospheric “contamination” by intercepting air CO₂. The trap is placed on a PVC receiver (typically 10-cm ID, or 20-cm ID for added sensitivity as it can increase the covered surface area [22]) which penetrates the soil a few centimeters (≤18 cm). To protect against precipitation, a cover made of a larger, perforated pipe (e.g., thin-wall PVC pipe with an ID of 15 cm) is integrated into the trap design [47].

The Keith and Wong [30] model is simpler than McCoy’s design. It is based on the configuration previously used in forestry studies [106]. This model is composed of a large, bowl-shape chamber, with dimensions shown in Figure 4b and a surface area of 0.08 m². However, based on the researcher’s recommendations, dimensions are flexible as long as the chamber size is large enough to capture soil spatial heterogeneity. Granular soda lime (50 g, 2.4–4.8 mm diameter per granule, mesh size of 4–8) is exposed to CO₂ efflux by placing it in a sufficiently wide, inert (e.g., glass) dish which provides a surface density of 0.06 g soda lime per cm² of soil surface in chamber. A white color should be opted for the trap body, to reflect sunlight and minimize overheating. Reflective coating or a shield can be used instead as a more effective strategy. The chamber should penetrate the soil to 1 cm [30].

Both designs should be optimized to improve the method reliability. Based on preliminary experiments, Tracy [47] provided recommendations for optimizing the McCoy model, about sorbent mass, CO₂ breakthrough from the bottom to the top element and height of the bottom sorbent element from the soil surface and insertion depth [47]. Method verification with respect to different influencing parameters was performed on the model developed by Keith and Wong [30]. However, in both cases, conclusive results were not achieved to optimize the design properly and comprehensively. Additionally, relative accuracy of these configurations is not known. This warrants further research, especially given that the Keith and Wong model has a closed top and a 10-time larger surface area than that of the McCoy model. These differences can impose negative or positive influence on precision and accuracy of CO₂ efflux measurements (see Section 4.2). Method verification and optimization are a must, considering that many factors which are discussed in Section 2.3 influence surface fluxes. Chamber design and installation can cause deviating from natural environmental conditions.

4.2. Accuracy and Precision

Many researchers attempted to determine the precision of CO₂ trap, equipped with either caustic solution [103,110,111,112] or soda lime [11,46,47,113], relative to other NSZD quantification techniques, such as DCC [30,47,103,110,111,112,114], thermal monitoring [11,46,113], CG [47,114] and absorbent Solvita gel [112], which is a patented environmental measurement system and is sensitive to certain gases [112].

Linear regression analysis and coefficient of determination (R²) indicated that the data of alkali CO₂ trap collected from different literature sources agreed well with those obtained by other methods, especially DCC (R² = 93%) and CG (R² = 94%) (Figure 5). R² for CO₂ trap versus thermal monitoring varied 11–51% (Figure 5c). A low R² value of 11% was due to a large discrepancy between December 2014 data points. At this time of the year, traps were unable to capture CO₂ efflux most likely due to frozen surfaces. Upon removing December 2014 data points, R² value bounced to 51%. In all the cases, regression curves diverged from 1:1 line, with a tendency toward the trap data (i.e., x-axis), indicating that trap tends to overestimate NSZD rate compared to the other techniques. However, this finding should be validated with more data points, especially with those concentrated at higher CO₂ effluxes and that these supplementary data points (and generally any precision study) should be collected systematically at a given time and location to improve the comparability of the results.

However, studying the literature sources individually revealed contradictions. At a geologically heterogeneous site, the most consistent NSZD data was collected using carbon traps and thermal method [11]. Apparent NSZD rates estimated using these two techniques showed a difference of 8.4%, indicating validity of both methods [46]. Comparing DCC and trap, R² was 0.61–0.82 for field data [30,110,111] and 95% under controlled laboratory conditions [112]. In Rochette et al. [110], a negative exponential relation (R² = 0.61−0.74, depending on the soil texture) was established between CO₂ fluxes captured by DCC and trap. Paired t-test indicated that the values captured by alkali traps were on average 0.008–0.016 mg CO₂/m²/s higher than those determined using DCC, which was occasionally statistically insignificant (p > 0.05) [111]. On the other hand, Freijer and Bouten [103] reported an underestimation of CO₂ efflux by alkali traps and explained it by low absorption velocity in traps, dictated by slow diffusion due to lack of proper contact between the gas and sorbent. In Rochette et al. [111], trap occasionally underestimated CO₂ efflux by about 22%, compared to DCC. Rochette et al. [110,111] identified that changes in soil temperature and CO₂ concentration inside trap did not establish any relation with lower respiration captured by alkali traps, compared to DCC.

Precision does not necessarily guarantee accuracy [115]. All these techniques can introduce bias into NSZD rate measurements. This argument can be supported by the findings of a controlled lab study performed by Tracy [47], where DCC and alkali trap almost accurately captured the imposed flux (with an error of ±7%). Whereas, CG underestimated the imposed flux within 38% (see Figure 6). It is because almost all these methods are more or less disruptive and their measurements are affected by different factors. As discussed earlier (see Section 2.3), these factors can affect microbial processes which advance the NSZD and/or gas transport mechanisms in the subsurface. Tracy [47] explains that underestimation of the CG method (which can be understood from Figure 6b) was due to: (i) the CG results largely depending on accurate estimation of effective diffusion coefficients; (ii) the CG method only incorporating diffusive gas transport, overlooking the effect of advection [47]. In the field study of Eichert et al. [114], CG showed the highest sensitivity to input parameters compared to CO₂ trap and DCC, while CO₂ trap indicated relatively high intersample variability [114]. It is argued that DCC, thermal monitoring and CO₂ trap are likely similarly biased by water table fluctuations and soil moisture. However, the available CO₂ trap data in the literature does not necessarily reflect this, partly due to scarcity of CO₂ trap data [46].

It is proposed that continuous monitoring of NSZD rate (e.g., at least a year) can average out positive and negative biases and converge the apparent and actual NSZD rates [46]. This approach can specifically become important in using CO₂ traps, as NSZD rates and gas transport processes substantially vary seasonally. Seasonal variations resulted in a lag time of 5–7 months between peak NSZD-driven subsurface CO₂ concentration (happened in fall–winter) and peak CO₂ efflux at grade (observed in summer). In cold regions like Canada, long-term periods of frozen soil can inhibit gas emissions and exacerbate the situation [26].

4.3. Field Deployment and Grid Design

Use of CO₂ trap for quantifying NSZD rate (as opposed to its usage in agriculture and forestry research) is quite recent [11,30,37,46,47,56,108]. In the last two decades (from 2005 to present), there have been numerous studies about NSZD concept and mechanisms [25,27,28,29,30,31,32,33,34,35,36,37,38]. Based on these studies, Interstate Technology and Regulatory Council (ITRC) has updated its guideline with NSZD, where CO₂ trap is incorporated as one of the frequently-used NSZD assessment methods [22]. The ITRC is a state-led coalition in the U.S. that works closely with the legislative bodies (e.g., U.S. Environmental Protection Agency (US EPA)) and other institutions (e.g., Environmental Council of the States (ECOS)) to reduce barriers in adopting innovative air, water, waste and remediation environmental technologies and processes.

4.3.1. Field Deployment

The main part of trap body or its receiver (depending on the design; see Section 4.1) should be installed prior to trap deployment (i.e., commencing CO₂ absorption into sorbent). For the Keith and Wong model, it is recommended that CO₂ trap is installed and left uncovered several weeks before starting the measurements [30], while for the McCoy model, this time has been limited to a few days (≥1 day) [56]. Regardless of this discrepancy in the recommended duration, this procedure is to offset adverse effects of field installation, by allowing soil to recover from disturbances. The mechanical disturbance of soil may alter the diffusion coefficient and can also cause a flush of CO₂ due to provoked microbial activity [30,56]. During installation, care should be taken to have minimal deviation from natural conditions. It is recommended that the soil within and around the receiver (with the McCoy model) is recompacted back to its original level [22,56]. This strategy can be adopted when the Keith and Wong model is used, which includes recompacting soil around the edges of trap body penetrated 1 cm into the soil (see Section 4.1).

4.3.2. Deployment Grid Design: Number and Distribution Pattern

The number and distribution pattern of CO₂ traps at NSZD fields (hereinafter, the deployment grid) are not defined to this date. In the field studies conducted so far, CO₂ traps seemed to be randomly distributed over the anticipated LNAPL area, with almost no clarification about the decision-making process and considerations for the development grid design [11,37,46,56]. In a multiple-site study (six sites), 10–23 traps were deployed per site, with the traps placed at 1 to 3 m distance from each other and a few traps (2–3 traps) reserved for background location (non-impacted areas) [56]. As part of this review, we performed an area estimation analysis on the examined sites in McCoy [56] using ImageJ 1.52 v. The LNAPL impacted areas varied in a wide range of 6.2–188.2 ha (15.3–465.2 acres). This area was 28.5, 188.2 and 66.6 ha at those sites where similar number of traps (i.e., 23 traps) were deployed (Sites A, E, F, respectively in McCoy [56]), indicating that the area covered by each trap was not chosen equally. In each site, there was at least one spot with replications (e.g., at Site A, three traps were placed 2 m apart) to capture CO₂ efflux variability at a single spot [37,56]. In other studies, a limited number of traps (i.e., three to five CO₂ traps) were only deployed at the LNAPL zones, with sporadic trap distribution pattern and no replication. Limited information was provided about LNAPL distribution in the subsurface and the impacted area [11,46]. Thus, an urgent need exists to develop a grid design protocol for this technique if the method is to be used frequently and is expected to provide accurate representation of the field conditions.

In the relevant studies [11,37,46,56], we anticipate that the researchers have most likely considered a few factors to define the grid design: (i) historical LNAPL contamination data (gas and liquid analyses); (ii) anticipated shape of LNAPL body; (iii) relative distance in each pair of deployed traps (1–3 m). These have not been explicitly confirmed in those studies. In these sources, relevant information about the deployment grid is merely an informative piece of methodology, rather than a detailed description of a systematic experimental design for determining/optimizing the deployment grid design [11,37,46,56]. Hence, it is not feasible to develop a grid design protocol based on the information in the current literature. It is hypothesized that the density of CO₂ traps (i.e., the number of traps per unit LNAPL impacted area) and the grid pattern (the relative distance and placement of the traps) can alter the obtained NSZD rates. Additionally, the required number of CO₂ traps might depend on not only distribution of the LNAPL body (the magnitude and shape of LNAPL area) but also the LNAPL composition (e.g., light vs. heavy LNAPL). Therefore, the outstanding question is “what an optimized grid design protocol is to resolve NSZD rate most accurately and precisely”. Further research on this topic is critical for amending the current guidelines to include the most comprehensive information about this site assessment technique and a consistent protocol for the deployment grid design. The outcome of such research will help practitioners to adopt this technique with utmost certainty and provides researchers with a consistent protocol which will ultimately increase the comparability of the results obtained in different trap studies.

4.3.3. Preliminary Site Delineation

CO₂ traps should be deployed over the LNAPL affected zone to capture CO₂ efflux originated from the contamination source, leading to a necessity for the spatial extent of LNAPL body to be known. Site characterization and LNAPL mapping, also referred to as subsurface investigation or site delineation, can be performed through analyzing historical data. Different data collection methods exist, which are generally referred to as direct push technologies [116]. These methods include: (i) rapid optical screening tool (ROST) laser-induced fluorescence (LIF) system (in short, the LIF) [117,118,119,120,121,122]; (ii) geophysical, non-destructive techniques, such as electrical resistivity (EC), electromagnetic induction and ground penetrating radar (GPR), which all take advantage of perturbations in subsurface physical characteristics (e.g., EC and permittivity) due to the presence of LNAPL—the LNAPL act as an insulator (low EC of ~0.001 μS/m, as opposed to ~0.1–1 μS/m of formation fluids) and have low permittivity (2–3) compared to clean water (~80) and water-saturated sands (20–30) [10,123]; (iii) soil vapor surveys [29,31] which can detect the hotspots of LNAPL biodegradations via in situ CO₂, CH₄, organic vapor readings; (iv) soil borings and soil sample collection [124] which can reveal contaminant composition and distribution and microbial community structure in a soil column; (v) installation of groundwater monitoring wells, followed by groundwater (and possibly gas) sample collection and analysis [125].

Site delineation in the majority of NSZD studies was performed using LIF, either solely or in combination with soil core and monitoring well data [37,46,56,64]. Adopting LIF technique, McCoy refined the extent of LNAPL body which was primarily approximated using other methods. When only limited monitoring well data was utilized for site delineation [11,56], a “blob shape” could only be distinguished for LNAPL impacted areas [56].

The LIF is adopted not only for determining NAPL mobility (migration and distribution) [118], but also for volume quantification (saturation of NAPL in soil porous structure) [117]. The LIF response has shown to correlate to LNAPL transmissivity (T_n) up to a greater extent than LNAPL saturation (S_n) [118]. Teramoto et al. [117] ascertained that relation between LIF response and LNAPL saturation followed an exponential function, where the empirical coefficients highly depended on the composition and the weathering stage of LNAPL, suggesting that these coefficients were site-specific [117]. The traditional site delineation techniques may not quantify actual NAPL volume appropriately, which can in turn result in false interpretation of the effectiveness of remediation techniques [117], such as NSZD. Use of LIF is a step forward in effective site management, providing accurate LNAPL volume and distribution, even at sites with complex commingled NAPL contamination [117,118].

Contamination with chlorinated solvents is often classified as dense non-aqueous phase liquid (DNAPL) [126]. Presence of DNAPL in groundwater is associated with historical release due to storage or use of these chemicals for metals degreasing, dry cleaning and other uses. These solvents are denser than water and often penetrate deep into permeable soil. They penetrate deep, sinking to the bottom of the water body. When capillary pressure does not allow further penetration, they move laterally along the top of finer-grained soil layers. Under the water table, they slowly dissolve in water. Due to aging, DNAPL dissolution reduces over time, decreasing DNAPL saturations and the likely disappearance of part of the residual DNAPL in groundwater [127]. The term commingled NAPL is also used to refer to DNAPL contamination. Commingled NAPL is composed of LNAPL detected under the water table along with solvent DNAPL constituents, where it is atypical to detect concentrated LNAPL. Commingled NAPL is very common as most halogenated solvent discharges are not pure. These releases mainly consist of trichloroethylene/tetrachloroethylene (TCE/PCE) combined with polyaromatic petroleum constituents released into the commingled NAPL due to industrial processes, disposal and/or transport [119].

A new variant of LIF “DyeLIF^TM” can provide high-resolution, three-dimensional subsurface mapping of all types of NAPL (LNAPL, DNAPL and commingled NAPL), including chlorinated solvents and other dense, non-naturally fluorescing NAPL [120,127]. The DyeLIF^TM method showed excellent (98%) agreement with the results obtained through aboveground soil core analysis. This method could detect DNAPL even in samples with relatively low DNAPL pore saturation (e.g., 0.7%) [127].

4.4. Downstream Laboratory Analyses and Result Interpretation

4.4.1. Recovery and Quantification of Absorbed CO₂

CO₂ absorption is proportional to weight gain of the sorbent [30,105]. Acid-base titration is utilized for estimating CO₂ absorption in caustic solutions [103] and dry weight change analysis has been mainly adopted when soda lime is used [30,104,105]. However, in most literature sources focused on the use of CO₂ traps for NSZD rate resolution, the concept of titration has been incorporated into weight-change analysis [37,47,56,108]. A strong acid, like HCl, has been utilized to convert solid CO_2(s) (carbonate) to gaseous CO_2(aq). This is an adoption of “gravimetric acid analysis” developed by Bauer et al. [128] for determination of carbonates in soil samples. A weight difference determines the amount of CO₂ released through acid recovery [128].

The weight-change method with no acid addition is simpler and safer, compared to gravimetric acid analysis. However, it has the risk of exposure to atmospheric CO₂ contamination during drying and weighing. Drying is recommended to be performed at 104 °C for 14 h in a fan-forced oven before and after deployment [30]. For shorter drying duration, we recommend that weight measurement is repeated every 1–2 h until no weight change is observed between two consecutive weighings. It is hypothesized that soda lime’s water content completely evaporates within the first few hours. Based on our lab experience (unpublished data), drying can be performed in a microwave oven within 3–4 min only, which eliminates the risk of exposure to atmospheric CO₂. Fan-forced drying can particularly increase the possibility of sample contamination in the first few hours before water content is substantially removed. As illustrated in Equation (1), water initiates the cascade of reactions required for CO₂ absorption to soda lime. Correct understanding of the error associated with the laboratory procedure of dry weight change analysis requires further quality assurance/quality control (QA/QC) analysis.

Dry weight change analysis requires an important correction. According to Equation (4), which is the sum of Equations (1)–(3), the absorption process generates 1 mole water (H₂O) per 1 mole CO₂ absorbed, which was originally bound to soda lime through chemical bonds in the form of two hydroxyl (-OH) groups [104]. Evaporation of this water during the drying process causes underestimation of CO₂ efflux. To account for this error, a correction factor of 1.69 is recommended to be applied on the weight difference [104].

$${\mathrm{CO}}_{2\left(\mathrm{g}\right)}+\mathrm{Ca}{\left(\mathrm{OH}\right)}_2 \to \mathrm{CaC}{\mathrm{O}}_{3\left(\mathrm{s}\right)}+{\mathrm{H}}_2\mathrm{O}$$

(11)

Field, travel and laboratory blanks should be incorporated into the protocol, to respectively correct for any error during field deployment/retrieving, transfer process and downstream laboratory analyses [30,104,108].

Thus far, there is no study to validate the relative accuracy of these two analytical procedures (i.e., dry weight-change method and gravimetric acid analysis). A brief literature review of Freijer and Bouten [103] demonstrated that mass balance error resulted from using acid titration was −9.1% to −0.6% [129,130,131], while it was −33% with weight change (without acid addition) technique [132]. However, in this source, weight change analysis was applied to NaOH fixed on carriers [132], making it difficult to extrapolate this observation to the weight change analysis of soda lime.

4.4.2. Origin of Absorbed CO₂

As depicted in Figure 1, a portion of emitted CO₂ at grade is associated with modern carbon (modern-C) which can lead to overestimation of NSZD rate if not accounted for appropriately. Modern-C is due to respiration of plant roots, rhizosphere, microbes and fauna residing in the surficial soil zone [34,133]. To distinguish between the carbon emitted from LNAPL biodegradation (ancient-C) and modern-C, background subtraction, isotopic analysis or both should be incorporated into alkaline trap technique [34,35,38,56].

Background subtraction is based on the assumption that total soil respiration (TSR) is the sum of contaminant-driven respiration (CSR) and natural soil respiration (NSR) [34]. Background CO₂ efflux measurements are collected by traps deployed over “unaffected” areas. Thus far, background correction is the simplest and least expensive approach and is most suitable for sites where CSR >> NSR, either due to high CSR rates or low NSR rates. Otherwise, a similarity between CO₂ effluxes of background locations and LNAPL areas can be seen which can lead to ambiguity about the occurrence of NSZD at sites with relatively low LNAPL loss rates. In addition, spatial variability between background readings (due to soil heterogeneity and/or variable NSR) can result in additional uncertainty in determining NSZD rates [35,56]. Moreover, due to uncertainty about the mobility and areal extend of LNAPL in many active and decommissioned industrial sites, selecting a true background location poses a great challenge [56].

Analyzing stable (δ¹²C and δ¹³C) and radioisotope (¹⁴C) carbon signatures can significantly improve the method accuracy through identifying the origin of captured CO₂ most appropriately. Prior to its application to NSZD rate resolution [34,35,38,56], isotopic analysis was used to: (i) monitor natural attenuation (MNA) at hydrocarbon and chlorinated solvent contaminated sites [133,134,135]; (ii) study weathering of petroleum reservoirs [136]; (iii) determine contribution of anthropogenic and natural sources to atmospheric gases, including CO and CO₂ [137,138,139]. In MNA, changes in δ¹³C of CO₂ could indicate the dominant microbial (e.g., methanogenesis vs. aerobic methane oxidation) and metabolic (e.g., acetotrophic vs. hydrogenotrophic methanogenesis) pathways [35,133,135,140]. Differences in δ¹³C content of organic matter species are used as a tracer to understand the contribution of PHC to metabolic end products found at MNA sites [135].

Nevertheless, interpretation of stable isotopic data is often complicated due to, most likely than not, simultaneous influence of multiple reaction pathways [35]. Using stable carbon isotope technique requires knowing fractionation factors for the possible metabolic pathways (e.g., simultaneous conversion of CO₂/H₂ and acetate to CH₄ in methanogenesis). These factors are not constant and sometimes largely differ from site to site and condition to condition. This effect highlights that site-specific circumstances can greatly impact the structure of natural microbiomes and that these factors should be determined for each individual site independently [141]. In addition, isotopic fractionation during acetate turnover poses an additional challenge as acetate can enter multiple competing reactions as either reactant or product [141]. Acetate can be produced from organic matter (i.e., heterotrophic acetogenesis [36]) and CO₂ (i.e., autotrophic acetogenesis [141,142]). It can be utilized for biomass generation as well as CH₄ and/or CO₂ production [36,141]. Another issue associated with stable isotope technique is that a significant overlap exists between δ¹³C contents of PHC, some indigenous soil carbon sources (e.g., plant root respiration and soil carbonates) and metabolites of methanogenesis [134,135]. In biochemical processes, the ratio of stable isotope of metabolic end-products per substrate materials changes. This fractionation effect is much smaller in aerobic metabolism (<6‰) than in anaerobic processes. Even in anaerobic processes, this ratio can be significantly different considering different by-products. For instance, δ¹³C of CO₂ and CH₄ produced through acetotrophic methanogenesis is respectively 20–30% higher and 20–30% lower than the substrate, i.e., acetate (see Conrad et al. [135] and the references therein). Thus, the use of stable isotope δ¹³C technique without supporting information can produce ambiguous results. For example, high δ¹³C CO₂ produced through methanogenesis can be mistakenly taken for high δ¹³C atmospheric CO₂. As such, radioisotope ¹⁴C technique is used as a complementary analysis [135].

Use of radiocarbon as a tracer relies on the fact that modern-C is rich in ¹⁴C. Whereas, ancient-C is entirely depleted of ¹⁴C due to long-term weathering. This results in a distinguishing feature between carbons with different ages (in this case, biogenic vs. fossil sources) [11,35,37,56,134,135,138,140,143]. Radiocarbon is a cosmogenically produced radionuclide, with a half-life of 5730 ± 40 yr. It is produced in the upper layers of the atmosphere (troposphere and stratosphere) through bombardment of nitrogen atoms. Nitrogen atoms are converted to ¹⁴C isotope (n + ¹⁴₇N → ¹⁴₆C+ p) when they react with thermal neutrons (n) produced by cosmic rays interfacing with Earth’s atmosphere. Radiocarbon enters the nutrient cycle of Earth through photosynthesis afterwards [138,143].

Correction based on ¹⁴C can identify NSR outliers and resolve the upper boundary of CSR more accurately, compared to background correction [35]. Sitewide average NSZD rate estimated using location-specific ¹⁴C correction was 2.5 times less than that achieved using background correction (1.2 vs. 3.0 μmol CO₂ m⁻² s⁻¹, respectively) [114]. A radiotracer study using [¹⁴C]-bicarbonate and [2–¹⁴C]-acetate could provide an understanding about methanogenic biodegradation of crude oil under high temperature and pressure, demonstrating that syntrophic acetate oxidation coupled with hydrogenotrophic methanogenesis was the dominant metabolic pathway under this circumstance [97]. Nonetheless, ¹⁴C analysis is a costly and labor-intensive method and requires special laboratory equipment [35]. CO₂ should be physically recovered from soda lime using an acid extraction method, which requires an additional laboratory step [47,56]. Furthermore, this method tends to overpredict CSR because the gas transport distance between contaminant source and emission at grade is not accounted for in the underlying calculations; it is assumed that the percent of modern-C (pMC) does not change as gas travels in the subsurface. To minimize the associated error, it is recommended that only the pMC values which were measured most closely to the surface are incorporated into data analysis [35]. Isotopic fractionation can influence ¹⁴C as well. However, unlike δ¹³C, this effect is so small (within the analytical precision of the measurement, e.g., ±0.01 times pMC) that it can be safely overlooked [135].

4.4.3. Conversion of CO₂ Efflux to NSZD Rate

NSZD rate is estimated stoichiometrically using the quantified CO₂ efflux. ITRC recommends using octane (C₈H₁₈) as the representative PHC (Equation (12)) [22]:

$$2 {\mathrm{C}}_8{\mathrm{H}}_{18}+25 {\mathrm{O}}_2 \to 16 \mathrm{C}{\mathrm{O}}_2+18 {\mathrm{H}}_2\mathrm{O}$$

(12)

$$2 {\mathrm{C}}_{10}{\mathrm{H}}_{22}+31 {\mathrm{O}}_2 \to 20 {\mathrm{CO}}_2+22 {\mathrm{H}}_2\mathrm{O}$$

(13)

However, ITRC clarifies that a different PHC can (and in fact should) be selected to suit the site-specific hydrocarbon compositional profile [22]. Decane (C₁₀H₂₂) has been widely used as the representative PHC in NSZD studies [11,13,29,34,37,46,47]. Ten moles of CO₂ are produced per one mole of decane degraded (see Equation (8)) [29,37]. Given decane density of 0.73 g/mL, 1 μmol CO₂/m²/s is equivalent to 6138 L LNAPL/Ha/yr [37]. NSZD rate is commonly reported as L LNAPL degraded/Ha/year or U.S. gallon LNAPL degraded/acre/yr [25].

5. Thermal Monitoring: A Recent Technique

Thermal monitoring is an emerging technique that uses subsurface biogenic heat as a surrogate for NSZD rate measurement [12,22,26,46,78,81,144,145,146]. Temperature was first used in 2014 as a tool to evaluate aerobic degradation of PHC within a subsurface [146]. However, due to the similarity between enthalpies (ΔH_r^o) of aerobic PHC degradation and methane oxidation (see Table 1), it is likely that the temperature anomalies captured by Sweeney and Ririe [146] attributed to aerobic PHC decomposition were in fact due to NSZD by methanogenesis that progressed through methane oxidation consecutively (see Section 2.1).

The primary source of heat in both saturated and unsaturated zones of an LNAPL contaminated subsurface is methane oxidation [12]. It can be understood from NSZD thermodynamics discussed in Section 2.4.2. The enthalpy of methane oxidation is 32 to 750 times higher (with decane and hexadecane as representative LNAPL, respectively [12,22,78]) than the enthalpy of methanogenesis, rendering the heat generated from the latter negligible (0.2% of the overall heat generated from LNAPL degradation [12]). It is therefore thought that thermal groundwater plume which overlaps the center of the oil body is most likely transported from the overhead unsaturated zone where methane oxidation occurs [12]. Using Gaussian process regression (GPR), Warren and Bekins [12] showed that a thermal plume forms near the center of oil body and expands along the groundwater flow path, without the temperature retrieving background values (6.3–6.5 °C) even up to 250 m downgradient.

This technique was advanced by a number of researchers [46,78,81,144] and patented in 2015 [145]. Other researchers noticed subsurface temperature anomalies in LNAPL contaminated areas as well [12,26,146,147,148]. Biogeochemical processes cause marked effects in a contaminated subsurface including alteration in temperature, among many other signatures (e.g., changes in petrophysical properties, mineralogy, solute concentration of pore fluids, etc.) [10]. At LNAPL spill sites, the vadose zone temperature establishes a correlation with contaminant respiration [26], which was about 2.7 °C greater than background temperature [12]. A thermal plume forms at the contamination source zone in the saturated zone, with the highest temperature overlapping the center of oil body A temperature difference of 1.2–2.9 °C between the thermal plume and background location has been recorded [12]. The temperature of the smear zone is consistently higher (2–2.5 °C) than that out of the LNAPL footprint [146].

A variety of devices are used for monitoring subsurface temperature, including thermocouples [12,26,46,78,81], thermistors [12,146], resistance temperature detectors (RTDs) [149] and bandgap-based digital temperature detectors [150], with thermocouples being used most frequently in NSZD studies. However, a standard approach or protocol to select for a preferred device does not exist [22]. Analog thermistors are suitable for preliminary assessments to have a snapshot of site conditions. The use of thermocouples, RTDs and band-gap sensors allows simple remote data acquisition, less human error and yearlong measurements [150]. Data collected using thermistors and bandgap sensors agreed well, with a difference of 0.02–0.03 °C which was within the sensors’ resolution limits [150]. Comparing RTDs and thermocouples, the former is preferred especially in cold climates, because: (1) they result linear data over wide temperature ranges, hence linearizing will be less challenging; (2) they do not need cold junction compensation due to their electric characteristic features [151].

Data reduction and analysis steps include: (1) background correction—correcting the recorded temperatures by removing background thermal effects; (2) applying energy balance to the background-corrected temperature profile of a reference volume in the monitored zone—this reference volume usually has the ground surface as the upper boundary and the lower boundary is chosen below the LNAPL source zone; (3) combining the results of energy balance and LNAPL degradation stoichiometry using Equation (8), to resolve the NSZD rates [46,145]. Key assumptions of the energy balance are outlined by Karimi Askarani et al. [46]. Similar to CO₂ trap method (Section 4.4.3), decane is commonly selected as generic hydrocarbon in this technique [22,46,78,81]. Although, other sources are used as the representative PHC as well, such as hexadecane in Warren and Bekins [12]. Finding the most appropriate PHC representative for this substitution is necessary to establish an appropriate understanding of site-specific LNAPL composition and realistic NSZD rates. Advantages and limitations of this method are discussed in Section 3.

6. Conclusions and Outlook

NSZD is an emerging sustainable and cost-effective in situ technique for monitoring and regulatory compliance of LNAPL bioremediation in decommissioning and reclamation industry. Depending on restoration objectives, NSZD can be used as either the primary or sole LNAPL management technique. To achieve this goal, the NSZD rate should be quantified accurately and precisely.

Based on the characteristic features of NSZD sites, four principal methods are developed so far to resolve NSZD rates, including CG, DCC, CO₂ trap and thermal monitoring. Each of these methods has a number of advantages and limitations, which are either common between them or unique to them. For instance, both CG and DCC can provide instantaneous site readings. However, data reduction and analysis are time consuming and labor intensive in both methods. While CG can only account for diffusion, DCC can capture gas effluxes caused by diffusive and advective gas transport mechanisms. Thus, CG data are not representative of site conditions where advection is the dominant gas transport mechanism. Refer to Table 2 to learn more about advantages and limitations of each method.

Part of the limitations are caused by the equipment or the method protocol. However, a significant portion of disadvantages lies in environmental conditions, such as wind, barometric pumping, precipitation, solar radiation, soil texture and moisture, water table fluctuations, etc. The subsurface constantly interacts with these environmental factors, resulting in substantial changes in the subsurface characteristics even at high-resolution timescales. To that end, these factors introduce high degrees of temporospatial heterogeneity into the site conditions. This introduces considerable variations in NSZD rate measurements temporally and spatially. This can be interpreted as an inherent characteristic of NSZD phenomenon rather than shortcoming of the quantification techniques per se.

Investigations performed on relative precision and accuracy of these four techniques indicated that CO₂ trap data agreed well with those obtained by DCC and CG (linear regression analysis: R² > 90%). It fluctuated considerably for thermal monitoring (R² = 10–50%), depending on the time of the year these methods were used. It is because CO₂ trap data is more sensitive to seasonal circumstances, especially those affecting ground surface characteristics. Trap seems to overestimate NSZD rate compared to all other techniques. However, this should be validated with more data points, especially at higher CO₂ efflux ranges. In the future, precision studies on all these techniques should be performed systematically. All four methods should be investigated simultaneously at the same location to improve comparability of the results. All these techniques can introduce bias into NSZD rate measurements, as all of them are affected by certain factors (outlined in Section 2.3 and discussed in Appendix B) one way or another. It is proposed that continuous monitoring of NSZD rate (e.g., at least a year) can average out positive and negative biases and converge the apparent and actual NSZD rates [46]. This is the point where thermal monitoring becomes an interesting technique, as this method is deemed to introduce fewer biases due to climatic factors into NSZD rate measurements, reduces field work and is cost-effective for yearlong studies.

Certain outstanding questions and concerns include the following, which can provide directions for future research:

• The role of predatory organisms on biotic processes at NSZD sites has not been studied and comprehended so far. Predation has been categorized as a less critical factor on biotic processes at NSZD sites [39]. This assumption can be due to difficulties in culturing strictly anaerobic protozoa, which may have caused lack of experimental proof for the speculations about protist roles in anaerobic ecosystems.
• Nutrient amendment as a site management strategy should be further investigated. These questions should be answered first: does nutrient supplementation increase NSZD rates in the subsurface? What are the most effective nutrient compositions to improve natural biogeochemical processes at NSZD sites? Does optimal nutrient composition vary from site to site and does it depend on LNAPL composition? Secondly, it should be deduced if nutrient addition will have any adverse effect on NSZD rate, either through inhibiting the microbial processes or interfering with measurements.
• Methanotrophs are seen to thrive in micro-aerobic environments (0.5–2% O₂) rather than highly oxygenated environments. Oxygen inhibition at higher O₂ concentration is suspected and is recommended to be quantified [152]. This should be validated through further investigation.
• Anaerobic methanotrophy occurs at considerable rates in anoxic coastal and wetland sediments, presumably via syntrophic relationship between methanotrophic archaea and sulfate-reducing bacteria with elemental sulfur transfer [153,154,155]. The question is whether anaerobic methane oxidation (AOM) can be happen at NSZD sites. Additionally, it should be further investigated whether it can be responsible for significant methane oxidation, given that the bottom layers of the unsaturated zone, confined between the aerobic vadose zone and the water table, constantly experience anaerobic condition.
• Systematic, comprehensive studies should be planned to investigate not only the main effect of influencing factors, but also the interactions between them. The existence and the degree of these interactions should be validated to improve the accuracy of NSZD rate quantification protocols (i.e., actual NSZD rate ≌ apparent NSZD rate).
• These systematic studies (explained in the previous numbered point) should attempt to study the influencing factors methodologically and analyze the main and interaction effects quantitively and qualitatively. The current literature requires more data to validate and make solid conclusions about the effect of each influencing factor.
• As described in the Introduction section, critical review of the underlying factors and providing thorough explanation of the underlying phenomena was out of the scope of this literature review article. There is a need in the current literature for critical review articles regarding different aspects of NSZD, which can explain the underlying phenomenon relevant to each factor.
• How significant is the contribution of “signal shredding” to high temporospatial variability of NSZD rates? This question merits further research given that no study has been dedicated to this phenomenon at NSZD sites to date.
• It is hypothesized that surface properties can influence actual NSZD rate. Tracy [47] argues that based on Le Chatelier’s principle, accumulation of reaction by-products (i.e., gas vapors) due to surface capping can establish an equilibrium at lower rates, resulting in decreased rates of all forward reactions (i.e., biodegradation, dissolution and volatilization) [47]. Although this hypothesis is theoretically valid, it requires further investigation to determine whether and how subsurface heterogeneity and complexity can determine the importance of this factor in regulating actual NSZD rates.
• What is the magnitude of the effect of precipitation on CO₂ trap readings? It is anticipated that CO₂ trap deployment during precipitation can form preferential gas transport pathways, resulting in NSZD rate overestimation. Deployment immediately after precipitation can cause underestimation, as soil saturation can prevent gas emission at grade.
• What is the relative accuracy of CO₂ trap models, given that the two available models (i.e., the Keith and Wong model and the McCoy model) have some differences in scale (the former has a 10-time larger surface area than the latter) and configuration (closed top in the Keith and Wong model)?
• Based on what criteria, protocol or standard should CO₂ traps be distributed at site to resolve the most accurate NSZD rate and also result in consistency in the use of trap and an improvement in comparability of the results of different studies? Do the density of CO₂ traps (i.e., number of traps per unit LNAPL impacted area) and the grid pattern (relative distance and placement of the traps) alter the results (i.e., the obtained NSZD rates)? Does the required number of traps depend on distribution (i.e., the magnitude and shape of LNAPL area) and composition (e.g., light vs. heavy LNAPL) of the LNAPL body? An urgent need exists to develop a grid design protocol and standard for this technique if the method is to be used frequently and is to be expected to provide accurate representation of the field conditions.
• Accuracy of the downstream analyses of the retrieved CO₂ trap from the field (e.g., dry weight change analysis and radiocarbon isotope analysis of the soda lime) merits further research. This is to establish a correct understanding of the error introduced into data during laboratory procedures.
• In thermal monitoring technique, a variety of devices can be used for monitoring the subsurface temperature, including thermocouples, thermistors, RTDs and bandgap-based digital temperature detectors. The thermocouples are used most frequently. A standard approach or protocol to select for a preferred device does not exist [22].
• Further research should be dedicated to eliminating or reducing the dependency of all methods on background correction. Background correction is the simplest and least expensive approach that resolves the origin of CO₂ efflux or biogenic heat used as NSZD rate surrogates. Nonetheless, finding the most reliable background location can be challenging, due to a number of reasons: (1) uncertainty about the mobility and areal extent of LNAPL in many active and decommissioned industrial sites; (2) similarity between CO₂ effluxes of background location and LNAPL impacted area. The latter can lead to ambiguity about the occurrence of NSZD at sites with relatively low LNAPL loss rates. It can result in negative background-corrected values making the impression that the LNAPL compounds are being produced.
• The NSZD is a bioremediation and reclamation technique at an embryonic stage. Therefore, enough information was not available in the current literature about capital cost, operating cost and revenue based on the technical and financial input parameters to feed into a techno-economic analysis of NSZD. It can be fit into the scope of the future literature review analysis of NSZD. In addition, it is recommended that the future critical literature review articles include a bibliographic analysis of the topic, to determine the new knowledge and the direction of the future research closer to the most recent updates that are being released to the public every day.