Dense nonaqueous phase liquid (DNAPL) source zones comprise persistent sources of groundwater contamination that are recalcitrant to complete remediation using conventional (e.g., pump and treat) or emerging (e.g., surfactant flushing) technologies. Increased attention to the assessment of the benefits of partial mass removal from such contaminant source zones has intensified efforts to model multiphase flow and transport behavior. This paper describes the simulated recovery of a tetrachloroethene (PCE) spill in a statistically homogeneous but nonuniform aquifer, incorporating nonuniformity in both nonaqueous phase liquid saturation and pore velocities. We developed a ganglia-to-pool metric to quantify DNAPL source-zone architecture, and explored the correlation of this metric with dissolved mass flux behavior in response to partial DNAPL mass removal. Dissolution of 20%–70% of PCE mass from models exhibiting low ganglia-to-pool ratios resulted in a larger predicted reduction of dissolved contaminant mass flux than models with high ganglia-to-pool ratios. Results of this study suggest that DNAPL source-zone characterization at field sites with homogeneous, nonuniform aquifers would benefit from inclusion of an estimate of the overall ganglia-to-pool ratio. Simulations demonstrate that flux reduction behavior depends on the source-zone architecture, which is not readily predictable using a priori assumptions about the spatial correlation of physical aquifer parameters. Model results further suggest that stochastic investigations of DNAPL source remediation at field sites should avoid reliance upon Leverett scaling of capillary entry pressures to permeability fields, which can artificially narrow the range of simulated behaviors.

Following the accidental release of a dense nonaqueous phase liquid, spatial variability of physical and chemical soil properties can exert a controlling influence on infiltration pathways and the distribution of entrapped organic liquids (Kueper et al., 1993; Lemke and Abriola, 2003; Phelan et al., 2004). Organic liquid spreading, fingering, and pooling typically result in dense nonaqueous phase liquid (DNAPL) source zones characterized by contaminated regions of irregular organic liquid saturation with complex boundaries (Anderson et al., 1992; Feenstra et al., 1996; Kueper and Frind, 1991). Spatial variability in aquifer properties further influences groundwaterflow pathways and pore velocities that affect non-aqueous phase liquid dissolution dynamics and the transport of contaminants dissolved in the aqueous phase (Miller et al., 1998).

The inherent complexity of DNAPL source zones contributes to their recalcitrance to natural attenuation and engineered remediation efforts. Field, laboratory, and modeling efforts have demonstrated the difficulty in achieving complete removal of entrapped DNAPLs from nonuniform subsurface formations (EPA, 2003; National Research Council, 1994, 2004; SERDP/ESTCP, 2002; Stroo et al., 2003). Consequently, investigations have begun to explore the benefits of partial source-zone remediation, such as the reduction of DNAPL mobility, reduced source longevity, and the enhanced efficiency and effectiveness of complementary remediation technologies (e.g., Christ et al., 2005; Rao et al., 2002; Saenton and Illangasekare, 2003). The need to quantify the potential benefits of partial mass removal from DNAPL source zones has intensified efforts to model multiphase flow and transport behavior.

The inherent complexity of DNAPL source zones also makes their behavior challenging to model. Recent efforts have relied on both analytical and computational models. Numerical models assuming various idealized conditions, such as nonuniform DNAPL distributions within uniform flow fields (Sale and McWhorter, 2001) or uniform DNAPL distributions within nonuniform flow fields (Rao and Jawitz, 2003), have yielded conflicting predictions of contaminant mass flux reduction in response to partial DNAPL mass recovery. Simplified models assuming positive and negative correlations between nonaqueous phase liquid saturation and permeability have also been proposed (e.g., Berglund, 1997; Rao et al., 2002). This paper describes the simulation of a tetrachloroethene (PCE) spill in a statistically homogeneous but nonuniform aquifer, incorporating nonuniformity in both nonaqueous phase liquid saturation and pore velocities. We analyzed predicted dissolved phase mass flux as a function of PCE recovery from simulated DNAPL source zones for realizations constructed using different assumptions concerning the spatial correlation of aquifer properties. We quantified DNAPL source-zone architecture in terms of organic liquid mass distributed in pools and zones of residual saturation, or ganglia, and explored the correlation of this architecture metric to dissolution mass flux.

