Multimethod Analysis of NSZD and Enhanced SZD by Solar‐Powered Bioventing at the Guadalupe Restoration Project

Active remediation at sites with light non‐aqueous phase liquid (LNAPL) often leaves residual hydrocarbons in the subsurface, necessitating long‐term management. While much of the effort in recent years has focused on natural source zone depletion (NSZD) as the primary method for demonstrating continued hydrocarbon removal, the same data collection methods can quantify biodegradation enhancements that can sustainably increase the rate of SZD. This approach has been applied at the Guadalupe Restoration Project, one of the first sites at which NSZD measurements and monitoring technology were demonstrated. Sitewide NSZD quantification was conducted using CO2 efflux measurements and subsurface temperature profiling. The results fell within the range of previously reported estimates that were based on soil‐gas profiling in the early 2000s, demonstrating the viability of the new methods at this site. The data collection methods were then deployed during pilot testing of solar‐powered bioventing. The system used seven 400‐W solar panels to power a regenerative blower that delivered approximately 0.85 cubic meter per minute (30 cfm) air to the LNAPL‐impacted vadose soil near the interface with the groundwater table. Soil‐gas data indicated an upward fanning of injected air toward ground surface. Elevated temperature due to hydrocarbon oxidation yielded an approximate 10.2 kg day−1 source depletion rate above the baseline NSZD mass removal rate over an approximate 30 m (100 ft) radius of influence, which aligned well with a 8.2 kg day−1 rate estimated from CO2 efflux measurements. Introduction of O2 via bioventing substantially increased the LNAPL biodegradation rate from baseline NSZD processes by almost an order of magnitude. The results demonstrate that site management can proceed along a sequenced program that began with aggressive hydraulic recovery of hydrocarbon product, transitions to enhanced SZD in areas with poorly recoverable LNAPL, and then to NSZD without intervention to address residual LNAPL across the full footprint of the LNAPL bodies.


