Natural Source Zone Depletion (NSZD): from process understanding to effective implementation at LNAPL-impacted sites

: This paper describes the proceedings of a Special Session on Natural Source Zone Depletion (NSZD) at the AquaConSoil conference, held virtually in June 2021. It spans research across Europe, the United States of America and Australia.NSZDhas been described as the ‘ combination of processes that reduce the mass of LNAPL light non-aqueous phase liquid) in the subsurface ’ . LNAPL NSZD research and investigations have been focused on a range of hydrocarbon products, such as gasoline, diesel, jet fuel, as well as crude oil. Key NSZD processes include aerobic biodegradation, fermentation and methanogenesis of LNAPL constituents, dissolution of LNAPL constituents into groundwater and volatilization of LNAPL constituents into the unsaturated zone. In turn the generated methane, carbon dioxide and dissolved/volatilized constituents can be cycled and biodegraded in the unsaturated and saturated zones. Importantly these physical, chemical, and biological processes can act without human intervention to reduce the NAPL mass and toxicity. Over time NSZD can both reduce LNAPL mass, and change its chemical composition resulting in risk reduction, and ultimately source depletion. Methanogenesis of organicmaterials haslong been recognized in municipallandfills andnatural anoxicenvironments, such as peat and wetlands. Recognition of similar processes in LNAPL source zones in the past decade along with high rates of aerobic biodegradation observed in unsaturated zones above LNAPL-impacted areas, has significantly revised the conceptual model of LNAPL source zone behaviour and persistence. Several NSZD monitoring approaches have been developed and are being applied in field studies. While the quantitative NSZD rates derived can vary between techniques, they all demonstrate that NSZD LNAPL removal can exceed that delivered by engineered LNAPL recovery techniques, particularly for mature LNAPL bodies.

The ubiquitous use of high volumes of crude oil and liquid hydrocarbon fuels represents a hazard to water quality, including to groundwater resources. Their potential accidental release to the subsurface during transport, storage or use can result in the light non-aqueous phase liquid (LNAPL) product/crude oil reaching the watertable, partially penetrating and laterally spreading there due to its immiscible, less-dense-than-water nature. Conceptual models of LNAPL fate in the subsurface are sensitive to the hydrogeological setting and site-specific case detail (CL:AIRE 2014; Tomlinson et al. 2017). Even so, there are broad commonalities. A radially extensive subsurface 'source zone' of LNAPL invariably results, which contains significant mass and a diversity of hydrocarbon constituents. This source zone comprises both mobile LNAPL, dynamically moving up-and-down with a fluctuating watertable and perhaps continuing to spread laterally, but also immobile LNAPL permanently or temporarily entrapped by capillary forces in the subsurface pores or fractures present (Lenhard et al. 2017). Mobile LNAPL may migrate and pose risks to nearby boreholes and surface waters, and discharge as LNAPL oily sheens. The LNAPL source zone mass depletes with time but may generate additional risks that require management. Volatile constituents can pose vapour intrusion risks to nearby buildings and their occupants (Patterson and Davis 2009;Rivett et al. 2010). The LNAPL slowly dissolves into water forming dissolved-phase plumes in groundwater, especially of the more water-soluble aromatic hydrocarbon components such as BTEX (benzene, toluene, ethylbenzene, xylenes) that could pose risks to drinking water supply wells (Thornton et al. 2013). Such risks are often considerably reduced by virtue of many petroleum hydrocarbons being relatively biodegradable in the subsurface leading to the 'natural attenuation' of dissolved plumes that may form the basis of 'monitored natural attenuation' (MNA) remediation strategies (Wiedemeier et al. 1999;Rivett and Thornton 2008). This is often in conjunction with source zone intervention technologies that aim to reduce LNAPL source zone mass and longevity of site management (e.g. CL:AIRE 2014).
A widely accepted paradigm has been that the source zone LNAPL mass is large and, as many hydrocarbon constituents are not highly soluble, natural depletion of the LNAPL body via dissolution into a soil water or groundwater phase is slow. Increasingly though, site observations are challenging this paradigm. The pioneering work of  and Lundegard and Johnson (2006) corroborated by site observations internationally (e.g. Rayner et al. 2020) endorse that the natural depletion of LNAPL source zones may be much higher than hitherto assumed. The primary contentions being that the significance of methanogenesis-based biodegradation of hydrocarbon LNAPL mass and associated aerobic biodegradation of vapour-based losses have been underestimated. The incumbent paradigm shift in understanding now offers: the potential for so-called 'Natural Source Zone Depletion' to be deployed as a source zone management approach that uses rate measurements and monitoring to verify source depletion (Garg et al. 2017; CL:AIRE 2019); reduced operational timeframes or needs for direct intervention engineered LNAPL recovery techniques; a baseline to quantify the relative efficacy of those techniques against; a typically more sustainable approach (CL:AIRE 2010; Bardos et al. 2011;Dottridge and Smith 2021) especially for mature LNAPL bodies; and better constraints estimation of timeframes required for LNAPL site management thereby improving long-term decision making (Sookhak Lari et al. 2019a).
This paper describes the proceedings of a Special Session on Natural Source Zone Depletion (NSZD) at the AquaConSoil 2021 conference, held virtually in June 2021. Presentations made by the practitioner, problem holder and researcher authorship herein drawn from their North American, European, Australasian, and wider global experiences at LNAPL-impacted sites aimed to provide an upto-date overview of NSZD under the Session title 'NSZD: Developing from Principles to Effective Implementation'. Under the umbrella of problem conceptualization and risk-based LNAPL management, the principles of NSZD were set out highlighting key processes and controlling factors followed by a review of the five main methods that have gained, or are beginning to gain, traction for monitoring and assessing NSZD. The most recent of these methods LNAPL Compositional Analysis being described in detail. The Session concluded with a set of implementing Australian case studies that allowed assessment of the efficacy of NSZD versus active engineered recovery of LNAPL. A similar format is adopted herein with the above aims inherent. Our overall ambition being to communicate the nature, current main approaches, and the potential significance of NSZD to the management of LNAPL impacted sites.