An unconsolidated, unconfined aquifer located at the site of a former dry cleaning business in Oscoda, Michigan, United States, provided the basis for alternative aquifer characterization models used in this study. The aquifer is composed of relatively homogeneous glacial outwash sands and is underlain by an extensive clay layer ∼8 m below the ground surface. A PCE DNAPL source zone beneath a building at the site was the target of a pilot-scale surfactant-enhanced aquifer remediation (SEAR) test conducted in the summer of 2000 (Abriola et al., 2005). The spatial distribution of aquifer properties and the magnitude of hydraulic gradients utilized in this investigation are representative of, but not meant to explicitly reproduce, the conditions found at the field site.

Aquifer characterization data collected from continuous core samples at the site are described by Lemke (2003) and Lemke et al. (2004a). The mean porosity, ϕ, was 0.36. Nonuniform permeability, k, values estimated from 167 grain-size distributions ranged from 1.3 × 10−8 to 6.5 × 10−7 cm2. Variance of the log-normal transformation of the estimated k values (σ2 ln[k]) was 0.29. Capillary pressure–saturation relationships were modeled using the Brooks and Corey (1964) formulation. Values for the capillary entry pressure, Pb, and the pore size distribution factor, λ, were measured in the laboratory for representative aquifer samples and estimated directly from grain-size distribution curves and porosity values using the Haverkamp and Parlange (1986) method (Lemke et al., 2004a).

The influence of alternative treatment of the correlation between spatially variable porosity, permeability, and capillary retention parameters on PCE entrapment and dissolution was evaluated by Lemke et al. (2004a, 2004b), who found that the method of scaling of capillary pressure–saturation parameters to the permeability fields created the most significant differences in ensemble model results. Capillary pressure–saturation parameters were either scaled to the ϕ and k fields using Leverett (1941) scaling, or were estimated independently of the k field using the Haverkamp and Parlange (1986) method. For ensembles that did not employ Leverett scaling, PCE infiltration from a hypothetical spill exhibited increased depth of penetration, decreased spreading, and an increased tendency for DNAPL pooling, reflected in maximum organic saturations exceeding 0.60, even when compared to Leverett-scaled ensembles with identical ϕ and k distributions (Lemke et al., 2004a). Subsequent modeling of PCE dissolution and recovery under natural and engineered conditions predicted one- to two-order reductions in dissolved phase mass flux that were linked to the correlation of physical aquifer properties and associated with the evolution of the DNAPL source-zone architecture (Lemke et al., 2004b).

This paper extends the analysis of Lemke et al. (2004a, 2004b) to investigate DNAPL persistence and dissolved phase mass flux in relation to DNAPL source-zone architecture quantified using a new ganglia-to-pool (GTP) ratio metric. Lemke et al. (2004b) adopted a working definition for pools as model cells with organic liquid saturations exceeding the maximum residual organic saturation, Sor, and ganglia as cells containing organic liquid saturations at or below Sor. This definition is applied here to compare ensembles and individual realizations on the basis of ganglia-to-pool ratios calculated using the PCE mass fraction in cells with saturations above and below the modeled Sor value.

Two aquifer model ensembles are considered for this analysis. The ensembles shared identical ϕ and k fields but employed different methods of assigning capillary pressure–saturation parameters. The first ensemble employed Leverett scaling. Consequently, aquifer realizations in this ensemble contain perfectly correlated Pb and k fields. In contrast, the second ensemble utilized the Haverkamp and Parlange approach to estimate λ and Pb values without consideration of the k field. As a result, λ and values in the second ensemble Pb are strongly associated with the spatial distribution of independently simulated geostatistical indicator classes corresponding to sediment grain-size distributions (Lemke et al., 2004a). Based upon this difference, realizations in the ensemble employing Leverett scaling are herein referred to collectively as correlated realizations, while those from the ensemble generated using the Haverkamp and Parlange method are referred to as independent realizations.

