Petroleum NAPL Depletion Estimates and Selection of Marker Constituents from Compositional Analysis

Spills and releases of hydrocarbons may result in zones of nonaqueous phase liquids (NAPL) within soils and groundwater. The NAPL will change, or “weather” over time due to a range of physical, chemical, and biological processes. Hydrocarbon constituents in environmental samples collected from NAPL‐impacted groundwater wells, sediments, or soils vary in composition over time due to weathering. The changing composition can be used to estimate mass depletion rates and trends for the bulk NAPL and for individual constituent chemicals relative to marker constituents which are less susceptible to weathering. Methods for the selection of marker constituents and for quantitatively and conservatively estimating NAPL depletion rates and trends over time are shown. Estimates are included for two sites with different NAPL mixtures present (crude oil and gasoline/diesel‐range product), for which depletion of half the initial total NAPL is estimated at 13.6 ± 2.9 years and 7.3 ± 1.8 years respectively, and with no active NAPL remediation at either site. Similar methods for oil or NAPL depletion estimates have often relied on a prior‐identified suite of presumed‐conserved marker constituents. The method presented here includes steps which identify the best set of analyzed candidate marker constituents in a NAPL mixture. This can confirm prior‐selected markers but is particularly useful for NAPL mixtures in which no prior‐identified marker constituents are present.