Principles of NSZD: Processes, controlling factors and composition change
In recent years NSZD has emerged as an important remediation approach for LNAPL sites. NSZD has its roots in Monitored Natural Attenuation (MNA) but with a different focus. MNA, originally developed in the 1990s, focuses on 'how far will a plume migrate' by identifying key attenuation fate and transport processes in the dissolved phase, determining plume stability (whether expanding, stable, shrinking, or exhausted), and evaluating indirect indicators of attenuation such as groundwater geochemical conditions (Prommer et al. 1999(Prommer et al. , 2002. NSZD, a much newer remediation approach, focuses on 'low long will a plume/source zone persist' and is designed to measure the rate that LNAPL itself is depleted, either directly by LNAPL chemistry changes, or indirectly through soil gas, temperature or other suitable measurements. The LNAPL conceptual model has evolved dramatically since 2010, when a realization occurred that NSZD rates were one to two orders of magnitude higher than estimated only from aqueous dissolution rates. The emergent NSZD conceptual model identified why this was happening: biodegradation reactions via methanogenesis are present in most organic-rich LNAPL source zones. Key analogues which support the concept of methanogenic-dominated systems include the biogeochemistry of wetlands, landfills (Williams 1998), and peat lands and engineered anaerobic digestors at wastewater treatment plants. At LNAPL sites, these important anaerobic and methanogenic NSZD processes are expressed quite differently from conventional MNA biodegradation processes in plumes, and therefore require different approaches to verify the presence of NSZD processes and to measure NSZD rates. Additionally, it was realized that aerobic biodegradation was concomitantly large within the unsaturated zone. This was found for volatile compounds (Franzmann et al. 1999;Davis et al. 2005), for methane degassing from the LNAPL zone  and for direct LNAPL degradation (Davis et al. 1998(Davis et al. , 2013. The overall NSZD conceptual model is built around methanogenic biodegradation within the LNAPL and aerobic biodegradation in the unsaturated zone above the LNAPL but includes a wide range of ancillary biological and physical processes. Figure 1 shows: (i) microbial production of methane and carbon dioxide directly in the LNAPL source zone; (ii) movement of methane and other gases from the saturated zone into the unsaturated zone via diffusion and ebullition of these gases as gas bubbles; (iii) vertical transport of volatile organic compounds (VOCs), methane and carbon dioxide from both the saturated and unsaturated zones vertically upward toward the surface via gas-phase diffusion and in some cases advection; and (iv) biological conversion of VOCs and methane to carbon dioxide by methane and VOC oxidizing microbes at an aerobic/anaerobic interface in the unsaturated zone (e.g. Davis et al. 2021).
Earlier conceptual models of LNAPL attenuation relied on initial dissolution of the hydrocarbon into water before biodegradation of the molecule within a microbial cell could occur. By contrast the production of methane and carbon dioxide by methanogenic microbes in the LNAPL source zone, and degassing of CH 4 , CO 2 and VOCs from the LNAPL zone into the overlying unsaturated zone where aerobic biodegradation can occur, results is significantly greater mass transfer and subsequent biodegradation than dissolved phase processes alone can facilitate. This revised NSZD conceptual model has several important implications for how to measure NSZD (described in the next section) and how NSZD fits into current LNAPL regulatory frameworks.
As described by the U.S. Interstate Technology and Regulatory Council's (ITRC) LNAPL Site Management Guidance (ITRC 2018), LNAPL at impacted sites is typically managed using two broad metrics: LNAPL composition goals (e.g. concentration-based goals for individual LNAPL constituents) and LNAPL saturation goals (e.g. the actual amount of the bulk LNAPL present in the pore space in the subsurface). The former for instance may relate to the reduction of benzene concentrations posing potential risks to people exposed via impacted drinking water consumption or vapour inhalation. The latter address LNAPL migration risks posed by high saturations. While MNA has traditional focused on managing several key site constituents in dissolved plumes, NSZD measurement techniques indicate the rate that bulk LNAPL is being destroyed by NSZD processes, typically in a mass flux type formulation with units of the volume of LNAPL destroyed per area per time (e.g. litres of LNAPL per hectare per year). Therefore, NSZD projects can help understand the longevity of the bulk petroleum impacts such as 'how long will impacts at the site remain?'. By integrating NSZD rates at several locations across a site (in units of volume of LNAPL destroyed per area per year) and the corresponding LNAPL volume per area estimates (in units of volume of LNAPL per area), simple calculations can be made to provide planning-level estimates of the LNAPL longevity. While there is some level of uncertainty in these longevity forecasts, they can be useful to emphasize that: (i) the LNAPL is not a permanent presence at impacted sites, and (ii) natural processes will eventually remove the bulk LNAPL.
Estimates of LNAPL longevity are not the only way NSZD rates are utilized. Another important application is to use NSZD to determine the relative efficacy of on-going or proposed active LNAPL recovery projects. For example, as an LNAPL recovery system approaches the point of diminishing returns (Sookhak Lari et al. 2020), the LNAPL removal rate from NSZD may be significantly higher than the removal rate from the active recovery system. In that event, shutting down the active recovery system can provide significant sustainability benefits (e.g. reduced energy use, carbon dioxide emissions, neighbourhood disruption and cost) without extending the LNAPL bulk longevity significantly. For sustainable remediation assessment at LNAPL sites NSZD is often an appropriate 'base case' against which to evaluate the merits of alternative remediation approaches (CL:AIRE 2010, 2014).

Methods for monitoring and assessing NSZD
The study of the natural attenuation processes of dissolved constituents in plumes led to the development of MNA as an accepted groundwater remediation technique. The onus being to present complementary lines of evidence demonstrating plumes were sufficiently attenuated by natural process prior to reaching receptors. However, the key MNA protocols and regulatory guidance neither addressed nor extended MNA concepts to the NAPL source zones themselves. In general, NAPL source zones were thought to only attenuate slowly and primarily deplete via dissolution of soluble constituents caused by groundwater flowing horizontally through a NAPL source zone. Two commonly used MNA models (BIOSCREEN and BIOCHLOR;Newell et al. 1996;Aziz et al. 2000) utilized this dissolution-only source attenuation approach, leading to predictions of relatively long source remediation timeframes at many sites. Other simple hydrogeological risk assessment models (e.g. Environment Agency Remedial Targets worksheet, 'P20'. Carey et al. 2006) assume a quasi-steady state source term. While this is acceptable for managing impacted sites based on estimates of maximum soluble groundwater contaminant plume length, it does not adequately address the duration of the depleting LNAPL source. Some models including ConSim (Environment Agency 2011) and PHT3D (Prommer et al. 2002) simulate depleting sources and more complex influences, but few couple depleting source mass with time-variant compositional changes such as those described by Thornton et al. (2013).
Then in 2006, two key publications Lundegard and Johnson 2006) reported on a remarkable discovery where much higher attenuation (depletion) rates of LNAPL source zone mass were observed by measuring the vertical downward flux of oxygen from the surface and through the unsaturated zone. Overall, the authors concluded 'Depending on the depth of contamination and other factors, the mass loss rate per unit surface area for those processes varies from c. 0.1 to 1.0 kg total petroleum hydrocarbons per m 2 a −1 (annum). In comparison, mass losses from the submerged part of the source zone and involving ground water transport processes (i.e. source dissolution and plume biodegradation) were estimated to be about approximately two orders of magnitude lower' . The serendipitous discovery was triggered using vertical carbon dioxide fluxes near the ground surface to search for shallow hydrocarboncontaining pits, leading to their puzzling ubiquitous observations of carbon dioxide efflux from the subsurface over the entire LNAPL source zone (as opposed to a few discrete pits) and, the very high rate of this carbon dioxide efflux. This led to soil gas studies that showed significant oxygen consumption within the unsaturated zone profile and carbon dioxide efflux at the surface, attributed to NSZD biodegradation processes (e.g. McHugh et al. 2010).
The two authors named this process 'Source Zone Natural Attenuation (SZNA)', and by 2009 the idea that LNAPL biodegradation rates were 10 or 100 times faster than once thought captured the imagination of the remediation community. Three years later ITRC (2009), a U.S. consortium of industry, researchers, and regulators issued a detailed overview of the technique, but with a new name: 'Natural Source Zone Depletion', appropriately emphasizing the pivotal significance of mass depletion to attenuation of LNAPL source zones. Methods for measurement and estimation of NSZD rates have been collated and promoted (ITRC 2009;API 2017;CRC CARE 2018).
By 2014, the term 'Gradient Method' had come into common use (Adamson and Newell 2014) to describe the method where soil gas concentration gradients are multiplied by an effective diffusion coefficient to obtain the flux of oxygen being consumed or the flux of carbon dioxide being generated by biodegradation processes from LNAPL source zone zones, and then using this flux to stoichiometrically establish the NSZD rate in units of mass biodegraded per area per time.
Ground surface carbon dioxide efflux measurements using a 'Dynamic Closed Chamber' (DCC) then appeared in the technical literature around 2011 (Sihota et al. 2011;Sihota and Mayer 2012) and was then extended to include seasonal temporal NSZD measurements (Sihota et al. 2016). The Dynamic Closed Chamber approach (sometimes called the Li-Cor method from a specific the manufacturer's DCC) was originally developed for soil respiration and climate change studies. For the NSZD use case, a DCC is used to measure surface carbon dioxide efflux in units of moles of carbon dioxide flux per area per time [typically µM m −2 s −1 ] over an LNAPL body. Similar efflux measurements from background locations remote from the LNAPL impact are then subtracted to remove the non-LNAPL CO 2 surface emission flux contribution from degrading natural organic material in soil. Converting the resulting CO 2 flux to petroleum hydrocarbons (based on an assumed 1:1 molar reaction stoichiometry of carbon to carbon dioxide) yields a 'total petroleum hydrocarbon mineralization rate' in units of mass per area per time [typically kg m −2 a −1 or L ha −1 a −1 ] (the NSZD rate).
In 2011, McCoy et al. published a second surface efflux method where a 'CO 2 Trap' was used to measure carbon dioxide efflux rates from LNAPL source zones. This new device, comprised of a small canister containing a permeable carbon dioxide scrubbing sorbent, is driven a few centimetres into the soil profile to capture carbon dioxide leaving the subsurface. The device was designed to integrate the CO 2 efflux signal over an extended period, typically two weeks compared to just a few minutes for the DCC. After being deployed, the CO 2 Trap is sent back to a laboratory where the accumulated mass of sorbed carbon dioxide is measured gravimetrically. To remove the background signal from the respiration of 'modern organic carbon', a novel 14 C isotope analysis was incorporated in the method (McCoy et al. 2011). CO 2 originating from petroleum sources is depleted in 14 C. A final innovation was to report results as the volume of biodegradation per area per year such as gallons per acre per year in the United States, units that were particularly meaningful to practitioners and regulators. The CO 2 Trap method is now supported by a vendor who sells the traps and performs the gravimetric and isotope analysis to provide NSZD rates (E-Flux 2021).
About the same time as the 2014 CO 2 trap paper was published, several groups began work on 'Temperature methods' (Sweeney and Ririe 2014;Warren and Bekins 2015;Askarani et al. 2018;Sale et al. 2018;Askarani and Sale 2020). With this method, the temperature gradient is measured and multiplied by the soil thermal conductivity to yield a heat flux from the biodegrading LNAPL body upwards to the surface and downwards into soil underlying the LNAPL. With the heat flux (generated by exothermic microbial biodegradation) and known heat of reaction for the LNAPL biodegradation processes, the NSZD rate can be calculated. One vendor, ThermalNSZD.com, continuously collects temperature data remotely for telemetry to a central location, and then reports NSZD rates on a web-based dashboard.
The recently highlighted method applied for NSZD estimates, 'LNAPL Compositional Analysis', uses the change in marker constituents in the LNAPL itself over several years to estimate NSZD of LNAPL bodies (DeVaull et al. 2020). This method is discussed in detail in the next section of this paper.
In summary, there were five primary methods used to estimate NSZD rates from LNAPL bodies that were addressed in the workshop: (1) Gradient Method; (2) Dynamic Closed Chambers; (3) CO 2 Traps; (4) Temperature (or 'biogenic heat') method; and (5) LNAPL Compositional Analysis. Several studies have compared some of these methods in field situations, such as Kulkarni et al. (2020), andSmith et al. (2021). A series of multiple site studies show NSZD rates at LNAPL sites ranging from 10s to 100s of gallons per acre per year (100s to 1000s of litres per hectare per year) (Garg et al. 2017;Rayner et al. 2020). Table 1 provides a subjective comparison of the advantages and disadvantages of each method.

Field data analysis of total and constituent NAPL depletion rates
The most recently introduced NSZD assessment method (Table 1), 'LNAPL Compositional Analysis', uses the composition change in the LNAPL itself to estimate NSZD of LNAPL bodies. Petroleum in the environment weathers with a differential change in constituent composition over time (e.g. Brown et al. 2017). Quantitative estimates of total petroleum depletion between two samples from the same LNAPL body separated in time (t 1 . t 0 ) may be estimated from an inverse ratio of measured concentrations (H 0 , H 1 ) for an identified conserved marker in the LNAPL (Douglas et al. 1996).
Estimates of depletion for other individual constituents in the LNAPL include measured concentrations (C 0 , C 1 ) for that constituent in the LNAPL.
The method of Douglas et al. (1996) depends on the presence of an identified presumed-conservative marker compound in the oil. Viable persistent markers have been identified in some crude oils  (Prince, et al. 1994) and in diesel (Johnston et al. 2007), but these markers are not always measurable or conserved, may be absent from some petroleum products, and are often absent from nonpetroleum derived LNAPLs. A method for selecting the best available marker or markers in any petroleum LNAPL is described by DeVaull et al. (2020). LNAPL samples were collected for several years over increments of time and analysed for multiple constituent concentrations in the oil. Regression trends were fitted for each constituent. The constituent rates were rank ordered, and the constituent concentration trends in oil increasing at the greatest rates indicated the best available near-conserved markers for the specific LNAPL body under evaluation. The selected marker concentration trends were used in estimating the total LNAPL depletion. Depletion of specific constituents in the LNAPL body were estimated from the total LNAPL depletion and the constituent-specific concentration trends in the LNAPL.
Application of the compositional analysis method in two case examples (DeVaull et al. 2020) showed estimated depletion of half the initial total LNAPL in 13.6 years for a crude oil release and 7.3 years for a gasoline/diesel release. In addition to expected depletion of volatile and soluble constituents of the LNAPL, significant depletion of higher molecular weight LNAPL constituents was also illustrated by the method and is likely attributable to source-zone methanogenesis.
NSZD estimates from composition change are a measure of the primary depletion of the initial petroleum constituents. This depletion will include LNAPL losses from constituent volatilization, constituent water dissolution, and direct constituent degradation. NSZD estimates from CO 2 surface flux are a measure of ultimate degradation of petroleum to mineral products including contributions from subsurface methanogenesis (CO 2 , CH 4 ) of LNAPL and unsaturated zone aerobic biodegradation of volatilized LNAPL constituents and evolved CH 4 to CO 2 and H 2 O.
Implementing NSZD: Field data on effectiveness of NSZD v. Active LNAPL Recovery

Background and Approach
As part of an Australian approach to NSZD application to LNAPL management (CRC CARE 2018, 2020; Rayner et al. 2020) sitespecific measurements and case studies were undertaken across six sites in Australia. The sites had multiple LNAPL product types and age of spill or release, and different geologies, unsaturated zone depths and climatic conditions (Table 2). For these sites NSZD rates were estimated for crude oil, diesel, jet fuel, gasoline and gas condensate. The NSZD rates were compared to the most representative active product recovery rates from each of these sites to provide context to the scale of the NSZD rates.
LNAPL NSZD measurement methods used at the six sites included a selection of (i) fluxes of carbon dioxide exiting at ground surface using E-flux canisters and Li-Cor chambers, (ii) major gas composition (oxygen, carbon dioxide) depth profiles across the unsaturated zone and within open existing wells, (iii) biogenic heat production measured by temperature profiles in existing boreholes (manual thermistor and iButton) and buried thermistors. A general description of these methods can be found in CRC CARE (2018). Some NSZD rates were estimated twice over a year (May and September). LNAPL recovery occurred at sites by a range of methods including skimming, skimming with drawdown, skimming with drawdown and vacuum enhanced recovery, total fluids recovery and multiphase extraction (MPE). Table 3 provides a summary of the ranges of NSZD rates across all methods adopted across each site. Not all methods were able to be deployed at all sites (see Table 4 for which were used at each site). Estimated NSZD rates for Site A ranged from zero or low values for some methods up to 34 000 L ha −1 a −1 based on E-Flux measurements and as high as 36 000 L ha −1 a −1 based on the gradient method for CO 2 flux using vapour multiport monitoring data. The upper estimates of the range of NSZD rates from temperature data (maximum of 4471 L ha −1 a −1 ) were, in general, lower than those determined by other methods. Reasons for the discrepancies between temperature and other methods might be resolvable, for example, by measuring the subsurface thermal conductivity of different strata in the profile.

Summary of Rates across the Sites
For Site B, NSZD rate estimates varied between 1000 and 17 000 L ha −1 a −1 across both crude oil and diesel NAPLs (Table 3). Temperature data were consistent across multiple thermistor strings  . Note that the back calculated NSZD rates in L ha −1 a −1 are averages across the area unlike the raw NSZD rate ranges seen in Table 3. *Site B estimates presented relate to diesel impacted zone (not crude oil). † Veronoi polygons are used to calculate a site-wide spatial average). Mass is in metric tonnes assuming a density of 0.703 kg L −1 . ‡ Attempted historically but deemed not successful. § Attempted during this study but background corrected rates could not be obtained. || May. ¶ September. and measurement periods. Li-Cor carbon dioxide fluxes gave NSZD rate estimates that supported temperature estimates. Gas-gradient based NSZD rate estimates were highest. Crude oil NSZD estimates had a larger range (with the lowest and highest rates) than estimates for diesel.
For Site C, higher estimates of NSZD rates resulted from using temperature and surface CO 2 flux methods. Li-Cor surface CO 2 flux measurements were recorded at only one monitoring site, and the average NSZD rate was 866 L ha −1 a −1 in May and 12 700 L ha −1 a −1 in September of the same year. E-Flux estimates ranged from 1141 to 10 185 L ha −1 a −1 , with an average of 3599 L ha −1 a −1 . Temperature estimates were more than three times higher in May 2019 (maximum of 12 531 L ha −1 a −1 ) compared to September 2019 and the reason for this has not yet been determined. For Site C, compositional changes from 2011 to 2019 suggest estimated mass losses (due to volatilization, dissolution, and biodegradation) of between 3 and 21% over the eight-year period.
For site D, NSZD rates were determined using temperature profiles, three E-Flux deployments of four traps, and from in-well sampling of oxygen. In well oxygen data yielded the lowest NSZD rates (280 to 1110 L ha −1 a −1 ) and temperature profiling gave the highest NSZD rates (2837 to 5140 L ha −1 a −1 ). LNAPL compositional changes (due to volatilization, dissolution, and biodegradation) from 2013 to 2018 suggest estimated relative mass losses of between 1 and 24% over that five-year period.
For Site E, very high estimates of NSZD rates were determined using the Li-Cor and E-Flux carbon dioxide surface flux methods. The rates ranged from 47 339 to 164 634 L ha −1 a −1 and the ranges of Li-Cor and E-Flux results were similar.
For Site F, in well oxygen data yielded the lowest NSZD rates (0 to 1600 L ha −1 a −1 ) and temperature profiling gave the highest NSZD rates (0 to 15 000 L ha −1 a −1 ). There was good agreement between manual, short term (60 min) thermistor measurements and longer term deployed, monthly averaged iButton data. LNAPL compositional changes indicate potentially large mass losses (due to volatilization, dissolution, and biodegradation) of between 28 and 99% (since release). Table 4 has estimates of NSZD rates spatially averaged across the area of LNAPL for each site over which NSZD is estimated to occur, and these are listed along with active LNAPL recovery rate estimates.

Comparison between NSZD Estimates and LNAPL Recovery
For Site A, active recovery of LNAPL on site was a maximum of 3.7 t a −1 [960 L ha −1 a −1 ]. Estimated NSZD rates were between c. 1 and c. 40 times higher than active recovery for temperature methods and vapour port gas compositions respectively. It is noted that active recovery was restricted by operational infrastructure limiting access to the entire extent of the subsurface LNAPL plume. For Site B, at the diesel location, overall NSZD rates were estimated between 120 and 280 t a −1 for an impacted area of 28 ha. When compared with the active diesel recovery rates of about 17.5 t a −1 , NSZD rates were around 7-15 times higher. For Site C, across the 18 ha of LNAPL impacted area, overall NSZD rate estimates were 11-69 t a −1 . The maximum active LNAPL recovery rate from a range of remedial efforts was 3.8 t a −1 . Overall the NSZD rate is 3-20 times larger than this maximum active LNAPL recovery rate.
For Site D, active remedial efforts recovered around 5.8 t a −1 over 6 months of operation; during MPE of gasoline, however, these estimates can be variable due to the difficulty of measuring mass fluxes of volatile compounds at relatively high flow rates. In-well oxygen data yielded the lowest overall NSZD mass depletion of about 1.4 t a −1 with temperature profiling giving the highest overall NSZD mass depletion of 6.9 t a −1 over the assumed impacted zone at the site. Overall, the NSZD rate estimates were c. 5 times lower than LNAPL recovery estimates if based on oxygen data alone, but NSZD rate estimates were similar or higher compared to LNAPL recovery rates if based on temperature profile data. For Site E, active remedial efforts removed c. 0.2 t a −1 of LNAPL. Overall, the NSZD mass loss estimate range for a 175 m 2 area of the site was 1.1-1.2 t a −1 which were approximately seven times the LNAPL recovery rate. If all the surface flux data across the site within the plume boundary (0.48 ha) is considered, the mass flux estimates of NSZD are 14 t a −1 for LI-COR and 20 t a −1 for E-Flux. For Site F, recent remedial efforts had recovered c. 1.7 t a −1 of LNAPL. In well oxygen data yielded the lowest overall NSZD mass depletion rate of 0.40 to 0.68 t a −1 and temperature profiling gave the highest overall NSZD mass depletion rate of 3.4 to 4.3 t a −1 over the assumed impacted zone at the site. Overall, NSZD rates from gas gradient analysis were 2.5 to four times lower than LNAPL recovery rates but estimates from temperature data suggest NSZD rates were 2.0-2.5 times higher than LNAPL recovery rates.

Discussion and Conclusions
This study demonstrates different techniques for measuring rates of NSZD were found to be variable both spatially and temporally, but that ranges were similar to those found elsewhere (e.g. Garg et al. 2017;Davis et al. 2022). The variability may be due to the nonuniform distribution of LNAPL across an LNAPL source zone and seasonal impacts of variations in parameters affecting NSZD processes such as water table elevation and unsaturated zone moisture contents. As all sites demonstrated some level of NSZD, zero or very low rates may be real but may also be the result of the application of an inappropriate method, such as using a surface flux approach over a heavy clay profile, or the selection of nonrepresentative background locations for the temperature profiling method for example. The challenge of finding and applying suitable background sites for NSZD rate correction in all methods except E-Flux was evident across five out of six sites and this can have a profound effect on estimated rates. Additionally, having site specific measurements of relevant parameters needed for NSZD rate calculation (e.g. thermal properties for temperature methods and diffusion coefficients for gases) increases confidence in estimates and allows for detailed comparisons to be made between methods.
Despite spatial, temporal and method variability, Sites A, B, C and E all measured NSZD rates that were 2 to 80 times higher than LNAPL recovery rates. At Site D, all in-well soil gas and E-Flux NSZD rates were below LNAPL recovery rates by a factor of two to three; however, thermistor temperature methods were c. 20% higher than LNAPL recovery rates. The primary unknown uncertainties in determining these rates are the thermal conductivity and gas diffusion properties of the subsurface which were estimated from literature values. For Site F, NSZD rates determined from temperature profiles were 2-3 times higher than LNAPL recovery rates, while NSZD rates estimated from in-well oxygen data were 2.5-4 times lower than LNAPL recovery rates. Davis et al. (2022) present data that shows gasoline NSZD rates were higher than those for diesel for fresher releases, but that diesel NSZD rates were comparable and higher than those for gasoline after 10-20 years of LNAPL weathering. The transition in rates and processes is subject to on-going research.
Overall, the case studies investigated here support the growing body of evidence that in many instances, especially at mature LNAPL sites, LNAPL mass removal through NSZD processes exceeds those achieved through active remediation efforts, and NSZD is likely to represent the most sustainable remediation approach under these circumstances, where potential risks due to exposure are under control, or can be mitigated using institutional controls (ITRC 2016).

Way Forward
To create greater confidence and adoption of NSZD as an LNAPL management option, several aspects warrant future investigation: (1) That NSZD rates will be sustained or 'sufficiently sustained' over possibly decades to warrant site closure under a NSZD management option. Endpoint estimates for active remediation of LNAPL are now available (Sookhak Lari 2018, 2019a, b, 2020. Also, recently (Sookhak Lari et al. 2021), it has been shown to be feasible to establish a computational platform for quantifying the longevity and trends in NSZD rates as LNAPL mass is depleted over time. Such a platform promises to help define and optimize the critical transition from active remediation to NSZD management at sites and the cumulative mass depletion and risk-mitigation benefits of both.
(2) To systematically evaluate different NSZD estimation methods to reduce uncertainties including a focus on the use and value of pre-existing wells and infrastructure at sites v. specially installed infrastructure. Much has been done to advance and compare a variety of methods (e.g. Sihota et al. 2011;Ririe 2014, 2017;Sookhak Lari et al. 2017;Sale et al. 2018;Kulkarni et al. 2020, but systematic side by side comparison of in well v. buried multi-depth gas sampling and temperature measurement infrastructure is needed. Additionally, the contributions to NSZD estimate uncertainty need to be understood and reduced; for example, defining the value of site-specific measures of soil thermal and diffusive properties. (3) Somewhat aligned to 2., trends in rate estimates due to temporal influences, including daily, seasonal and other episodic events such as rainfall need to be understood. Sihota et al. (2016), evaluated seasonal changes via manual sampling, however more rapid changes may also occur. For example, Figure 2 shows, at 1 m below the ground in a 3.5 m unsaturated zone overlying aviation gasoline fuel, that oxygen concentrations determined from buried oxygen probes (Patterson and Davis 2008) vary between 0 and 10% several times over a two-month period. Where manual measurements of oxygen depth profiles might be used to calculate NSZD rates, such fluctuations would lead to twofold changes in rate estimates (for example using data from end of October compared to 4 November, and at other times). (4) To further LNAPL compositional changes as a holistic and direct measure of LNAPL mass losses; a better codification of standard procedures in evaluation of LNAPL composition, including ensuring consistency in laboratory data analysis for samples separated by long elapsed time periods.

Conclusions and recommendations
A summary of presentations and discussions from a Special Session of the AquaConSoil 2021 conference on NSZD are presented. The background on LNAPL risk-management, development of NSZD principles and good practice, overview of NSZD measurement techniques, and cases studies were presented. The recognition of NSZD processes, and in particular LNAPL methanogenesis, methane ebullition, and the high rate of aerobic biodegradation in soil unsaturated zones above LNAPL zones has materially changed the conceptual model of sites with LNAPL in the shallow subsurface. At many sites where the presence of LNAPL does not pose a direct risk to human health or the wider environment, NSZD has the potential to be a sustainable riskmanagement solution.
The results of comparison of different NSZD monitoring approaches have shown significant variability between NSZD methods, but NSZD rates and ranges in variability have also been shown to be similar in the USA, Australia and Europe.
Mass loss estimates from changes in LNAPL composition have proved very insightful as the technique involves direct measurement and interpretation of the LNAPL itself and is therefore a more direct measure than other methods. This has been shown in the Australian case studies  and in a systematic way by DeVaull et al. (2020). In addition, it can distinguish different processes by the relative change in different compounds representing different properties and this can be correlated with the subsurface conditions encountered on site. In general, LNAPL composition changes are only reported as relative percent mass loss and more work is needed to fulfil the potential of this technique. Comparison of percent loss estimates from composition change with the depletion rate estimates of the other methods requires a measure or estimate of initial spill mass or mass in place at some point in time. In this study absolute mass losses were not calculated as the mass of LNAPL in the subsurface relating to the relative mass loss period was not known.
Further research and illustrative case studies are recommended to further support and encourage appropriate use of NSZD to deplete residual LNAPL sources.
Acknowledgements The views expressed are those of the authors and may not reflect the policy or views of their employing organizations. The authors would like to acknowledge the contributions of Poonam Kulkarni of GSI Fig. 2. September to November 2020 time series of oxygen concentrations from oxygen probes buried at a range of depths below ground (0.1, 0.5, 1.0, 1.5, 2.0 and 2.5 m) in a 3.5 m unsaturated zone overlying aviation gasoline LNAPL distributed from 3.0 to 4.0 m below ground. The site location is denoted D9.