PCE infiltration and entrapment simulations for 50 realizations from each of the two ensembles were used as the basis of this work. These simulations were produced using Michigan Vertical and Lateral Organic Redistribution (MVALOR) (Abriola et al., 1992; Rath-felder and Abriola, 1998), a two-dimensional immiscible fluid-flow model. Some simulation results and a complete description of the PCE spill modeling procedure are given in Lemke et al. (2004a). In the present study, 50 new independent MVALOR simulations were completed to correct for a recently discovered coding error reported by Lemke et al. (2005). DNAPL entrapment behavior predicted by MVALOR was quantified for each ensemble using PCE saturations, the ganglia-to-pool ratio, and vertical and horizontal spatial moments.

DNAPL persistence and dissolved phase mass flux were evaluated using MVALOR PCE distributions as initial conditions. Rate-limited PCE dissolution was simulated using MISER (Michigan Subsurface Environmental Remediation) (Abriola et al., 1997; Rathfelder et al., 2000), a two-dimensional compositional flow and transport model. Details of the MISER modeling approach are provided in Lemke et al. (2004b). Two sets of dissolution simulations are considered here. In the first, a 2 d surfactant flush, followed by a 10 d water flood was simulated for each of the 50 realizations in both ensembles. In the second, the duration of the surfactant flush was varied between 0.001 and 25 d over a series of runs to predict differences in downgradient dissolved phase mass flux (i.e., aqueous phase contaminant concentration) as a function of PCE mass removal from the DNAPL source zone. Modeled remediation efficiencies in these simulations range from 0% to 99% PCE recovery.

Because of excessive computational demands for the second set of simulations, a subset of realizations representing the range of predicted source-zone configurations was selected for analysis from each ensemble 01. Realizations from the correlated ensemble were chosen to encompass minimum, mean, and maximum spreading, infiltration, and saturation metrics for the initial DNAPL source-zone condition. Initial ganglia-to-pool ratios for these realizations deviated by less than one standard deviation from the mean ganglia-to-pool ratio for the correlated ensemble. Comparative realizations from the independent ensemble were selected on the basis of the observed range of ganglia-to-pool ratios defined using an Sor value of 0.151, the default in the MVALOR simulator. These realizations spanned the range of initial ganglia-to-pool ratios (1–66) observed in the independent ensemble.

MVALOR PCE spill simulations revealed considerable differences in predicted source-zone architecture that were related to the correlation of physical aquifer properties. PCE saturation statistics for the correlated and independent ensembles are given in 02 and 03, respectively. Typical distributions of pool and ganglia cells are illustrated in Figure 1, and PCE saturation distributions for the first realization in each ensemble are plotted in Figure 2. Greater variability in ganglia-to-pool ratios was found for the independent ensemble compared to the correlated ensemble (Fig. 1). Predicted ganglia-to-pool ratios ranged from 0.9 to 65.6 for the independent ensemble, extending above and below the narrower range of 2.1–14.3 observed for the correlated ensemble. PCE saturation metrics reported in 03 and depicted on Figure 2 also indicate increased depth of penetration, decreased spreading, and higher organic saturations for the independent ensemble.

Calculated ganglia-to-pool ratios for independent ensemble realizations were found to be insensitive to the selected Sor values within a reasonable range (between 0.10 and 0.20) because most of the simulated PCE mass was distributed in model cells with very high or very low saturations. Correlated realizations, which displayed a more continuous saturation distribution in model cells, were comparatively more sensitive to the selected Sor cutoff value. However, the variation of calculated ganglia-to-pool ratios for each realization was small compared to the total range of ganglia-to-pool estimates among realizations (GTP = 0.91–87.9). Furthermore, none of the realizations examined in this study would have been re-classified from a high to low ganglia-to-pool ratio (defined based on the criteria discussed below) or vice versa across the stated range of Sor values.