Introduction
Sites contaminated with petroleum in the form of light non-aqueous phase liquid (LNAPL) near the water table are a continuing challenge for cleanup due to the residual contaminant masses that remain after hydraulic recovery reaches a practicable endpoint (Huntley and Beckett 2002;CRC CARE 2015;ITRC 2018).In many cases, natural source zone depletion (NSZD) can be an ultimate management technique, but remediation timelines may not be acceptable to all stakeholders.This motivates possible tran-sition technologies between aggressive hydraulic LNAPL recovery and NSZD.
Much of the recent research for such late-stage LNAPL sites (e.g., those where lateral migration is at equilibrium [ITRC 2018]) has quantified NSZD rates for the remaining hydrocarbon mass using subsurface methods including temperature monitoring (Warren and Bekins 2015;Karimi Askarani et al. 2018), nested vapor well sampling (Johnson et al. 2006), and groundwater geochemical sampling (Johnson et al. 2006;Mackay et al. 2018b), as well as surface techniques such as CO 2 efflux measurements (Sihota et al. 2011).Each of these methods rely on the observation that much of the NSZD occurs via microbial biodegradation, such that heat fluxes, or biodegradation reactant/product fluxes, can be converted to hydrocarbon removal rates.Multiple studies have indicated that NSZD rates are generally in the 0.3 to 5.4 hydrocarbons per square meter per day (g HC m −2 day −1 ) range for late-stage LNAPL sites (Sihota et al. 2016(Sihota et al. , 2018;;Eichert et al. 2017;Mackay et al. 2018b), with higher values observed at sites with transmissive soils and where LNAPL mass in place is relatively high (Davis et al. 2022).Because many LNAPL sites have footprints that are multiple hectares, this rate can equate to large quantities of hydrocarbon removal over time.
NSZD studies have provided evidence that natural hydrocarbon biodegradation is a significant natural attenuation process and can be an important component of LNAPL site management.However, stakeholders may find the pace of NSZD insufficient in some cases.In such situations, it may be valuable to quantify the degree of additional hydrocarbon removal that can be achieved by enhancement of natural attenuation processes.One such opportunity is by introducing additional O 2 to hydrocarbon-impacted soil via bioventing.The approach uses either positive-pressure injection of air (Abu-El-Sha'r and Al-Zou'by 2005) or lowflow soil vapor extraction (Rathfelder et al. 1995) to deliver O 2 to LNAPL bearing vadose zone soils to increase aerobic biodegradation rates.While soil water availability may limit bioventing applicability in some settings (Hinchee et al. 1995), LNAPL is typically present in unsaturated soil near the water table where moisture is not limiting airflow.Thus, so long as the soil permeability is greater than 0.1 Darcy (Leeson and Hinchee 1997) and the contaminant is in contact with the vadose zone (three-phase saturation conditions), the technology is potentially viable for enhancing hydrocarbon removal.LNAPL impacts in soil may include some fraction of mass in the capillary fringe and vadose soil above the water table due to natural water table fluctuation.
Typical characterization of bioventing performance calculates an in situ respiration rate based on O 2 utilization at test wells (EPA 1995) or by estimating a soil hydrocarbon concentration-based attenuation rate (Eyvazi and Zytner 2009).The former requires that both the lateral and vertical zone of treatment be estimated or assumed in order to scale up to a bulk hydrocarbon removal rate.The threedimensional treatment zone is difficult to estimate for most settings, so a scale-up calculation has significant uncertainty.Estimating a hydrocarbon concentration-based attenuation rate requires sequential resampling of soil to scale up to a bulk hydrocarbon removal rate, which is typically not practicable except perhaps as a low-frequency component of a monitoring program.In contrast, CO 2 efflux measurements at ground surface can potentially provide a direct indicator of in-situ hydrocarbon destruction, and subsurface temperature measurements report directly comparable information.Large numbers of samples in space (multiple CO 2 efflux measurements over a large footprint) or time (subsurface temperature measurements repeated at a given location) can increase confidence in rates generalized to larger scales.In addition, such NSZD methods, when combined with soilgas sampling, are likely to provide a more comprehensive picture of enhanced hydrocarbon removal with limited uncertainty.
We investigated how enhancements to hydrocarbon biodegradation can be quantified by traditional NSZD methods, and how the measured rates compare to natural conditions.This study applies NSZD methods to evaluate the influence of introducing O 2 into the subsurface soil via a solar-powered bioventing system at the Guadalupe Restoration Project, where large-scale natural hydrocarbon removal rates from LNAPL were first quantified (Lundegard and Johnson 2006).The study objectives include: (1) characterizing NSZD broadly at the site with updated methods and comparing current rates to those historically measured, (2) quantifying "green remediation" hydrocarbon removal at a localized solar-powered bioventing test area, and (3) evaluating potential strategies for deployment of bioventing and NSZD across the site to optimize site-wide mass removal efforts.

Site Description
The Guadalupe Restoration Project is located at the former 3000-acre Guadalupe Oil Field in central California (Figure 1).During the operating life of the oil field from 1947 to 1990, a refined hydrocarbon (approximate C 10 to C 25 ) LNAPL was used in the field as a diluent to reduce the viscosity of the heavy crude to improve crude transport (Lundegard and Johnson 2006).Leaks from diluent piping and tankage went undetected and resulted in LNAPL bodies forming near the water table underlying the field.While many localized LNAPL bodies have been identified at the field, there are three main bodies (DT, TB9, CP) that comprise the largest footprint, at nearly 50 ha.These three bodies generate dissolved phase plumes measured as total petroleum hydrocarbons (TPH) in the C 10 to C 25 range.The DT, TB9, and CP LNAPL bodies have been the focus of automated hydraulic LNAPL recovery systems since 1998.The systems have alternately run as "enhanced skimming" operations with pneumatic pump intakes set to the LNAPL layer depth in wells, and as total fluids recovery systems with pump intakes set lower into the water column.In recent years, LNAPL recoverability assessments have identified many automated pumping locations that have reached a practical recovery endpoint, with some remaining pockets of higher recoverability LNAPL.Those remaining areas with recoverable LNAPL are currently pumping in a dualphase liquid extraction regime with intentional water table drawdown to increase LNAPL recovery rates in extraction wells (Trihydro 2021).
NSZD rate estimates for this site were first reported in the mid-2000s (Lundegard and Johnson 2006) in a seminal evaluation that became the basis for subsequent guidance on integrating site data to quantify hydrocarbon removal rates by natural processes (ITRC 2009).The initial study employed 10 nested vapor wells to characterize vertical profiles of O 2 , CH 4 , and CO 2 .These profiles, combined with estimated vapor diffusion coefficients, were used to calculate O 2 transport downward through the vadose zone to support biodegradation of hydrocarbons and CH 4 conversion to CO 2 .The study also estimated saturated zone NSZD rates by quantifying fluxes of biodegradation reactants (dissolved oxygen, nitrate, and sulfate) and products (dissolved iron, manganese, and methane) into and out of the LNAPL bodies.The NSZD rates calculated for soil-gas transport in the main LNAPL bodies at the site ranged from 0.4 to 3.0 g HC m −2 day −1 , and the NSZD rates for groundwater transport were approximately two orders of magnitude lower.However, some fraction of LNAPL biodegradation occurring below the water table can result in, for example, CH 4 gas ebullition from the saturated zone into the vadose zone and is thereby accounted for via vadose zone NSZD measurements (Ng et al. 2015).Thus, the attenuation rates observed from groundwater monitoring and vadose zone monitoring should be considered a representation of the total hydrocarbon removal rate from saturated and unsaturated soil zones.

Methods
Data collection was conducted in 2019 to characterize NSZD at the three main LNAPL bodies by CO 2 efflux measurements and subsurface temperature profiles to assess the viability of both, followed by quantification of hydrocarbon removal rates by solar-powered bioventing in a portion of the LNAPL body in the DT area.Data collection methods followed the framework outlined in ASTM (2022).

Sitewide NSZD Characterization
For the initial sitewide effort, CO 2 effluxes were measured at 10 grids across the site (Figure 1): two grids above background areas, six grids above the three LNAPL bodies (CP, TB9, DT), and two grids above dissolved phase plumes.Each grid was comprised of five rows and five columns of CO 2 efflux measurement points, for a total of 25 points per grid.The spacing between grid points was approximately 30 m, for a grid area of approximately 1.5 ha.Soil collars were installed by hand in the dune sands at least 24 h prior to measurements.CO 2 efflux measurements were taken when shallow soils were dry in April and October 2019 following methods detailed in Eichert et al. (2017) with "background correction" (i.e., no radiocarbon sampling).The CO 2 efflux was measured by dynamic closed chamber (LI8100A; LI-COR Biosciences, Lincoln, Nebraska) at each location once in the morning and once in the afternoon for replicates at different times of day.In addition, because two instruments were used for the field effort to decrease overall deployment time, the measurements in morning versus afternoon at a given soil collar were taken with different instruments.Consistent with the early work at this site (Lundegard and Johnson 2006), decane was used as the representative hydrocarbon when converting CO 2 effluxes to hydrocarbon removal rates.
One groundwater monitoring well within each CO 2 efflux grid was fitted for subsurface temperature measurements.Each well was 4-in.diameter and extended below the current water table.In each well, a 2-in.diameter blank polyvinyl chloride (PVC) pipe with a sealed bottom cap was inserted to the total depth.Water from on-site well PAP3-1WW (screened in the Principal Aquifer, unaffected by petroleum hydrocarbons) was then poured into the PVC pipe and allowed to equilibrate with the subsurface temperature.The timeframe for this equilibration period was determined to be less than 24 h based on preliminary testing that confirmed minimal changes to the temperature profile for several days afterwards.The water was placed into the well to increase the representativeness of temperature measurements.Specifically, water has a higher viscosity than air, so the water minimized convection that might occur due to a temperature gradient in the subsurface, such as due to heating at depth (Warren and Bekins 2015).
Temperature measurements were taken with a Fish Hawk TD probe.The Fish Hawk TD has an epoxy-coated thermistor (Betacurve series GA100K6) for taking temperature measurements.The thermistor has a less than 1 s response time, a 0.1-°C tolerance, and range of 0 to 70 °C.The Fish Hawk TD also has a pressure transducer for measuring depth of the instrument below the top of water.The transducer is 15.6 mm (MAP 36A series), has a 1 millisecond response time, and a 0.43 m accuracy.The Fish Hawk TD is programmed to measure water temperature at depth incre- ments of 1.5 m.The data are stored in the probe software while the probe is lowered.Upon reaching the total depth, the measurements cease, and when the probe is retrieved, the user can cycle through each of the measurements and record the results on field forms.The probe's internal memory is then wiped when the next measurement is initiated.For each measurement, the probe was slowly lowered with a fishing line down the water-filled blank pipes.The slow downward velocity, and the small volume of the fishing line, minimized disturbance and displacement of water in the pipes while measurements were taken.Replicate measurements were taken in immediate succession, reviewed to confirm the results were indeed similar, and then values for each depth were averaged.This provided one temperature profile for each well and measurement event.Temperature peaks were identified as the maximum value for a given profile.The peaks are considered to be accurately measured at the 1.5-m depth measurement increments (i.e., 0.43 m accuracy of probe is smaller than the increments) and sufficient for thermal gradient and energy storage change calculations described below.
Temperature measurements were taken in April and October 2019.Background correction used a "thermal anomaly" (i.e., difference from background) for each location.Energy flux calculations accounted for the heat transport along the thermal gradient on the thermal anomaly profile during each measurement event.Net soil heating or cooling between events may theoretically be incorporated in the energy flux calculation (Karimi Askarani et al. 2018), but temperature profiles were relatively stable over time for a given well, so this component was not included.The energy flux for the thermal gradient was calculated from the temperature peak, typically near the water table, upwards through the unsaturated soil interval, as follows: where q H,grad is the energy flux along the thermal gradient (J m −2 s −1 or W m −2 ), K T is the thermal conductivity (J m −1 s −1 K −1 ), and ΔT (K) is the change in temperature along vertical distance Δz (meters; m) after correcting for background temperature.A K T value of 1.17 J m −1 s −1 K −1 was used based on average porosity (n = 0.411) and residual water content (S wr = 0.077) for site soils (Padre Associates, Inc. 2015) and the following empirical relationship for sandy soils (Chen 2008): The calculated thermal conductivity value is within the 0.17 to 3.98 J m −1 s −1 K −1 range reported for soil (Sweeney and Ririe 2014).The vertical interval Δz was set from the temperature peak to shallower depth across which the temperature gradient was approximately linear, such that R 2 was optimized on a linear fit.All R 2 values for the linear fit were 0.97 or greater.The Δz values ranged from 4.2 to 9.1 m.This is at least one order of magnitude greater than the 0.43 m measurement accuracy on the temperature probe, so the calculations are considered representative of the subsurface gradient depth interval.Ideally, the thermal gradient would be calculated not just for upward heat transport from the temperature peak, but also downwards.Because the site wells generally have a small water column (and therefore limited depth below the temperature peak), the thermal gradient was not characterized for downward energy flux.Temperature profiles are typically approximately symmetrical immediately above and below the temperature peak (see examples in Sweeney and Ririe 2014;Warren and Bekins 2015) and expected to be so in this dune sand environment.Therefore, the upward energy flux was doubled to account for the downward energy flux due to heat transport along that thermal gradient.
The energy flux values were converted to hydrocarbon biodegradation rates by the following equation (Warren and Bekins 2015;Karimi Askarani et al. 2018): where ΔH is the enthalpy of formation for aerobic CH 4 oxidation (865.4 kJ per mole CH 4 ).As for the CO 2 efflux measurements and consistent with previous work at the site (Lundegard and Johnson 2006), decane (C 10 H 22 ) was assumed as the representative hydrocarbon molecule.Also consistent with the CO 2 efflux measurements, J HC is in units of g HC m −2 day −1 .

Enhanced Hydrocarbon Removal in Bioventing Area
In July and August 2021, a solar-powered bioventing system was installed at the KI11-3 grid (see Figure 1).One biovent well was installed in the testing area, installed by a Rhino M5T hollow stem auger rig with 8-in.diameter augers.Hydrocarbon odor and elevated photoionization detector (PID) readings were first encountered at 24 m below-ground surface (bgs).These hydrocarbon indicators were observed continuously to the water table at 26 m bgs and below this to the total borehole depth of 27 m bgs.
The biovent well was constructed with solid 4 in.diameterSchedule 40 PVC pipe from ground surface to 23.5 m bgs.The well screen was slotted PVC and set from approximately 23.5 to 27 m bgs.The bottom of the screen was set below the current water table, as a contingency to help address hydrocarbon mass exposed to unsaturated soil conditions during potential future falling water table periods.
A bottom PVC cap was placed below the screen.The borehole annulus around the well screen and above it to approximately 23 m bgs was filled with #3 Monterey sand.From 23 to 22 m bgs, bentonite chips were placed and hydrated, and from 22 to 2 m bgs, the annulus was filled with Portland cement as a grout.From 2 m bgs to ground surface, the borehole was filled with native sand.
The power supply for the biovent system was established as a set of seven 400-W photovoltaic solar panels.Solar power was used to decrease the remediation environmental footprint (ASTM International 2013) by adding an energy source outside of the California electricity grid and eliminating the need for new power lines that cross the sensitive dune habitat.The solar panels were affixed to a mounting rack that was secured to ground surface with ver- (2) K T = 7.5 1−n 0.61 n (1 − 0.0022)S wr + 0.0022 0.78n (3) tical posts and stabilized with sand bags.The solar panels were oriented to maximize capture of solar energy flux at the site location.Wiring from the solar panels was connected to a DC-to-AC inverter and variable frequency drive (VFD) controller (Pico Cell 2000, Thunderbird Solar Supply, Frederick, Colorado).An AC power cord was connected from the controller to the blower.A regenerative blower (Republic Manufacturing KP4RC210-H16) with a 0.6 kW rating was used for generation of slightly compressed air to the biovent well.The target injection rate was 0.85 cubic meters per minute (30 cfm).The blower was placed on a platform within an enclosure for protection from weather and to minimize noise.
The biovent blower outlet was piped to the biovent well using 3.2 to 5.1 cm (1.25 to 2 in.) diameter Schedule 40 PVC pipe and connected to a pressure/temperature datalogger (Monarch 5396-0338) and flow meter.This biovent outlet pipe transitioned to 5.1 cm (2 in.) diameter LDPE and was fitted with a packer (Baski) that blocked the top 1 m of screen so that the deep vadose zone was targeted by air injection.Sample tubing was run from ground surface through the packer, and fitted at the bottom with an inverted funnel.The funnel was lowered below the packer bottom during biovent operations and was raised to seal against the bottom of the packer so that gases below the packer could be isolated and sampled when in situ respiration testing was conducting during temporary shutdown of the biovent well.
Biovent monitoring infrastructure was installed in radial orientation and even spacing (Figure 2).CO 2 efflux soil collars were installed at 7.6 m incremental distances and measurements followed methods described above.Subsurface temperature monitoring locations were installed as blank 5.1 cm diameter PVC pipes at 7.6 m incremental distances, with the exception of KI11-3 which was an existing monitoring well converted to a temperature monitoring point.Temperature data processing followed methods described above, with an additional calculation of the heating between measurement events (Karimi Askarani et al. 2018).
Nested vapor wells were installed starting at 7.6 m from the biovent well and at 15.2 m incremental distances thereafter.Similar to nested vapor wells used in the 2006 NSZD assessment (Lundegard and Johnson 2006), the soilgas probes were set with relatively tight vertical spacing immediately above the water table and wider spacing in the shallower portions of the vadose zone.Soil-gas probes were installed as 15 cm long stainless steel screens and fitted with a weighted drive point on the bottom.Nylaflow tubing (0.64 cm outer diameter) was connected to the top of each probe and strung to ground surface.Flexible silicone tubing with a nylon insert plug was fitted at the top ends of the Nylaflow tubing for easy interface with sampling equipment.The borehole annulus across the vapor probes was filled with #3 Monterey sand.The borehole immediately above and below these probes was sealed with benseal and a bentonite grout.For wide (3 m or more) zones between well screens, the intervals between seals were filled with native sand.Sampling of the nested vapor wells included field measurements of O 2 , CO 2 , CH 4 , and volatile organic compounds; samples were extracted in Tedlar bags using a lung box.The pre-pilot samples were collected under diffusion conditions that are assumed for the "concentration gradient method" (Johnson et al. 2006): (4) where J HC is the NSZD rate (g HC m −2 day −1 ), S i is a stoichiometric conversion coefficient for gas i (g HC per g gas i), D eff,i is the vapor diffusion coefficient for gas i (m 2 day −1 ), and ΔC i /Δz is the concentration gradient for gas i over vertical interval z (g m −4 ).The selected gas for the concentration gradient method was CO 2 , to align with the CO 2 efflux measurements.The S i value for CO 2 is 0.32 g HC per gram CO 2 for conversion from decane.The vapor diffusion coefficient across the vadose zone was measured by helium injection to each sample port and periodic measurement as suggested in (Johnson et al. 1998).
Biovent operations were conducted from October 2021 through May 2022, with one approximate 3-week temporary shutdown in late January/early February 2022 due to an electrical short circuit of an activator switch.In May 2022, the system was intentionally shut down for an in situ respiration test that tracked decreasing O 2 concentrations versus time (EPA 1995).

Sitewide NSZD Characterization
CO 2 effluxes for all locations and across both April and October 2019, prior to the biovent pilot test, ranged from an average of 0.4 to 2.1 μmol CO 2 m −2 s −1 (Table 1 and Figure 3).Background CO 2 effluxes at the H-8 and ITA-5 locations were higher in April 2019 (average = 0.8 to 0.9 μmol CO 2 m −2 s −1 , SE = 0.05) than October 2019 (average = 0.4 μmol CO 2 m −2 s −1 , SE = 0.02 to 0.06).Time lags in the range of 1 month to multiple months for hydrocarbon-related CO 2 vadose zone transport have been observed elsewhere (Sihota et al. 2016;McAlexander and Sihota 2019).For these relatively permeable sandy soils, the elevated April CO 2 effluxes are consistent with higher precipitation in the months preceding the spring and higher natural soil respiration.Average CO 2 effluxes at the LNAPL bodies ranged from 1.5 to 2.1 (±0.06 to 0.14) μmol CO 2 m −2 s −1 in April 2019 and 0.9 to 1.7 (±0.05 to 0.16) μmol CO 2 m −2 s −1 in October 2019.All LNAPL body CO 2 effluxes were significantly higher than background (Wilcoxon rank sum test, α = 0.05, sample size for each LNAPL body and background = 25).The higher values in April versus October were essentially attributable to higher natural soil respiration, such that the NSZD rates by background correction were similar, generally in the 0.6 to 1.6 g HC m −2 day −1 range (Table 2).CO 2 effluxes above the dissolved phase plumes ranged from 1.0 to 1.1 μmol CO 2 m −2 s −1 in April 2019 and were 0.5 μmol CO 2 m −2 s −1 at both measurement grids (206R-D and I6-1) in October 2019.The dissolved phase plume (206R-D and I6-1) CO 2 effluxes were not significantly different from background (Wilcoxon rank sum test, α = 0.05; sample size for each dissolved phase plume and background = 25).
Subsurface temperature measurements in the background areas were stable at depth (Figure 4), with the ITA-5 location reporting slightly higher values than the H8-1 area.The variation might be due to the ITA-5 being located adjacent to a hillside with a high surface area for contact with the atmosphere.The H8-1 location topography is more indicative of the other monitored locations, so it was used as background for temperature anomaly calculations.All six temperature profiles above LNAPL bodies demonstrated a heated interval in the deep vadose zone, likely evidence of hydrocarbon or CH 4 oxidation by aerobic microorganism biodegradation reactions (Warren and Bekins 2015).This was indicated as temperatures consistently higher than background.Converted to hydrocarbon removal via equation (3) above, the NSZD rates ranged from 0.4 to 1.7 g HC m −2 day −1 (Table 2).These values showed minimal variability between the April and October 2019 measurement events, and were similar to those estimated by CO 2 efflux measurements (Table 2).Given the close agreement between rate estimates, the temperature data support the conceptual model that elevated CO 2 efflux at ground surface is due to hydrocarbon/CH 4 oxidation in the deep vadose zone, followed by CO 2 transport upwards.The temperature profile for the I6-1 dissolved phase plume was not significantly different from background (Wilcoxon pairwise test matched to depth, α = 0.05; sample size for each profile = 12), but the temperature profile for the 208R-D dissolved phase plume was significant (sample size for each profile = 13).The tem- Values in parentheses are one standard error of the mean.Sample size is 25 measurements for CO 2 effluxes at each location and varies from 5 to 23 measurements for temperature, depending on depth of water table (see Figure 4).perature peak at this location was in the shallow vadose zone (3 m) and less than 1.0 °C (Figure 4).Hydrocarbon biodegradation cannot be ruled out for the temperature anomaly at this location, but the shallow peak suggests a weaker heating source not associated with contaminant attenuation.The 2006 soil-gas profiling reported NSZD rates of 0.4 to 3.0 g HC m −2 day −1 (Lundegard and Johnson 2006) above the main LNAPL bodies.This is comparable to the 0.6 to 1.6 g HC m −2 day −1 range for the current study's CO 2 efflux measurements over these same LNAPL bodies.
The wider range in NSZD rate estimates for the 2006 work is attributable to uncertainty in the actual vapor diffusion coefficient at any given soil-gas sampling location due to potential soil moisture variation.The CO 2 efflux measurements are a direct measure of gas transport over a gridded area such that an average value for the footprint can be calculated.Thus, the CO 2 efflux's "error" term is due to actual measured variability rather than uncertainty for an input value.The NSZD rate above LNAPL bodies (0.4 to 1.7 g HC m −2 day −1 ) for subsurface temperature measurements was also narrower than the 2006 soil-gas profiling.While the subsurface measurements require input parameters for the calculation (e.g., soil thermal conductivity and heat capacity), they are generally less sensitive to variations in soil moisture than the vapor diffusion coefficient.Soilgas profiling provides information about subsurface transport processes (e.g., whether methanogenesis is significant), but this study demonstrates the tighter estimates on NSZD rates that can be achieved by the other two methods.
The stability of NSZD rates 6 months apart contrasts with seasonal variation observed at sites in the middle portion of the U.S. (Sihota et al. 2016;Eichert et al. 2017;McAlexander and Sihota 2019) and aligns with steady rates observed in southern California (Mackay et al. 2018a).The primary weather component that varies at the site is precipitation, with rainfall (30 cm, on average; 45 cm in 2019) typically in the winter months and dry or foggy conditions the remainder of the year.CO 2 efflux measurements were taken Dissolved phase plumes Diluent bodies Background when shallow soils were dry to provide data representative of prevailing conditions, avoiding short term gas transport limited conditions.NSZD rates have been correlated with precipitation rates at a semi-arid site (McAlexander and Sihota 2019), but the depth to groundwater was less than 5 m.At the current site, the water table was deeper than 15 m at all studied locations (Figure 4).The large transport distance for percolating water to reach the microbially active zone may have dampened a detectable seasonal effect.Alternatively, with the studied portions of the site only containing LNAPL near the water table where moisture supply may be sufficient regardless of weather conditions, influence from precipitation on NSZD may be minimal.The 2019 NSZD rates were within the range of values reported for the 2006 event (Lundegard and Johnson 2006).Given the different data collection method for that initial effort, and the approximate order-of-magnitude range in values, it is not possible to assess the trend in NSZD rate over time based on the two datasets.Hydrocarbon removal from an LNAPL body is continually evolving as individual hydrocarbons are removed by different mechanisms (Essaid et al. 2011).Davis et al. (2022) observed generally decreasing NSZD rates for a large petroleum refinery over approximately 2 decades.Both the CO 2 effluxes and subsurface temperature measurements provided narrow NSZD rate estimates despite characterizing LNAPL bodies over a wide footprint, suggesting promise for long-term tracking.

Enhanced Hydrocarbon Removal in Bioventing Area
Air injection at the solar-powered biovent pilot test area had demonstrable effects on soil-gas profiles (Figure 5).At baseline, the profiles aligned with previous sitewide characterization (Lundegard and Johnson 2006), with depleted O 2 near the water table.By November 2021 when the biovent system had been operating for approximately 1 month, O 2 concentrations greater than 20% were observed near the water table 7.6 m from the injection well, and elevated O 2 (greater than 5%) was observed at a 22.9 m distance approximately 6.1 m above the water table.The O 2 profiles indicated a potential fanning of injected air upwards toward ground surface.The high O 2 conditions were observed in all other months when the system was operating, demonstrating persistence of the elevated O 2 in the deep vadose zone.For the February 2022 event, when the system was not operating, the soil-gas conditions were intermediate between baseline and operating conditions.This indicates a detectable response in terms of biogenic O 2 consumption during temporary suspension of bioventing operations.O 2 utilization to support enhanced biodegradation was observed by development of a deep vadose zone heated above baseline conditions (Figure 6).At baseline, the maximum temperature immediately above the water table was 21.4 °C.By the final monitoring event more than 200 days later, the maximum temperature in this zone had increased to 24.1 °C, while the background and reference locations had stable subsurface temperatures (Supplemental Information).For the late January 2022 measurement event (panel d of Figure 6) when the biovent system was not operating due to the weather-induced shutdown, the subsurface temperature remained elevated.This suggests that the temperature parameter is relatively insensitive to short-term changes to bioventing operations.For all events, the temperature peak was observed immediately above the water table.The lateral distance of heating was 30 m, providing an estimate of the radius of influence (ROI) for enhanced source zone depletion.
The increased O 2 utilization in the deep vadose zone generated an elevated CO 2 efflux signal at ground surface (Figure 7).For all events when the biovent system was operating, median biovent CO 2 effluxes were significantly greater than average control area CO 2 effluxes (Wilcoxon rank sum test at α = 0.05; treatment n = 26, control n = 18).For the January 2022 event when the biovent system was not operating, median biovent and control area CO 2 effluxes were statistically similar.Thus, similar to subsurface O 2 concentrations, CO 2 effluxes exhibited a rapid response to temporary suspension of biovent operations.The lateral distance of elevated CO 2 effluxes was at least 38 m (Figure 8).However, the 30 m ROI estimated by subsurface temperatures is favored as more representative of the actual subsurface biodegradation zone, as soil-gas data (Figure 5) suggested upward fanning of injected air that could have resulted in a wider apparent footprint at ground surface.While the CO 2 efflux measurement was demonstrated as robust for hydrocarbon removal rate estimates, it is not as useful for estimating the subsurface ROI under advective vapor transport conditions because the gases likely spread laterally before reaching ground surface.
The baseline soil-gas profiles can be used to estimate NSZD rates in the pilot test and control areas under baseline conditions prior to advection induced by the system operation.The baseline NSZD rates ranged from 0.2 to 0.3 g HC m −2 day −1 , comparing well with the 0.2 to 0.7 g HC m −2 day −1 range for temperature profiles and 0.2 to 0.6 g HC m −2 day −1 range for CO 2 effluxes (Table 3).Similar to the sitewide NSZD rates reported in the previous section, these values were consistent with those measured by the earlier concentration gradient method (Lundegard and Johnson 2006).This demonstrates that all three methods provide similar hydrocarbon removal rate estimates and confirms the conceptual model that each of them reasonably quantifies hydrocarbon biodegradation reactions occurring in the vadose zone and that all methods are suitable for informing site management decisions.
During treatment, estimated mean NSZD rates for the control area, where only natural hydrocarbon degradation occurred, were steadier for temperature profiles (0.7 g HC m −2 day −1 for baseline versus 0.9 g HC m −2 day −1 for treatment; Table 3) than for CO 2 effluxes (increase from 0.3 to 0.6 g HC m −2 day −1 ).This may indicate the higher responsiveness of CO 2 efflux data to detecting subtle changes in NSZD rates over time.During treatment, source depletion  rates were considerably higher than baseline in the pilot test area, including higher than the 0.3 g HC m −2 day −1 change observed in the control area for CO 2 effluxes.Temperature profiles estimated 1.5 to 5.4 g HC m −2 day −1 for distances from 0 m to 30.4 m from the biovent well (see treatment values for temperature profile on Table 3), with an average value of 3.5 g HC m −2 day −1 .At 38.0 m from the biovent well, the estimated source depletion rate for temperature profiles was 0.6 g HC m −2 day −1 , which was similar to baseline and indicative of a smaller (nominally 30.4 m) ROI.Estimated CO 2 effluxes were 1.4 to 2.8 g HC m −2 day −1 from 0 to 38.0 m from the biovent well (see Treatment values for CO 2 effluxes on Table 3).Across this 38.0 m inferred ROI, the CO 2 efflux average value during treatment was 1.8 g HC m −2 day −1 .Integrating over the estimated ROIs to obtain total mass removal rates, the temperature profiles suggests an estimate of 10.2 kg day −1 and the CO 2 effluxes suggest 8.2 kg day −1 above the NSZD total mass removal rate for these footprints.These similar values support the inference of soil-gas fanning upwards, such that the 30 m ROI in the deep vadose zone where the reaction occurred was spread over a 38 m ROI at ground surface.
In-situ respiration testing conducted upon shutdown of the system in May 2022 indicates linear decreasing O 2 con-centration with time, mirrored by increasing CO 2 concentration (Figure 9).This further corroborates the inference of O 2 consumption for hydrocarbon biodegradation to generate CO 2 .The 0.86% day −1 O 2 utilization rate corresponds to reaching O 2 limited conditions (5%; EPA 1995) at Days 18 to 19 after system shutdown.Based on soil-gas data, the biovent system required approximately 3 days to ramp up from NSZD (low O 2 ) conditions to atmosphere-level oxygenated conditions in the deep vadose zone, so the 0.85 cubic meter per minute (30 cfm) blower is oversized for one well.On a full-scale application design, one blower with this capacity could be used more efficiently by circulating air injection to multiple wells, each receiving a finite period of air injection before moving to the next well, such that O 2 concentrations are held above 5% over a wide footprint.

Summary
Following mechanical remediation at LNAPL sites, residual hydrocarbons may remain in the subsurface, necessitating long-term management.While much of the effort in recent years has focused solely on NSZD as the primary method for demonstrating continued hydrocarbon removal, the same data collection methods can quantify sustainable, Treatment CO 2 efflux events are for "System On" timeframes shown on Figure 8. Units are g HC m −2 day −1 .Values in parentheses are one standard error of the mean.low energy consumption methods to enhance biodegradation.This can provide a bridge between active mechanical mass removal systems and NSZD.
The first objective of this study was to measure NSZD rates broadly across the site over footprints that could not be achieved with the soil-gas monitoring method.Both subsurface temperature profiling and ground surface CO 2 effluxes were demonstrated to compare well with the early soil-gas concentration gradient method for estimating NSZD rates.Due to shorter times required for data collection, the temperature and CO 2 efflux methods provide a larger number of measurements at a reduced time and cost compared to the soil-gas gradient method and, consequently, provided greater statistical strength.
The NSZD rate estimates at the Site for the current study used these two newer methods, and corroborated the results reported in 2006 for the soil-gas gradient method.The NSZD rate estimates by both CO 2 effluxes (0.6 to 1.6 g HC m −2 day −1 ) and subsurface temperature measurements (0.4 to 1.7 g HC m −2 day −1 ) were similar to but in a narrower range than the earlier estimates (0.4 to 3.0 g HC m −2 day −1 ).The wide range for the 2006 estimates would render a parametric statistical hypothesis test comparison with the 2019 rates as uninformative.The 2019 results do indicate minimal seasonal variability for that year.All three methods are now demonstrated as viable for tracking NSZD, and the strengths of each method (ASTM International 2022) can be considered when designing a monitoring program.For instance, where a large footprint of measurements is desired, the CO 2 efflux measurements can be deployed.Where temporal variation is suspected, subsurface temperature measurements can be taken.Where knowledge about the relative contribution of methanogenesis and aerobic biodegradation is needed, the soil-gas gradient method can be used.Over time periods of multiple decades, the NSZD rate may decrease as the quantity, and potentially composition, of the remaining LNAPL changes (Davis et al. 2022), and now this site has multiple tools to track this progression.
This study demonstrated that dissolved phase plumes do not provide strong hydrocarbon removal signals for CO 2 efflux and subsurface temperature measurements.There are several potentially relevant factors.First, LNAPL bodies have high quantities of hydrocarbons undergoing biodegradation, imparting significant demand for electron acceptors and having methanogenesis as a dominant attenuation mechanism (Lundegard and Johnson 2006;Essaid et al. 2011).The methane is ultimately oxidized in the vadose zone, generating significant quantities of CO 2 and heat, both of which are detected by the methods employed in this study.Second, the dissolved phase plumes at this site over the monitoring timeframe have had low concentrations of benzene, toluene, ethylbenzene, and xylenes and are detected generally in the semivolatile ranges (C 10 to C 25 ) hydrocarbon mass range.These higher molecular weight compounds tend to have low volatility, limiting one attenuation mechanism.Third, silica gel cleanup on groundwater samples in the dissolved phase plumes prior to analysis tends to remove a large fraction of the C 10 -C 25 compounds, indicating that the primary components of the dissolved phase "TPH" plumes are oxygencontaining organic compounds and likely metabolites from hydrocarbon biodegradation (Zemo et al. 2013).These polar compounds tend to have high solubility and low volatility, resulting in minimal generation of CO 2 or heat in the vadose zone upon further biodegradation.
The second and third objectives of the study were to direct the new NSZD methods to quantify enhanced biodegradation of hydrocarbons by solar-powered bioventing, and to apply these learnings to LNAPL mass removal management for the site, respectively.A 0.85 cubic meter per minute (30 cfm) blower connected to a 2800 W photovoltaic supply achieved enhanced source zone depletion rates as evidenced by all three of the NSZD data collection methods.Soil-gas data indicated an upward fanning of injected air from the deep vadose zone to ground surface.This was corroborated by elevated CO 2 effluxes over a wider footprint than the identified deep treatment zone.Elevated temperature due to hydrocarbon oxidation yielded an approximate total mass removal rate of 10.2 kg day −1 source depletion rate above the NSZD total mass removal rate for a corresponding 0.29 ha footprint (30.4 m ROI), which aligned well with a 8.2 kg day −1 rate estimated via CO 2 efflux measurements for a 0.45 ha footprint that accounts for this biodegradation product gas spreading as it migrated upward to ground surface (i.e., 38 m ROI).The baseline NSZD rates for the pilot test area ranged from 0.2 to 0.7 g HC m −2 day −1 , which for a 0.29 ha footprint (i.e., at the source zone prior to gas spreading) corresponds to 0.6 to 2.0 kg day −1 .The bioventing source depletion rate was a substantial increase over baseline NSZD, approaching an order of magnitude.Separate evaluation of the active LNAPL recovery system at the site has determined that wells with low LNAPL transmissivity generally recover less than 5 kg day -1 each (Trihydro 2019).This is lower than the hydrocarbon removal rate achieved by bioventing for a comparable footprint.This indicates that areas with large footprints of LNAPL presence, but low LNAPL transmissivity, can benefit from a transition from hydraulic recovery to this active "green" remediation strategy (Smith et al. 2023) that seeks to minimize greenhouse gas emissions.The solar-powered biovent system also avoids new power lines that cross the sensitive dune environment.Endpoints for bioventing can be established as the NSZD rate, such that those operations can transition to NSZD without enhancement once the active remediation is providing minimal additional benefit.Site management can proceed along a sequenced program that starts with aggressive hydraulic recovery to address mobile recoverable LNAPL, transitions to engineered/enhanced NSZD in areas with poorly recoverable LNAPL, and ultimately transitions to NSZD without intervention across the full footprint of the LNAPL bodies.
At this site, water table elevation changes occur over relatively long periods, such that the influence on NSZD rates could not be assessed.Theoretically, a higher water table elevation could dampen vapor transport in the shallower portions of the smear zone, affecting hydrocarbon loss rates under both natural and induced high O 2 conditions and the ultimate remedial timeframe.However, such a theoretical condition has not been demonstrated elsewhere (McAlexander and Sihota 2019), possibly because initial methanogenesis within the low-O 2 smear zone is dominant (Ng et al. 2015) and occurs both above and below the water table.At this and other sites, such evaluations can be conducted during full scale enhanced SZD operations, as long as concurrent NSZD monitoring is conducted in untreated areas.Thus, remedy implementation can proceed once pilot testing establishes a benefit for enhanced SZD, and identification of controlling environmental factors can be incorporated to long-term performance monitoring.
Similar to many LNAPL affected sites, initial LNAPL recovery operations have successfully transitioned subsurface conditions to low LNAPL transmissivity in many locations.NSZD is one possible next step following such aggressive treatment.This study builds on the body of evidence that NSZD occurs over large footprints, and that data collection methods can be demonstrated as viable on a site-specific basis.However, these same methods also demonstrated solar-powered bioventing as one viable "green remediation" approach that can be used, if warranted, before transitioning to NSZD.The substantial, quantified increase in source depletion rates at this site indicates benefit for using enhancements to natural processes and continued engineered remediation for site management.

Figure 1 .
Figure 1.Site overview, showing the distribution of delineated LNAPL and dissolved total petroleum hydrocarbon (TPH; carbon range C 10 to C 25 ) plumes across the site, and the locations of NSZD study areas.

Figure 3 .
Figure 3. CO 2 effluxes for sitewide NSZD assessment in (a) April 2019 and (b) October 2019.Horizontal dashed line indicates average background value (H8-1).CO 2 effluxes for all three diluent bodies (six grids) were significantly different from background (α = 0.05) by Wilcoxon rank sum test.CO 2 effluxes for both dissolved phase plumes were not significantly different from background.

Figure 4 .
Figure 4. Temperature profiles for sitewide NSZD assessment (April 2019), with "BG" indicating background locations.(a) The two background plots, (b) TB9 and CP area profiles versus background, (c) DT area profiles versus background, and (d) dissolved phase plume areas versus background.All displayed depths are below-ground surface.The water table depth varies across the site and is near the bottom of each profile.

Figure 5 .Figure 6 .
Figure 5. Soil-gas O 2 profiles for bioventing area.Measurement points are displayed as black circles.Top of vadose zone is 27 m above water table, 3 m above the shallowest sample point.Biovent screen above water table is displayed as a gray rectangle at 0 m on each x axis.

Figure 8 .
Figure 8. CO 2 effluxes by distance for bioventing area.Plotted values are converted to NSZD rates in Table3.

Table 2 Average NSZD Rates for LNAPL Bodies Grid April 2019 October 2019 CO 2 Efflux Temperature CO 2 Efflux Temperature
Units are g HC m −2 day −1 .Values in parentheses are one standard error of the mean (i.e., for 25 measurements of CO 2 efflux in a 5 × 5 grid pattern corrected to an average background value).Subsurface temperatures do not have standard errors as the full temperature profile is used to estimate one NSZD rate.