Introduction
PHT3D is a computer code for general reactive transport calculations, coupling MODFLOW/MT3DMS for transport and PHREEQC for chemical reactions. It was developed by Henning Prommer in the 1990s and has been applied by him and his coworkers to various groundwater problems of practical interest. The resulting publications (http://www.pht3d.org/pht3d public.html) show an impressive applicability of the code and illustrate the underlying understanding of quite complicated interactions (e.g., Prommer and Stuyfzand 2005;Prommer et al. 2008Prommer et al. , 2009). In the original version, transport is calculated during a time step, an input file is written for PHREEQC for calculating reactions such as ion exchange and precipitation or dissolution of minerals, and these steps are repeated for subsequent time steps until finished. This loose coupling has the advantage that updates of the master programs can be installed without much effort. A disadvantage is that the calculation of the chemical reactions needs to be initialized time and again for each cell in the model, which adds another time-consuming step to calculations that are already computer-intensive. Another disadvantage is that surface complexation reactions need to be calculated first using the water composition from the previous time step and then reacted with the changed water concentrations. This procedure was not implemented in the original version of PHT3D, and surface complexation reactions could not be calculated.
Prommer and Post recently released the second version of PHT3D that resolves the shortcomings and works very well. The improvement is owing firstly to the implementation of total-variation-diminishing (TVD) scheme that MT3DMS uses for calculating advective and dispersive transport (Zheng and Wang 1999). Secondly, it is because PHREEQC is now being used for storing the chemical data of the model, including the chemical activities and the composition of surface complexes from the previous time step. In addition, the procedure to transport total oxygen and hydrogen has been adapted from PHAST (PHAST is the 3D reactive transport model developed by Parkhurst et al. 2004, based on HST3D and PHREEQC). This enables the user to obtain the redox state of the solution without having to transport individual redox concentrations of the elements (e.g., C being distributed over carbon-dioxide, C(4), and methane, C(-4)). The tighter coupling quickens the calculations twofold at least, but probably by an order of magnitude for the more interesting cases. In this review, the background of the new implementation is presented and illustrated with examples and compared with results from PHREEQC and PHAST.

Introduction
Spills and releases of hydrocarbons may result in zones of nonaqueous phase liquids (NAPL) in soils and groundwater (Tomlinson et al. 2017). There are significant challenges in assessing risk and remediation effectiveness for NAPL zones, including estimating the depletion of NAPL mass over time.
Direct evaluation of longevity may include measuring and interpolating the changes in total in situ NAPL mass or volume over time. The longevity of a subsurface hydrocarbon NAPL zone may also be evaluated indirectly through measurement of constituent volatilization and dissolution flux from the NAPL, measurement of biodegradation reactants (O 2 , other electron acceptors) or products (CH 4 , CO 2 ) entering or leaving the NAPL zone, measurement of conserved or near-conserved gases (Ar, N 2 ) near NAPL, or with estimates of generated heat flux from the NAPL zone during degrada-tion (Amos et al. 2005;Johnson et al. 2006;Lundegard and Johnson 2006;API 2017;Garg et al. 2017;CRC CARE 2018).
Another option for direct estimates of NAPL depletion is based on the change in mixture composition over time. Examples include forensic applications in which a selected marker constituent is presumed to be conserved. Depletion of the remaining constituents is evaluated relative to the marker constituent (Douglas et al. 1996;Arey et al. 2007;Yang et al. 2015). An ideal NAPL marker constituent would normally be relatively nonvolatile, insoluble in water, and resistant to chemical and biological degradation. Typically these marker chemicals, if present, may include hopanes, steranes, naphthenes, naphthenobenzenes, or highly branched alkanes. Persistence of a constituent in NAPL, however, may also depend on the relative significance of these depletion mechanisms (volatilization, dissolution, degradation) under site-specific conditions; a NAPL constituent may be a viable conservative marker at some sites and under some conditions, but not at other sites or conditions.
A quantitative analysis method is presented which includes mathematical relations between mass fractions in NAPL and chemical concentrations in soils or sediments, a discussion of application, and example depletion estimates. The method is consistent with constituent marker methods (Douglas et al. 1996), but also specifically includes steps for selecting the most appropriate conserved marker constituents based on measured composition data, parameter transformations, and regression analyses. As presented the method Additional Supporting Information may be found in the online version of this article. Supporting Information is generally not peer reviewed. requires a finite NAPL mass with no added hydrocarbons from additional spills or releases during the monitored period; site history and forensic evaluation methods (Wang et al. 2006;ASTM D3328 2013) will aid in this determination.

Methods
NAPL composition measurements (samples taken at specific times, from one or more sampling points) are of constituent mass fraction in NAPL, χ i (g-i/g-oil), and χ T = 1 (g-oil/g-oil) is the total sum of the N constituent mass fractions. For a colocated soil volume, V, containing the NAPL, the constituent mass concentrations, c i (g-i/g-soil) sum to the total NAPL concentration in soil, c T (g-oil/g-soil). Mass in soil is proportional to a volume-averaged concentration M i = V ⋅ ρ s ⋅ c i . The associated soil averaging volume, V, and soil bulk density, ρ s , are presumed constant over time and therefore NAPL constituent soil concentration, c i , is proportional to NAPL constituent mass, M i . The constituent mass concentrations in soil and the mass fractions of NAPL are related, For a finite NAPL release, we specify a constituent concentration trend c i (t) = c T (t) ⋅ χ i (t) as a continuous constant or decreasing monotonic function [c i (t) ≥ 0, dc i (t)/dt ≤ 0 for all t ≥ 0] and the time rate of change dc i (t)/dt varies between constituents. Rates for each constituent are found.
For a binary mixture with i = A as a presumed conserved constituent (or summed conserved constituents A, and B is the sum of the remaining constituents), ∂ln[c A (t)]/∂t = 0, and, from Equation 2: For conserved mass fraction A at times, t 0 and t 1 , with Equation 3, and for c T (t 0 ) = 1: Application of Equation 4 is illustrated in Figure 1 for specified values of conserved mass fractions, A, in a binary mixture with χ A + χ B = 1.
In a general case, we have a mixture of many NAPL constituents. The most conserved constituents in a mixture are identified as the maximum (+) values of ∂ln[χ i (t)]/∂t in a set of N constituents. If there are no known conserved marker constituents, and the A constituent (or summed constituents) actually deplete over time, ∂ln[c A (t)]/∂t ≤ 0, and the estimated total depletion rate is an underestimate of the actual depletion rate. The ordered ranking for the set of N constituents by ∂ln[χ i (t)]/∂t from highest (+) to lowest (−), is the same as the ordered ranking of constituent depletion rates, ∂ln[c i (t)]/∂t.

Application
Here we develop a regression method applicable for a general data set. NAPL depletion estimates are illustrated for a data set of measured mass fractions, χ i (t j ), for i = 1 to N constituents at j = 1 to M discrete sampling times. The decrease in total NAPL concentration, c T (t), is first estimated by evaluating each individual constituent i in the NAPL mixture as a potential conserved marker in a binary mixture.
For a relative depletion estimate we specify an arbitrary initial total concentration c 0, T(i) = c 0, A(i) + c 0, B(i) = 1. Concentration trends for A and B, are represented by c A(i) These functions meet the requirements for a continuous constant or decreasing monotonic function [c i (t) ≥ 0, dc i (t)/dt ≤ 0 for all t ≥ 0] but are not unique. The mass fraction of constituent A is: The best-fit coefficients φ 0, A(i) and κ A(i) in Equation 7 are determined for M data pairs in time, [t j , ln φ A(i) (t j )], with linear regression (Bendat and Piersol 2000) for each of i constituents. Mass fraction trend from Equation 7 is: The natural log rate of change for Equation 8 is: We rank order the constituents using The introduced error in the approximation is proportionate to −χ 0, A(i) . For individual candidate marker constituents with positive κ A(i) in the examples (to follow) all χ 0, A(i) < 0.10 and the approximation is biased low by, at most, −10%.
The maximum value of κ A(i) in the N constituents is identified as the most conserved constituent in a NAPL mixture. Other potentially conserved or depleting constituents are identified through the ordered ranking of constituents by the mean rates, κ A(i) , from high to low. This ordered ranking is identified by the i-th constituent index. Discussion of confidence limits is included in the Supporting Information.
Multiple potentially conserved marker constituents may be identified. Summing the marker constituents may improve the depletion estimate. In a second step we define a set of binary mixtures where χ A(q) (t) includes the sum of the q = 1 to i-th individual mass fractions ordered with the greatest rate of mass fraction increase.
These summed mass fractions χ A(q) are evaluated for trend similar to the individual mass fractions, χ A( i) , as in Equation 7. This yields, for each of the q = 1 to i-th summed constituents the best-fit coefficients, φ 0, A(q) and κ A(q) . For a potentially conserved marker constituent or summed constituents (q = 1 or q = 1 to i-th), the evaluation yields a summed rate, κ A(q) , with total NAPL concentration The natural logarithm rate of change for Equation 11 is The best estimate of total mean depletion for Equation 11 is the minimum value of c T(q) (t) for any q at a specified time, t, including an acceptable confidence limit. This best estimate may result from different q at varied times.
In the examples to follow we chose one optimum marker sum q for all times. Additional possible steps in evaluating and selecting marker constituents are included in the Supporting Information.
The trend in individual constituent concentrations c i (t) = c T (t) ⋅ χ i (t) is estimated as a product of estimated total concentration, Equation 11 and constituent mass fraction, If a sum of constituents (q > 1) is selected as the conserved marker, the included individual marker constituent concentrations are expected to be relatively constant over time.
The best estimate of constituent mean depletion for Equation 13 at a specified elapsed time, t, is the minimum value of c i (t)/χ 0, A(i) , including an acceptable confidence limit. The optimum sum q may vary over time, and for different constituents, i. In the examples we chose one optimum sum of marker constituents q for all estimated times, t, and for all constituents, i. Equations 11 and 13 for total and constituent depletion are consistent with estimates which presume the marker constituent is known (Douglas et al. 1996, Yang et al. 2015. With Equation 8 into Equation 11, for measured data values: Similarly, with Equation 14 into Equation 13 The ratio of measured mass fractions in Equation 15 is directly comparable to the constituent ratios often applied in NAPL forensics evaluations (Arey et al. 2007).

Results
The method has been applied to NAPL compositional data measured at two sites with different NAPL mixtures present (crude oil and gasoline/diesel-range product). The intent in these examples is demonstrated selection of viable marker constituents, estimation of trend in total and constituent NAPL depletion over time, and comparison of observed similarities and differences in the results. Other applications of the method are possible, and would vary depending on defined objectives, which may include: • How much NAPL is remaining now compared to the initial release? • How long for half (or most) of the remaining NAPL to deplete? • Are constituents of potential concern in the NAPL depleting? • How long for half (or most of) each constituent of potential concern to deplete?
Data requirements would, of course, vary depending on the stated objectives. The methods for sample preparation and chemical analysis applied for the example data sets are documented (Yang et al. 2015) and widely available. Further discussion on chemical analysis is included in the Supporting Information.

Example 1. Bemidji, Minnesota, USA Crude Oil Release Site
A crude oil release from a pipeline to land and the shallow subsurface at Bemidji, Minnesota, is described elsewhere (Essaid et al. 2011). NAPL data collection and analysis is reported and applied in estimates of total and constituent mass depletion using varied data sets (Lundy 2014;Baedecker et al. 2018). We evaluate marker selection and depletion estimates using the data reported by Lundy (2014). Lundy (2014) presents data and depletion estimates for the north pool area of the Bemidji site using pristane as the selected marker. Mass losses from weathering were estimated based on ratios of each constituent relative to pristane in the initial reference sample and constituent depletion was summed to estimate total mass loss over time. Trend was represented by a first-order rate of 0.0308/year, corresponding to 36.4% of the initial total oil mass remaining after a 32.8 year time period. The estimate is generally consistent with other measurement and estimation methods presented in Lundy (2014).
In evaluation of the Lundy data with the present method, we have N = 56 detected and quantified hydrocarbon constituents in the range of C4 to C40. Constituents present at nondetect levels, including low intermittent detection of light olefins were omitted. The summed olefin mass fraction ranges from 0 to 0.0004; the effect of omitting and renormalizing was minimal. Data values reported as mole fractions relative to pristane were renormalized to mass fractions in oil using constituent molecular weights. The data includes M = 39 samples collected over 32.8 years.
Regression with Equation 7 yields coefficients χ 0, A(i) and κ A(i) for each of the N constituents. All χ 0, A(i) values are small and Equation 9 can be simplified; the constituents have been rank ordered on κ A(i) from high (+) to low (−). Results are shown in Figure 2a and 2b. The top ranked chemical is pristane, same as identified by Lundy (2014). Trends for the top 15 ranked individual constituent rates are all increasing within the specified 95% two-sided confidence interval. The top 21 ranked individual constituents have positive mean κ A(i) and could be potential individual marker constituents.
Confidence intervals in the parameter estimates and the related depletion estimates in some cases may be improved by summing multiple marker constituents. The time series of the ranked q = 1 to i-th individual mass fractions are summed as in Equation 10 with regression following Equation 7. The resulting coefficients χ 0, A(q) and κ A(q) are also included in Figure 2a and 2b. Both the ranked summed mass fraction coefficients χ 0, A(q) and the rate parameters κ A(q) increase monotonically for greater q. For the summed constituents the mass fraction χ 0, A(q) is not necessarily small and Equation (9 cannot be generally simplified; the parameter κ A(q) increases for larger q but applies for a smaller fraction of presumed depleting constituents, (1 − χ 0, A(q) ).
The optimum selection of summed marker constituents is based on estimates of remaining total concentration (or fraction of initial total mass) at a final sampled time of t = 32.8 years using Equation 11 and the fitted χ 0, A(q) and κ A(q) parameters. This is shown in Figure 3a versus the rank order of κ A(i) for individual constituents. The minimum estimated final concentration, including the upper confidence interval, is relatively insensitive to the number of ranked summed marker constituents for q = 1 to 22. A plot of the remaining total concentration (or mass) trend for q = 12 is shown in Figure 3b, where the trend line is for Equation 11 and the data points are from Equation 14. The nonlinear trend in Figure 3b results from the fitted regression coefficients, the long data record duration, and the applied regression equation, Equation 11. Further possible steps in marker constituent selection are included in the Supporting Information.
An effective first-order rate k T, eff, q (t) = (∂[lnc T(q) (t)]/∂t) from Equation 12 is a function of time; the corresponding instantaneous half-life [t 1/2, eff = − ln(0.5)/k T, eff, q (t)] is also a function of time. An effective first-order total concentration rate at t = 0 from Equation 12 is: For q = 12, k T, eff, q (t = 0) = −0.056 ± 0.01/year, with an equivalent initial half-life of 12.3 (10.4 to 15.1) years. This half-life at t = 0 differs from the estimated time to reach c T(q) = 0.50 using Equation 11 of 13.6 ± 2.9 years.
Constituent-specific results are presented both in terms of final concentration at t = 32.8 years and natural logarithm rate at t = 0 in Figure 4. An effective constituent-specific rate at t = 0 from Equations 2, 9, and 12 is Constituent depletion is related to an effective half-life, t 1/2, eff = − ln(0.5)/k eff . Constituent depletion relative to total depletion at t = 0 is given by the ratio Relative constituent-specific depletion rates are shown in Figure 4a using Equation 18. Ratios of the time for constituent to time for total concentration to reach 50% of the initial values, using Equations 11 and 13, are also included. Toluene (initial half-life: 3.3 years) decreases at the fastest initial rate, or approximately 3.7 times the total mean depletion rate. The estimated remaining concentration (or fraction) for each constituent at the final sampled time of 32.8 years are shown in Figure 4b. The estimated mean total concentration remaining is 0.24; of the 57 constituents, 19 have a mean estimate of less than 0.05 remaining at 32.8 years. These values may be underestimates if the selected marker chemicals are also depleting. Nothing internal to the method (or any presumed conserved marker method) will identify if a marker itself is depleting.

Example 2. Petroleum Terminal Site, USA
Operations at a petroleum fuel products terminal site in the mid-continent USA have led to release of mixed refined oil products to land and the shallow subsurface. The released fuel types and composition vary across the site. We estimate mass depletion from NAPL compositional changes over time with data for a single monitoring well impacted by a gasoline-diesel product mixture. No additional spills occurred in the well vicinity during the measured time period. While the well was impacted by multiple fuel products the chromatogram fingerprints for this well remained substantially similar over the monitored time period. No sample of the unweathered released NAPL mixture was available.
Samples of NAPL were collected using a low flow peristaltic pump, stored and shipped in 40 mL VOA vials, and analyzed in our laboratory using several approaches. Direct injection gas chromatography with flame ionization detection (GC/FID) was used to characterize the product type and the mass distribution of selected diagnostic components (ASTM Method D3328-06 2013; USEPA Method 8015 Modified, USEPA 2003). Direct injection gas chromatography with mass spectrometry (GC/MS) was also used for the detection and quantitation of specific hydrocarbon components (USEPA Method 8260, USEPA 2018).
The evaluated well included M = 18 quarterly sampling events over 4.1 years and N = 17 initially identified and quantified hydrocarbon constituents ranging from C5 to C20. The NAPL samples contained both gasoline (approximately 28% < n-C10) and diesel-range material. The quantified constituents account for approximately 29% of the C5 to C20 range hydrocarbons, with the remaining mass within a "hump" of unresolved complex mixture (UCM) in the diesel range (≥ n-C10). A small fraction (<5%) of unresolved mixture from C20 to approximately C23 is evident in the analysis. There is no evidence to suggest that higher molecular weight hydrocarbons (>C23) are present in the NAPL samples. The identified and quantified constituents were normalized to a summed mass fraction of unity. The total depletion estimates apply for the quantified and summed constituents included in the analysis.
The evaluation of individual and summed constituents followed the same methodology as the prior example. In this case, isopentane (i-C5) was identified with the maximum increasing mean rate, κ A(i) = 0.21 ± 0.24 /year, but with a decreasing rate (−0.03/year) at the lower confidence limit. Because of the uncertainty in the estimated i-C5 increasing rate, this constituent was omitted from the mass fraction distribution and the remaining mass fractions were renormalized to unity. The mean mass fraction for i-C5 in the samples was 0.0084; the effect of omitting i-C5 and renormalizing the mass fractions was minimal.
Results for the revised data are included in Figure 5a and 5b for individual and summed coefficients and rates. The minimum estimate of remaining mass fraction, for both the mean and upper confidence limit, at the final sampled time of t = 4.1 years occurs for q = 1 as shown in Figure 6a. The confidence limit range on remaining mass fraction significantly overlaps for any of q = 1 to 7 (approximately). We select q = 4 including n-C10, n-C13, n-C5, and n-C6 as marker constituents, all of which are normal paraffins. A plot of the remaining total concentration (or mass) trend for q = 4 is shown in Figure 6b. We estimate the time to reach c T(q) = 0.50 at 7.3 ± 1.8 years. The linear trend in Figure 6b is a result of the small exponential argument in Equation 11 at the maximum time duration, with (κ A(q) ⋅ t) < 1.
Constituent-specific results are presented both in terms of natural logarithm rate and remaining concentration. Relative constituent-specific depletion rates are shown in Figure 7a, along with the ratios of time for constituent and time for total concentration to reach 50% of the initial values. The normal alkane n-C19 (initial half-life: 3.2 years) is estimated to decrease at the fastest initial rate, or approximately two times the total rate (with an initial half-life of 6.6 years). Remaining mean concentration for 18 of the 29 constituents is less than that of the remaining mean total concentration of 0.66 over the 4.1 year sampled time period, as in Figure 7b. These values may be an underestimate if the selected marker chemicals are also depleting. Nothing internal to the method (or any presumed conserved marker method) will identify if a marker itself is depleting.

Discussion
Methods and results for estimating total and constituent depletion trends based on compositional change are presented for two case examples. The time for depletion of half the initial total NAPL is estimated at 13.6 ± 2.9 years and 7.3 ± 1.8 years in the example cases, with respectively, crude oil (C4 to C40), and a gasoline/diesel mixture (C5 to C20) present, and no active NAPL remediation at either site. The illustrated methods may be useful in assessing similar petroleum impacted sites. The methods may also be useful in comparison or in conjunction with evaluation of active site remediation technologies.
In Example 1, for crude oil (C4 to C40), the selected summed marker chemicals have been identified predominately as branched alkanes, consistent with chemicals most often identified as persistent biomarkers (Wang et al. 2006), including pristane, same as identified by Lundy (2014). In Example 2, for a gasoline/diesel mixture (C5 to C20), some conventional biomarkers (of >C20 range) are not present. In this case, the method-identified markers are normal paraffins, C14 and lighter. These would not normally be identified as persistent. This range of normal paraffins (<C14) has, however, been observed to be more persistent than other higher molecular weight normal paraffins in another diesel release (Hostettler et al. 2013). In both included examples, the marker selection method only indicates these markers are more persistent than other analyzed constituents in the NAPL; the markers may also be depleting.
In both examples the constituent-specific depletion rates for the C12 to C20 range n-paraffins exceed the rates for many branched alkanes in the same C12 to C20 range, consistent with prior observations (Wang et al. 2006;Wentzel et al. 2007). The C12 to C20 n-paraffin depletion rates also exceed the rates for some more volatile and lower molecular weight n-paraffins; this has been observed previously for both the Bemidji site and a diesel product release site  (Hostettler et al. 2013). Indirect evaluation of NAPL depletion solely through measurement of water-soluble constituents in groundwater or volatile hydrocarbon constituents in soil gas would not capture the evident depletion of these C12 to C20 range n-paraffin constituents in NAPL.
We also expect that the observed depletion of the C12 to C20 n-paraffin range NAPL will also ultimately contribute to a portion of the methane and carbon dioxide produced in petroleum NAPL biodegradation. Biodegradation of petroleum NAPL in the example cases likely occurs under methanogenic geochemical conditions. Presuming no biomass growth and no persistent intermediate metabolites, complete degradation of petroleum in methanogenic conditions optimally produces methane and carbon dioxide (Shelton and Tiedje 1984 Biodegradation of n-paraffins (>C10 up to C50) by mixed microbial cultures in methanogenic conditions has been observed for oil-impacted sediments and oil reservoir fluids in microcosm studies (Zengler et al. 1999;Anderson and Lovley 2000;Jones et al. 2008;Siddique et al. 2011;Wawrik et al. 2016;Liu et al. 2020) with measured production of methane in the range of 75% to 94% of that predicted by Equation 19. In the example cases and in similar NAPL releases, we may also expect production of methane, with oxidation of some produced methane to carbon dioxide in aerobic soil zones surrounding the NAPL.
In addition to petroleum NAPL, the evaluation method may also be applied for measured petroleum composition in soil or sediment samples, as in Douglas et al. (1996), and may be useful for differentiating total-and constituentdecreases due to mixing and dilution from decreases due to degradation and depletion.
The duration and extent of groundwater and vapor impacts for volatile and water-soluble chemicals at petroleum release sites depends on the presence and concentration of those chemicals in remaining NAPL. Evaluation of constituent-specific trends in NAPL mass fraction (and mole fraction, with Raoult's law), may aid in evaluating and reconciling the corresponding measured trends in adjacent groundwater and soil vapor.
Other methods for NAPL depletion estimates are available. The presented method only estimates relative rates for total and constituent depletion; a complementary assessment which provides an estimate of total mass at some instant in time could be used to tie the relative estimate to an absolute mass depletion estimate, as in Lundy (2014).
The included examples apply linear regression; robust statistical methods could be applied and may be appropriate, depending on observed data distributions, outliers, and censoring due to quantitation limits.
Refinement of the examples could include evaluation of spatial trends. In the first example trend is evaluated for the entire data set and yields an average depletion estimate over time for the entire data set. Separation in the data points ( Figure 3b) is noted at later elapsed times. Evaluation of subsets of the data, as in Baedecker et al. (2018), could be applied to evaluate spatial dependence of trends across the site. In the second example, different petroleum products are present in different wells and use of the method is shown for one well. The evaluation method can be repeated wellby-well, along with forensic evaluation methods, such as in ASTM D3328, to ensure reasonable consistency of NAPL composition in single wells, in compiling a set of well-specific depletion estimates across the site. The consolidated set may illuminate differences (or similarities) in depletion across the site in different wells, and for different product compositions.
The authors view the presented methodology as more of an outline than a cookbook. A key idea is in better identification of the "most conserved" (if not totally conserved) marker constituents with compositional analysis and evaluation of constituent changes, rather than basing marker selection on prior observations for other petroleum compositions at different sites in different physical conditions. The methodology, for extensive data sets, could be applied in evaluation of trend over time, as shown in the example data sets. For less extensive data, the method can be simplified to an evaluation of as few as two data points (or groups of data at two times), as discussed in the Supporting Information. Confidence limits on the estimates would improve with more data.