MISER simulations revealed differences in mass recovery and effluent concentration predictions that were linked to variations in predicted source-zone architecture. Figures 3 and 4 show a comparison of correlated and independent model responses to a simulated 2 d surfactant flush for the first realization in each ensemble. In the correlated realization (Fig. 3), removal of 97.8% of the initial PCE mass distributed throughout the DNAPL source zone was accompanied by a decrease in dissolved phase PCE mass flux from 44.6 to 2.1 g/d under natural gradient flow conditions. This corresponded to a decrease in the flux-weighted PCE concentration from 122 to 5.8 ppm across the downgradient model boundary. In the complementary independent realization (Fig. 4), removal of 73.5% of the initial PCE mass was accompanied by a decrease in dissolved phase PCE mass flux from 32.1 to 1.0 g/d, with a corresponding decrease in the flux-weighted PCE concentration from 88.8 to 2.7 ppm.

The behavior of the realizations depicted in Figures 3 and 4 is representative of model behavior for the correlated and independent ensembles. In general, correlated and independent ensembles generated similar PCE mass flux and effluent PCE concentrations, both prior to and following simulated 2 d surfactant flushing 04. Slightly greater PCE mass removal was predicted from correlated realizations (96.2% on average) compared to independent realizations (91.2% on average). This degree of mass removal was accompanied by a one- to two-order decrease in predicted PCE effluent mass flux for both ensembles. For the example shown in Figure 3, simulated surfactant treatment of the correlated realization removed PCE entrapped at residual saturations from most of the DNAPL source-zone cells. The remaining PCE, distributed in several cells throughout the model (Fig. 3B), functioned as a continuing source of dissolved phase contaminant (Fig. 3C). Similarly, residual PCE was removed from the independent realization during surfactant flushing, leaving behind zones of pooled PCE (Fig. 4B) that continued to dissolve (Fig. 4C).

Simulation of downgradient dissolved phase mass flux in response to varying degrees of PCE mass removal from the DNAPL source zone suggests that flux versus percent removal curves can be segregated on the basis of the initial DNAPL ganglia-to-pool ratio. Flux-weighted effluent PCE concentrations for six representative correlated realizations 01 are plotted as a function of PCE removal in Figure 5. The six curves plotted in Figure 5 are clustered, demonstrating a consistent response to PCE mass removal among the correlated realizations.

Flux-weighted effluent PCE concentrations for seven independent realizations 01 are plotted as a function of PCE removal in Figure 6. The plotted curves encompass the range of variability defined by the six correlated realizations. Realizations with low ganglia-to-pool ratios (<2) are associated with an accelerated decrease in PCE mass flux in response to PCE recovery, whereas realizations with high ganglia-to-pool ratios (>30) exhibit a delayed response. The two realizations with ganglia-to-pool ratios closest to the correlated ensemble mean value (6.5) behave similarly to the correlated ensemble realizations.

Ensembles of realizations employing Leverett scaling of capillary entry pressure to k values behave differently from realizations with entry pressure distributions modeled independently of k. In independent ensembles, predicted organic liquid distributions generally exhibit decreased lateral spreading, increased depth of penetration, and higher maximum nonaqueous phase liquid saturations. For the hypothetical PCE spill modeled in this study, ganglia-to-pool ratios calculated based on organic liquid saturations relative to Sor spanned a considerably wider range for the independent ensemble. Variations in these mass distribution metrics were associated with contrasts in predicted DNAPL source-zone architectures and resulted in different responses to simulated DNAPL mass removal.

In this study, correlated and independent realizations with identical ϕ and k distributions differed in their response to varying degrees of PCE recovery. However, correlated and independent realizations with different ϕ and k distributions but similar ganglia-to-pool ratios displayed similar mass flux responses to simulated PCE removal. For example, independent realizations with ganglia-to-pool ratios close to the correlated ensemble mean ganglia-to-pool ratio behaved similarly to representative realizations from the correlated ensemble (Fig. 6).

