Incorporating the Soil Gas Gradient Method and Functional Genes to Assess the Natural Source Zone Depletion at a Petroleum-Hydrocarbon-Contaminated Site of a Purification Plant in Northwest China

An increasing number of studies have demonstrated that natural source zone depletion (NSZD) in the vadose zone accounts for the majority (90%~99%) of the natural attenuation of light non-aqueous phase liquid (LNAPL). Until now, 0.05 to 12 kg/a.m2 NSZD rates at tens of petroleum LNAPL source zones have been determined in the middle or late evolution stage of LNAPL release, in which limited volatile organic compounds (VOCs) and methane (CH4) were detected. NSZD rates are normally estimated by the gradient method, yet the associated functional microbial activity remains poorly investigated. Herein, the NSZD at an LNAPL-releasing site was studied using both soil gas gradient methods quantifying the O2, CO2, CH4, and VOCs concentrations and molecular biology methods quantifying the abundance of the pmoA and mcrA genes. The results showed that the methanogenesis rates were around 4 to 40 kg/a.m2. The values were greater than the rates calculated by the sum of CH4 escaping (0.3~1.2 kg/a.m2) and O2 consuming (3~13 kg/a.m2) or CO2 generating rates (2~4 kg/a.m2), suggesting that the generated CH4 was oxidized but not thoroughly to CO2. The functional gene quantification also supported the indication of this process. Therefore, the NSZD rates at the site roughly equaled the methanogenesis rates (4~40 kg/a.m2), which were greater than most of the previously studied sites with a 90th percentile value of 4 kg/a.m2. The study extended the current knowledge of the NSZD and has significant implications for LNAPL remediation management.


Introduction
Natural source zone depletion (NSZD) is the natural loss of light non-aqueous phase liquid (LNAPL) through collective, naturally occurring processes of volatilization, dissolution, and biodegradation [1][2][3]. A growing number of studies have demonstrated that NSZD occurs at petroleum-hydrocarbon-contaminated sites at depletion rates ranging from thousands to tens of thousands of liters per hectare per year. Approximately 90% to 99% of the natural attenuation of petroleum hydrocarbons occurs in the vadose zone, the area extending from the surface to the regional groundwater table [4,5]. Across the petroleum-hydrocarbon-contaminated aquifer to the vadose zone, three zones, i.e., the Figure 1. The characteristics of the study area and LNAPL contamination. Plot (a) shows the sampling location, groundwater flow direction, and the TPH contaminant concentrations (as represented by the concentration contour plot). The red rectangle near the storage tanks (i.e., block A, B, and C) represents the area where NSZD was evaluated. The three red triangles in the non-contaminated area, north of the storage tanks, represent the background sampling points. Plot (b) shows the vertical distribution of the LNAPLs and contamination plume.

Sampling, Gas Measurement, and Functional Gene Determination
The most likely contamination source zone to the north of the in-service storage tanks was chosen to assess the rate of natural source zone depletion. Regarding the complex pipeline underground and the presence of flammable and explosive hydrocarbons, it was a forbidden area for traditional invasive drilling and long-term retention of monitoring wells according to the strict rules of the plant. To gain the NSZD rate and data on the different physicochemical pathways of NSZD, a minimally invasive measurement was carried out in the study area. The area was divided into 3 one-meter square blocks ( Figure  1a). All blocks were treated as replicates. The concentrations of soil gases and functional genes were measured in each block. Meanwhile, the non-contaminated area was chosen as the background.
(1) Soil gas concentration measurement Using a 2-cm-diameter, 1.0-m-long copper drill rod hammered with a copper hammer, boreholes with different depths were drilled. First, a 10-cm-deep hole was drilled. The concentrations of VOCs, CO2, O2, H2, and CH4 were immediately measured in the drilled hole using a portable multi-parameter gas detector (MultiRAE 6208, USA). Afterward, in the same way, the hole was deepened to 20 cm and the gas concentrations were measured. The same procedure was applied to a depth of 90 cm. Triplicate measurements were conducted in each block.
Due to the fact that the gas measured in the borehole was a mixture of gas extracted from different soil depths, the measured concentrations in each layer might not fully represent their real concentrations. In the current study, to obtain the real concentration in each layer, we proposed a procedure to estimate the real gas concentrations in a certain depth layer.
The assumptions of the measurement were i) the soil column is homogeneous within 90 cm, i.e., the gas permeabilities are uniform; ii) at the same depth, the gas concentrations are the same anywhere; iii) during the course of gas extraction and measurement, the gas is equally extracted from different depths, i.e., soil gases from different depths have the same contribution; iv) when the measured gas concentration values are stabilized, the stabilized values represent the soil gas concentrations. Figure 1. The characteristics of the study area and LNAPL contamination. Plot (a) shows the sampling location, groundwater flow direction, and the TPH contaminant concentrations (as represented by the concentration contour plot). The red rectangle near the storage tanks (i.e., block A, B, and C) represents the area where NSZD was evaluated. The three red triangles in the non-contaminated area, north of the storage tanks, represent the background sampling points. Plot (b) shows the vertical distribution of the LNAPLs and contamination plume.