More striking was the behavior of independent realizations with low ganglia-to-pool ratios. For these realizations, decreases in aqueous phase effluent concentrations of one order of magnitude or more were predicted after only 60%–70% PCE removal (Fig. 6). Simulations based upon additional DNAPL source-zone configurations presented in Lemke et al. (2004b) suggested that realizations with still-lower ganglia-to-pool ratios responded in a similar fashion after removal of 20%–60% PCE (Fig. 7). These modeling results suggest that a low ganglia-to-pool ratio is a potentially important characteristic for DNAPL source zones because measurable decreases in contaminant mass flux can result from removal of less than 70% of the DNAPL mass. Field studies of surfactant-enhanced DNAPL solubilization have demonstrated the technical feasibility of such mass fraction removal levels when DNAPL source zones are accessible for treatment (e.g., Fountain et al., 1996; Londergan et al., 2001). Conversely, DNAPL source zones with high ganglia-to-pool ratios can be expected to require recovery of proportionately more DNAPL mass before decreases in dissolved contaminant mass flux are realized. Consequently, the ganglia-to-pool metric may be useful for distinguishing between DNAPL source zones that respond quickly in terms of mass removal per remediation effort applied (i.e., high GTP source zones) and those that yield more rapid reductions in aqueous phase mass flux per organic mass fraction removed (i.e., low GTP source zones). This observation is consistent with previously published simulation results, which illustrated the rapid removal of PCE from ganglia relative to pools (Lemke et al., 2004b, their Fig. 11) and the more rapid reduction of dissolved PCE mass flux for low ganglia-to-pool realizations in response to groundwater flow under natural gradient conditions (Lemke et al., 2004b, their Fig. 12).

For a wider range of release scenarios and levels of permeability variability than those considered here, it is likely that flux response curves could vary even within a specific ganglia-to-pool ratio category, however. Such variation would be due to realization-specific differences in the horizontal and vertical distribution of pools, reinforcing the need for detailed site characterization to engineer effective remediation designs (Fountain, 1997). Nevertheless, the model results presented here highlight the importance of characterizing the relative organic liquid mass fractions in pools and zones of residual saturation if the benefits of partial DNAPL source-zone remediation are to be assessed at specific field sites. Such characterization likely would require multiple organic liquid saturation measurements from soil cores distributed throughout the DNAPL source zone coupled with an appropriate geostatistical analysis of their spatial variability. The application of indicator geostatistics is a logical choice for ganglia-to-pool characterization because of its ability to represent non-Gaussian random cumulative distribution functions and to characterize the spatial variability of categorical (i.e., discrete valued) random-function variables (Goovaerts, 1997) that can be used to classify ganglia and pools.

Furthermore, the predicted differences in the response to partial mass removal for realizations with contrasting ganglia-to-pool ratios described in this study have important implications for stochastic modeling of DNAPL source zones in cases where detailed subsurface measurements are limited. The greater variability of predicted DNAPL distributions, expressed in the larger range of ganglia-to-pool ratios documented here for the independent ensemble, captures a wider range in the uncertainty of potential source-zone architectures. The use of correlated ensembles that rely on the application of Leverett scaling may limit model variability and thus fail to capture the full range of variability in potential responses to partial source-zone mass removal. Such variability is desirable in stochastic modeling if it is physically real. Lemke et al. (2004b) argued that independent realizations were more representative of the modeled aquifer than correlated realizations because they incorporated additional information embedded within geostatistical indicator classes derived from grain-size distribution curves. Consequently, the independent stochastic modeling approach used here is judged to be superior to the more conventional correlated approach.

Finally, the results of this study provide an opportunity to assess alternative conceptual models that postulate a priori correlations between DNAPL saturations and aquifer properties such as k or Pb. Model cells from the independent ensemble containing PCE were interrogated for organic liquid saturation, So, and aquifer properties. Cross plots of predicted organic liquid saturation versus k and Pb are presented in Figure 8. Saturation distributions do not exhibit a discernible trend that would support a positive or negative correlation between So and k or Pb. Plots of So versus vertical contrast in k and Pb values (overlying cell value minus underlying cell value) also fail to reveal a clear relationship that would allow a priori prediction of So on the basis of local k or distributions (Fig. 8). Organic Pb liquid saturations in Figure 8D only exceed Sor in cells with negative vertical Pb contrasts exceeding a discernible threshold. Nevertheless, PCE pooling cannot be predicted on the basis of such a vertical Pb contrast, because saturations accumulating in individual model cells following a simulated spill are highly path-dependent. Assumptions about the distribution of organic liquid saturations in DNAPL source zones that are tied to the distribution of k or Pb are therefore not supported by this analysis.