Sampling, Gas Measurement, and Functional Gene Determination
The most likely contamination source zone to the north of the in-service storage tanks was chosen to assess the rate of natural source zone depletion. Regarding the complex pipeline underground and the presence of flammable and explosive hydrocarbons, it was a forbidden area for traditional invasive drilling and long-term retention of monitoring wells according to the strict rules of the plant. To gain the NSZD rate and data on the different physicochemical pathways of NSZD, a minimally invasive measurement was carried out in the study area. The area was divided into 3 one-meter square blocks ( Figure 1a). All blocks were treated as replicates. The concentrations of soil gases and functional genes were measured in each block. Meanwhile, the non-contaminated area was chosen as the background.
(1) Soil gas concentration measurement Using a 2-cm-diameter, 1.0-m-long copper drill rod hammered with a copper hammer, boreholes with different depths were drilled. First, a 10-cm-deep hole was drilled. The concentrations of VOCs, CO 2 , O 2 , H 2 , and CH 4 were immediately measured in the drilled hole using a portable multi-parameter gas detector (MultiRAE 6208, USA). Afterward, in the same way, the hole was deepened to 20 cm and the gas concentrations were measured. The same procedure was applied to a depth of 90 cm. Triplicate measurements were conducted in each block.
Due to the fact that the gas measured in the borehole was a mixture of gas extracted from different soil depths, the measured concentrations in each layer might not fully represent their real concentrations. In the current study, to obtain the real concentration in each layer, we proposed a procedure to estimate the real gas concentrations in a certain depth layer.
The assumptions of the measurement were (i) the soil column is homogeneous within 90 cm, i.e., the gas permeabilities are uniform; (ii) at the same depth, the gas concentrations are the same anywhere; (iii) during the course of gas extraction and measurement, the gas is equally extracted from different depths, i.e., soil gases from different depths have the same contribution; (iv) when the measured gas concentration values are stabilized, the stabilized values represent the soil gas concentrations. Then, the amount of a certain gas in the layer between (x − d) and x depths ( Figure 2) can be expressed as where x is the depth below the ground surface, m; C (x) is the actual gas concentration at x depth, g/m 3 ; is the amount of the gas in the layer between x − d and x depths, g/m 2 . Assuming the measured concentration in the x-depth hole is C m (x), the amount of the gas in the hole can be expressed as The amount of the gas in the x-depth hole can be deemed as the amount of the gas in the x − d depth hole plus the amount of the gas in the layer between the x − d and x depths, and can therefore be expressed as Therefore, For assessing NSZD in practice, when d is small enough, between x − d and x depths, C a (x) can be deemed as a constant value. Then, Equation (4) can be transformed as C m (x) and C m (x − d) in Equation (5) are easily measured. Then, the actual gas concentration at x depth can be calculated as In the present study, d was set as 0.1 m.
Then, the amount of a certain gas in the layer between (x − d) and x depths (Figure 2 can be expressed as where x is the depth below the ground surface, m; C (x) is the actual gas concentration at x depth, g/m 3 ; is the amount of the gas in the layer between x − d and x depths, g/m 2 . Assuming the measured concentration in the x-depth hole is Cm(x), the amount of th gas in the hole can be expressed as The amount of the gas in the x-depth hole can be deemed as the amount of the gas in the x − d depth hole plus the amount of the gas in the layer between the x − d and x depths and can therefore be expressed as For assessing NSZD in practice, when d is small enough, between x − d and x depths Ca(x) can be deemed as a constant value. Then, Equation (4) can be transformed as Cm(x) and Cm(x − d) in Equation (5) are easily measured. Then, the actual gas concen tration at x depth can be calculated as In the present study, d was set as 0.1 m. (2) Functional gene determination To verify the crucial biochemical processes upon NSZD, e.g., methanogenesis and methane oxidation, the corresponding functional genes methyl coenzyme M reductase gene (mcrA) and methane monooxygenase gene (pmoA) were quantified. The detailed procedure was as follows: At each depth of the gas concentration measurement, an approximately 50 g soil sample was collected using the Luoyang shovel and sealed in a sterile bag. Then, the soil samples were stored in an incubator filled with dry ice and transported to the laboratory for DNA extraction. The DNA was extracted from 0.8 g soil of each sample using a FastDNA kit (Q-BIO gene Corp. Irvine, CA, USA) [21]. Then, the mcrA and pmoA gene abundance was quantified using a fluorescence quantitative PCR instrument (ABI Q5, USA). The PCR primers and detailed experimental procedures are available in the relevant references [22,23].

Assessment of the Natural Source Zone Depletion
After the C(x) values of O 2 , CH 4 , and VOCs were gained, the NSZD rate could be calculated according to Equation [19] represent the vertical concentration gradients of hydrocarbon, methane, and oxygen, respectively, at depth (x) of the horizontal plane in (g/m 3 )/m; S CH4 is the stoichiometric coefficient for methanogenesis, and the value is 1.1 g-HC/g-CH 4 ; S O2 is the stoichiometric coefficient for aerobic biodegradation, ranging from approximately 0.25 to 0.29 g hydrocarbon/g O 2 consumed, depending on the relative contributions of direct hydrocarbon aerobic oxidation (0.29 kg hydrocarbon/mg O 2 ) and indirect hydrocarbon oxidation (0.25 mg hydrocarbon/mg O 2 , assuming that methane production occurs first in the anaerobic source zone and then methane is subsequently biodegraded aerobically as it diffuses upward); D HC , D CH4 , and D O2 represent the effective vapor phase diffusion coefficients for hydrocarbon, methane, and oxygen, respectively, at depth (x) of the horizontal plane in m 2 /s. Effective diffusion coefficients were estimated by the most widely used Penman model [24].
where D is the gas diffusion coefficient in soil (m 2 /s); D 0 is the gas diffusion coefficient in free air (m 2 /s); ε is the soil air-filled porosity. According to Equation (8), the effective vapor phase diffusion coefficients for hydrocarbon (represented by benzene), CH 4 , O 2 , and CO 2 were calculated to be 1.27 × 10 −4 m 2 /s, 3.08 × 10 −5 m 2 /s, 2.47 × 10 −5 m 2 /s, and 2.16 × 10 −5 m 2 /s respectively. In the estimation, the air-filled porosity was measured to be 0.2 and the gas diffusion coefficients in free air were gained from the work of Massman [25].
In this study, the CO 2 and CH 4 gradients were also used to assess the CO 2 production rate and methanogenesis rate, respectively.

Soil Gas Profiles
Based on the minimally invasive measurement of soil gas concentrations, the actual soil gas concentrations in the study area were calculated using Equation (6), and they are shown in Figure 3. It can be observed that the O 2 concentrations rapidly decreased with depth above 60 cm in all blocks, while the VOCs, CO 2 , and CH 4 concentrations increased. Below 60 cm depth, the VOCs and O 2 concentrations did not vary with depth and exhibited Life 2023, 13, 114 6 of 12 some fluctuation, while the CO 2 and CH 4 concentrations continued to increase with depth in blocks A and B, and were greater than the upper limits of detection in block C.

Soil Gas Profiles
Based on the minimally invasive measurement of soil gas concentrations, the actual soil gas concentrations in the study area were calculated using Equation (6), and they are shown in Figure 3. It can be observed that the O2 concentrations rapidly decreased with depth above 60 cm in all blocks, while the VOCs, CO2, and CH4 concentrations increased. Below 60 cm depth, the VOCs and O2 concentrations did not vary with depth and exhibited some fluctuation, while the CO2 and CH4 concentrations continued to increase with depth in blocks A and B, and were greater than the upper limits of detection in block C. Although the VOCs were detected and basically increased with depth in all blocks, the concentrations of VOCs are about 2 to 4 orders lower than CO2 and CH4. The trend of CO2 concentration variations and O2 displayed a significant negative correlation (r = −0.967, p < 0.001). The depth-dependent variations in CH4 concentration were nonlinear and had different patterns. Above the 50cm, 40cm, and 30cm for blocks A, B, and C, the CH4 concentration was less than 1% with fewer variations, while below these depths, the CH4 concentration was greater than 1% and almost linearly increased with depth. The lower CH4 concentration in the upper layers (<50 cm, 40 cm, and 30 cm for blocks A, B, and C) indicated that there might be less methanogenesis and more methane oxidation, and the O2 concentration gradients were mainly formed by diffusion.

Soil Gas Gradients
To calculate the NSZD rate, only the O2 concentrations in diffusion-dominated areas were taken into consideration. As methane oxidation also generates CO2, the CO2 diffusion-dominant areas were the same as for O2. It was hard to distinguish the CH4 and VOCs diffusion-dominated areas because CH4 and VOCs might be oxidized in the space coexisting with O2. If the oxidation rates were constant at different depths, the concentration gradients would be only caused by diffusion. Upon this assumption, to estimate the rates of VOCs volatilization, the concentrations at the same depths as O2 diffusion-dominated depths were taken into account; and to estimate the rates of CH4 generation in the source zone, the concentrations in the soil below the mutational depths (>50 cm, 40 cm, and 30 cm for blocks A, B, and C) were taken into consideration. Herein, the linear curves were fitted with the gas profile data. The fitted depths, calculated gradient, and adj. R 2 are listed in Table 1. Although the VOCs were detected and basically increased with depth in all blocks, the concentrations of VOCs are about 2 to 4 orders lower than CO 2 and CH 4 . The trend of CO 2 concentration variations and O 2 displayed a significant negative correlation (r = −0.967, p < 0.001). The depth-dependent variations in CH 4 concentration were nonlinear and had different patterns. Above the 50cm, 40cm, and 30cm for blocks A, B, and C, the CH 4 concentration was less than 1% with fewer variations, while below these depths, the CH 4 concentration was greater than 1% and almost linearly increased with depth. The lower CH 4 concentration in the upper layers (<50 cm, 40 cm, and 30 cm for blocks A, B, and C) indicated that there might be less methanogenesis and more methane oxidation, and the O 2 concentration gradients were mainly formed by diffusion.

Soil Gas Gradients
To calculate the NSZD rate, only the O 2 concentrations in diffusion-dominated areas were taken into consideration. As methane oxidation also generates CO 2 , the CO 2 diffusiondominant areas were the same as for O 2 . It was hard to distinguish the CH 4 and VOCs diffusion-dominated areas because CH 4 and VOCs might be oxidized in the space coexisting with O 2 . If the oxidation rates were constant at different depths, the concentration gradients would be only caused by diffusion. Upon this assumption, to estimate the rates of VOCs volatilization, the concentrations at the same depths as O 2 diffusion-dominated depths were taken into account; and to estimate the rates of CH 4 generation in the source zone, the concentrations in the soil below the mutational depths (>50 cm, 40 cm, and 30 cm for blocks A, B, and C) were taken into consideration. Herein, the linear curves were fitted with the gas profile data. The fitted depths, calculated gradient, and adj. R 2 are listed in Table 1. As shown by the adj. R 2 values in Table 1, the concentrations of O 2 , CO 2 , and CH 4 were well fitted, whereas the concentrations of VOCs were not.

Estimating the NSZD Rates
In agreement with the prior assumption, in the O 2 diffusion-dominated area, all soil gas concentrations only varied by diffusion, and the NSZD rate in each block could be calculated according to Equation (7). The amounts of O 2 consumed hydrocarbons in blocks A, B, and C and the background were calculated to be 3.9 ± 0.4 kg/a.m 2 , 3.5 ± 0.7 kg/a.m 2 , 12.3 ± 0.4 kg/a.m 2 and 0.1 ± 0.0 kg/a.m 2 , respectively. The amounts of CH 4 escaping to the air corresponding to biodegraded hydrocarbons in blocks A, B, and C and the background were calculated to be 0.9 ± 0.3 kg/a.m 2 , 0.4 ± 0.2 kg/a.m 2 , 0.3 ± 0.0 kg/a.m 2 , and 0.0 kg/a.m 2 , respectively. The calculated VOCs volatilization rates were less than 0.001 kg/a.m 2 in all areas. Therefore, the total NSZD rates were 4.7 ± 0.7 kg/a.m 2 , 3.4 ± 0.9 kg/a.m 2 , and 12.5 ± 0.4 kg/a.m 2 in accordance with the O 2 estimation.

Functional Genes Evidence
The mcrA gene, encoding the methyl coenzyme M reductase (MCR) complex, remained at a low level in the upper layer until depths deeper than 60 cm, where it remained steady ( Figure 4). This trend was the opposite of the depth-dependent oxygen variation. The correlation analysis results showed that mcrA gene concentrations had a significantly negative correlation with O 2 concentrations (r = −0.740, p < 0.001) and a significantly positive correlation with CH 4 concentrations (r = −0.686, p < 0.001), but were not significantly correlated with VOCs (r = 0.133, p = 0.518).  The concentration of the pmoA gene, encoding methane monooxygenase, did not vary as dramatically as the mcrA gene ( Figure 4). The correlation analysis results showed that the pmoA gene had a significantly negative correlation with O2 (r = −0.526, p = 0.005) and a positive correlation with CH4 (r = 0.669, p < 0.001). As Figure 5 illustrates, CH4 and the pmoA gene were not linearly correlated. The relationships can be divided into three circumstances: (I) low CH4 and low-to-middle pmoA, corresponding to the samples collected from the surface or near the surface; (II) middle CH4 and middle pmoA, correspondng to the samples collected from the intermediate depth (about 30-80 cm); and (III) high CH4 and high pmoA, corresponding to the samples collected from the deep depth (greater than 60 cm). There was no significant difference between pmoA concentrations in group I and group II (p < 0.05), indicating that the CH4 oxidation rates in these areas were basically invariant, while the pmoA concentrations in group III were significantly higher.  The concentration of the pmoA gene, encoding methane monooxygenase, did not vary as dramatically as the mcrA gene ( Figure 4). The correlation analysis results showed that the pmoA gene had a significantly negative correlation with O 2 (r = −0.526, p = 0.005) and a positive correlation with CH 4 (r = 0.669, p < 0.001). As Figure 5 illustrates, CH 4 and the pmoA gene were not linearly correlated. The relationships can be divided into three circumstances: (I) low CH 4 and low-to-middle pmoA, corresponding to the samples collected from the surface or near the surface; (II) middle CH 4 and middle pmoA, correspondng to the samples collected from the intermediate depth (about 30-80 cm); and (III) high CH 4 and high pmoA, corresponding to the samples collected from the deep depth (greater than 60 cm). There was no significant difference between pmoA concentrations in group I and group II (p < 0.05), indicating that the CH 4 oxidation rates in these areas were basically invariant, while the pmoA concentrations in group III were significantly higher.
positive correlation with CH4 concentrations (r = −0.686, p < 0.001), but were not signifi cantly correlated with VOCs (r = 0.133, p = 0.518).  The concentration of the pmoA gene, encoding methane monooxygenase, did no vary as dramatically as the mcrA gene ( Figure 4). The correlation analysis results showed that the pmoA gene had a significantly negative correlation with O2 (r = −0.526, p = 0.005 and a positive correlation with CH4 (r = 0.669, p < 0.001). As Figure 5 illustrates, CH4 and the pmoA gene were not linearly correlated. The relationships can be divided into three circumstances: (I) low CH4 and low-to-middle pmoA, corresponding to the samples col lected from the surface or near the surface; (II) middle CH4 and middle pmoA, corre spondng to the samples collected from the intermediate depth (about 30-80 cm); and (III high CH4 and high pmoA, corresponding to the samples collected from the deep depth (greater than 60 cm). There was no significant difference between pmoA concentrations in group I and group II (p < 0.05), indicating that the CH4 oxidation rates in these areas were basically invariant, while the pmoA concentrations in group III were significantly higher.

Biogeochemical Process Process in NSZD
All A-C soil gas profiles displayed more apparent O 2 usage and CO 2 and CH 4 production in the deeper subsurface. The concentrations of VOCs were about 2 to 4 or-ders lower than CO 2 and CH 4 , which suggested that nonmethane hydrocarbon contaminants were scarcely depleted through volatilization, but were rather potentially biodegraded by functional microorganisms in the source zone, consistent with previous studies [4,[13][14][15][16][17]19]. For the background (uncontaminated) samples, no VOCs and CH 4 were generated. The calculated O 2 usage and CO 2 generation due to hydrocarbon depletion fluxes were both around 0.1 kg/a.m 2 , which may merely reflect the soil respiration rates [4] and only accounted for 1%~4% of the values in the contaminated samples (i.e., A, B, and C blocks). Taken together, natural depletion was in progress in the source zone.
Below the 50cm, 40cm, and 30cm depths for blocks A, B, and C, the CH 4 concentrations were greater than 1% and almost linearly increased with depth. This result was not consistent with most of the sites, where CH 4 was always undetected in the vadose zone [4,17,19]. The high concentrations of CH 4 may have been caused by the continuous leaking of hydrocarbons in the in-service storage tanks, which supplied enough substrate for methanogenesis [4]. In anaerobic environments contaminated by LNAPLs, fermenters would biodegrade hydrocarbons and form dissolved hydrogen and/or acetate, which would be further utilized by methanogens as a substrate to form methane [26,27]. The lower CH 4 concentrations in the upper layers (<50 cm, 40 cm, and 30 cm for block A, B, and C) indicated that there might be less methanogenesis and more methane oxidation, and the O 2 concentration gradients were mainly formed by diffusion. In the deeper layers, O 2 may also be consumed by methane oxidation. Taken together, the soil gas profiles in the study area were consistent with the NSZD model showing the presence of both aerobic and anaerobic biodegradation in the source zone of the site.
The CO 2 estimated NSZD rates (2.4~4.8 kg/a.m 2 ) were slightly lower than the O 2 estimated values (2.5~12.9 kg/a.m 2 ). CO 2 is the ultimate product of hydrocarbon biodegradation; therefore, some intermediate products are excluded when estimating NSZD using CO 2 . While O 2 is used, some hydrocarbons may not be thoroughly biodegraded to CO 2 but are included in the NSZD calculation. This may result in lower estimated values of CO 2 . According to the NSZD concept model, as for O 2 -consuming processes in the vadose zone, only CH 4 oxidation is taken into consideration. Therefore, it is essential to evaluate the methanogenesis fluxes.
The calculated NSZD rates caused by methanogenesis (4~40 kg/a.m 2 ) were greater than the NSZD rates estimated by both O 2 and CO 2 gradients. One plausible reason is that the generated CH 4 was not thoroughly oxidized to CO 2 . The stoichiometric relationships between O 2 consuming (or CO 2 generating) and hydrocarbon biodegrading processes in Equation (7) are based on the reaction that CH 4 is thoroughly oxidized to CO 2 . However, in the natural environment, there are several steps to convert CH 4 to CO 2 . CH 4 may firstly be oxidized to methanol, and then formaldehyde, etc. If CH 4 was not thoroughly oxidized to CO 2 , although it was oxidized, the O 2 consumption would be low, and there would be less-or even no-CO 2 generated. This speculation is consistent with the calculated results: CH 4 generating rates > (O 2 consuming + CH 4 escaping) rates > (CO 2 generating + CH 4 escaping) rates. Although the methanogenesis fluxes were different in different blocks, the CO 2 generating flux was almost the same. Methanogenesis may not be the rate-limiting process for the mineralization of hydrocarbons due to a suite of environmental constraints [28,29], although methanogenesis has been deemed to be the rate-limiting process for NSZD [11].
Another possible reason is the overestimated gradients of CH 4 . In the estimated area, besides diffusion, there may be the presence of CH 4 oxidation and/or methanogenesis, which would vary the diffusion gradients. If the net consumed CH 4 (the oxidized CH 4 minus generated CH 4 ) in the upper soils was greater than in the deeper soils, the gradient would be overestimated and vice versa. However, it was hard to calculate the net consumed CH 4 dependent upon the monitored soil gas concentrations. Functional gene quantification may help to differentiate the CH 4 oxidation or methanogenesis processes.
The significantly negative correlation between the mcrA gene and oxygen concentrations was consistent with the consensus that anaerobic conditions are necessary for methanogenesis [30]. Although the mcrA gene was expected to be positively correlated with CH 4 , it was not significantly correlated with VOCs. Considering the low VOCs concentrations, therefore, most of the methane might not be generated by the biodegradation of VOCs in the vadose zone but by the biodegradation of petroleum hydrocarbons in the LNAPLs zone, consistent with the hypothesis of the calculation model (Equation (7)). It should be noted that MCR complexes (mcrABG subunits) can also serve as a core gene for anaerobic methane oxidation, but only in a few methanogens [31].
As the pmoA product is an O 2 -dependent enzyme [32], the negative correlation with O 2 suggested that O 2 was not a limiting factor for methane oxidation at the study depth. The positive correlation with CH 4 suggested that the amount of CH 4 provided methane oxidation activity. The pmoA concentrations in group III were higher, indicating that the CH 4 oxidation rates in the deeper soil layers were greater than those in the upper layers. The results also indicated that the CH 4 oxidation zone occurred in the deeper soils, and the CH 4 concentration gradients in the shallow depth (<60 cm) were primarily formed by diffusion, which conformed to the application condition of the NSZD assessment model. Furthermore, pmoA abundance was several orders of magnitude higher than that of the mcrA gene. This finding suggested that methane at the study depth might be mostly consumed rather than generated by functional microorganisms, and the net consumed CH 4 in upper soils might be less than those in the lower soils. This would cause an underestimation of the gradients of CH 4 . Therefore, the second possible reason for the greater methanogenesis rates is invalid, and the estimated rates of methanogenesis above were conservative. Therefore, methanogenesis was the rate-limiting process for NSZD at the study site, and the actual NSZD rates roughly equaled the methanogenesis rates, which is consistent with the previous study [11].

Environmental Implications
The calculated NSZD rates differed using different gas gradients. The different rates pointed to the different physicochemical pathways. In this study, the methanogenesis rates might equal the total NSZD rates. The NSZD rate values of 4~40 kg/a.m 2 are greater than 90% of the 40 LNAPL sites with a 90th percentile value of 4 kg/a.m 2 [9]. As previously investigated, the greatest NSZD rate was about 12 kg/a.m 2 , which was gained from accidental releases of denatured fuel-grade ethanol-affected sites [33]. The higher NSZD rates in our site might be caused by ongoing leaking from the storage tanks. The generated CH 4 was almost oxidized beneath 50cm, 40cm, and 30 cm, respectively, in blocks A, B, and C, but more than half of the CH 4 was not thoroughly converted to CO 2 . Therefore, using the O 2 consuming or CO 2 generating rate to evaluate the NSZD as by some previous studies conducted would underestimate the NSZD rate when there was excessive generated CH 4 that failed to be thoroughly oxidized to CO 2 .
On the other hand, even though the measured blocks were closely located, the NSZD rates varied. This could be caused by the unevenly distributed contaminants, which is commonly found at petroleum-contaminated sites [34]. Therefore, to gain a more accurate NSZD rate, around the source zone, measurements should be carried out as much as possible.
Multiple gases concentrations measurement combined with functional gene qualification would be a powerful means for distinguishing the different physicochemical pathways in NSZD. In the present study, only shallow soil gases and nucleic acid information were collected. Gas concentrations measured at different depths may produce different gradient values. However, the gained gradients are always in the same order of magnitudes [4]. Therefore, using the shallow subsurface gas concentration gradient to assess NSZD is theoretically acceptable. The same methods can also be applied to other petroleum-contaminated sites such as oil fields, as well as petrochemical plants. To verify the proposed deductions and gain more knowledge on NSZD, deeper samples and other molecular biological analysis should be carried out in further studies.

Conclusions
In this study, NSZD in the vadose zone at a petroleum-hydrocarbon-contaminated site of a purification plant in Northwest China was assessed using the gradient method based on a proposed minimally invasive measurement of multiple gas concentrations. Combined with the quantification results of the pmoA and mcrA genes, the concentration of CH 4 was deduced as it was generated from the deep LNAPL-contaminated zone, and the methanogenesis rates were estimated to be 4 to 40 kg/a.m 2 , which can be deemed to represent the NSZD rates. The NSZD rates were greater than those at most of the previously studied sites. The high NSZD rate at the site may be caused by ongoing leaking from the storage tanks. The generated CH 4 was almost oxidized before escaping into the air, but more than half of the CH 4 was not thoroughly converted to CO 2 . This was different during NSZD in the middle-or late-evolution stage of the LNAPL release. Therefore, our study added new data to the NSZD rate database and updated current knowledge about NSZD, especially in China.