Differences in PCE infiltration and entrapment for a simulated PCE spill in a homogeneous, nonuniform aquifer were predicted for two ensembles of realizations modeled with identical ϕ and k distributions. The first ensemble employed Leverett scaling of capillary entry pressure to k, while the second modeled entry pressure independently of k. Decreased lateral spreading, increased vertical spreading, larger maximum organic liquid saturations, and a wider range of predicted ganglia-to-pool ratios were predicted for independent realizations. The simulated spill generated low ganglia-to-pool ratios in several of the 50 independent realizations.

Differences in predicted source-zone architecture led to a range of responses to varying degrees of PCE recovery in the models investigated. Removal of 90%–98% of entrapped PCE reduced dissolved concentrations and mass flux by approximately one order of magnitude for all of the correlated realizations and those independent realizations with similar or greater ganglia-to-pool ratios. On the other hand, independent realizations with lower ganglia-to-pool ratios (GTP < 2) displayed a more favorable flux reduction response. An approximate order of magnitude reduction of mass flux was predicted after only 20%–70% removal of the DNAPL source zone mass in these realizations. In almost every case, however, removal of 98%–99% of entrapped PCE was required to achieve a second order of magnitude reduction in dissolved contaminant mass flux.

Although the predicted one to two order of magnitude decrease in effluent concentrations is insufficient to meet maximum contaminant level guidelines for most chlorinated solvents, potential benefits of DNAPL mass reductions in this range may include decreased exposure risk to downgradient receptors or enhancement of the potential for natural attenuation. Such benefits will depend upon the DNAPL source-zone architecture, which is not readily predictable using a priori assumptions about the spatial correlation of physical aquifer parameters. Results of this study suggest that DNAPL source-zone characterization at field sites with homogeneous, nonuniform aquifers should include a determination of the overall ganglia-to-pool ratio, and that stochastic investigations of DNAPL source zones at such sites should not rely entirely upon the assignment of capillary entry pressures through Leverett scaling. Quantification of the ganglia-to-pool ratio may provide a means to predict, at least qualitatively, the ease with which contaminant mass can be recovered and the potential magnitude of mass flux reduction that can be expected in response to the initial stages of mass removal from DNAPL source zones. Additional investigations incorporating alternative DNAPL spill scenarios in aquifer systems with discrete heterogeneities, greater variance in permeability, or longer correlation lengths are required before these conclusions can be generalized to a wider range of aquifer and source-zone conditions.

The U.S. Environmental Protection Agency, Great Lakes and Mid-Atlantic Center for Hazardous Substance Research (GLMAC-HSRC) under grant no. R-825540, the Michigan Department of Environmental Quality (MDEQ) under contract no. Y80011, the Strategic Environmental Research and Development Program (SERDP) under contract no. CU-1293, and Wayne State University provided funding for this research. The content of this publication has not been subject to agency review and does not necessarily represent the views of agency sponsors. This paper was improved by the thoughtful comments and questions provided by two anonymous reviewers.

Attribution: You must attribute the work in the manner specified by the author or licensor (but no in any way that suggests that they endorse you or your use of the work).
Noncommercial ‒ you may not use this work for commercial purpose.
No Derivative works ‒ You may not alter, transform, or build upon this work.
Sharing ‒ Individual scientists are hereby granted permission, without fees or further requests to GSA, to use a single figure, a single table, and/or a brief paragraph of text in other subsequent works and to make unlimited photocopies of items in this journal for noncommercial use in classrooms to further education and science.