Unraveling the Fate and Transport of DNAPLs in Heterogeneous Aquifer Systems—A Critical Review and Bibliometric Analysis

Abstract

Dense nonaqueous phase liquid (DNAPL) in the subsurface environment beyond the permissible limit poses a threat to human health and a suite of ecological services. An accurate prediction of the concentration and mass fluxes of DNAPL at environmentally sensitive locations and their temporal variations can be obtained using robust and efficient fate and transport mathematical models. Thus, this study evaluated 412 articles published from 1990 to 2022 utilizing the Scopus^® database to provide a quantitative overview of the present trends and future perspectives of the DNAPL transport research field, especially fate and transport models via bibliometric analysis. The major findings of the published literature based on the scale of the study and type of modeling framework, relationships of governing parameters with a scale of study, and recent developments in the mathematical models were discussed. The country-citation analysis revealed the USA and Canada as leading countries in DNAPL fate and transport research field. The findings of this study uncovered a need for studies considering low-permeability and stagnant regions, as well as the variable nature of the flow, transport, and reaction parameters to understand the complex plume evolution dynamics of DNAPLs under field-scale conditions. The nonlinear interactions between various flow and transport phenomena should be quantified during a model’s development by applying spatial- and time-varying global sensitivity analyses. The outcomes of the bibliometric survey, visual analysis, and concise review presented in this study can provide a wide range of references, emerging topics, and prospects to emphasize less focused on topics of DNAPL transport research.

1. Introduction

Dense nonaqueous phase liquids (DNAPLs) are chemical compounds detected at several hazardous waste disposal sites [1]; for example, tetrachloroethene (PCE) and trichloroethene (TCE) have been found in Superfund sites in the United States of America [2]. Studies have reported that the complete isolation or removal of DNAPL from a source region is difficult and can still leave residual concentrations above the acceptable limit (e.g., 0.005 mg/L for PCE and 0.005 mg/L for TCE), which causes the long-term contamination of aquifer systems [3,4,5]. The accumulation of residual DNAPLs in complex heterogeneous systems, such as aquifers with low-permeability, porous media (LPPM), and/or the storage of DNAPLs in the sorbed and dissolved phases has been reported in several studies [3,6]. The fate and transport of DNAPL and its dissolved components in saturated porous systems are governed by physical, chemical, and (bio)geological processes, such as dissolution, degradation, dispersion, molecular diffusion, mass transfer among different phases, and sorption–desorption processes [7,8,9]. The complexity of DNAPL plume evolution in the dissolved phase increases with the inclusion of dead-end regions or LPPM (e.g., clay/silt lenses) in a heterogeneous aquifer system [10,11]. It is observed that the LPPM, or aquitard layer, behaves as a sink during the DNAPL’s loading period and as a source during the offloading or isolation periods [4,12,13]. An efficient remediation plan for any DNAPL-contaminated site can be achieved when factors and processes governing the plume’s evolution are known from the lab- to field-scale levels. Numerous studies have been carried out over the last 20 years to investigate the transport behavior of DNAPLs (chlorinated solvents) in heterogeneous porous systems. These studies were related to the architecture of the source zone, pools, and fingers of the DNAPL [14,15,16]; mass transfer from the source region to the downstream region [8,17]; and physical, chemical, and biological processes governing the plume evolution [13,18]. Several review studies addressing the various aspects of the fate and transport of organic contaminants in the porous media from a lab to field scale [19,20,21] have been published, along with comprehensive reports by government authorities [22,23].

It is reported that the asymmetric concentration profiles from contaminant (tracer-aqueous phase) transport experiments carried out for sandy aquifers were majorly dominated by the spatial variability of the horizontal and vertical hydraulic conductivities [24]. Furthermore, the distance and/or time dependence of the dispersion coefficient/dispersivity was observed in several laboratory- and field-scale studies [25,26,27]. Smoothed particle hydrodynamics simulations highlighted the dominance of mixing over channeling when the flow was parallel to the major axes of elliptical grains, thus showing the significance of microscale heterogeneities on the plume evolution, even for passive tracers [28]. Muniruzzaman and Rolle [29] concluded that microscale physical processes, such as molecular diffusion and dispersion, significantly affect the transport dynamics at the macroscale for homogeneous and heterogeneous porous media. A review highlighting the influence of heterogeneity and low-permeability, porous media on the transport of DNAPLs could provide a more convincing representation of plume evolution.

This study aims to provide a bibliometric analysis and an overview of the governing flow and transport mechanisms and recent developments in the modeling of DNAPL transport in heterogeneous porous systems. It (1) takes into account the advancements in the DNAPL transport research field from 1990 to 2022 via a bibliometric analysis, (2) identifies the leading affiliations and countries and their correlation with the spatial distribution of DNAPL-contaminated sites, and (3) identifies thematic areas and their temporal evolution and categorizes the studies based on the scale and type of modeling framework implemented and, furthermore, summarizes their major findings.

2. Methodology

The overall methodology adopted for the bibliometric survey and literature review is presented in Figure 1. In the first step, a bibliometric and visual analysis based on the Scopus database was carried out. In the next step, the studies were classified based on the scale and mathematical modeling. Further, major governing factors and mechanisms were analyzed and reviewed to interpret the DNAPL transport in the porous systems in the presence of low-permeability, porous media.

2.1. Bibliometric Analysis Steps

Bibliometric analysis is an effective method for the statistical quantification of the research impact of a particular research field [30]. Thus, a bibliometric survey was conducted based on the published literature on DNAPL transport in the subsurface research field.

2.1.1. Database Selection

The Scopus database was found to be the largest abstract and citation database [31] and offered 20% higher coverage than Web of Science [32]. Thus, research articles, review papers, and book chapters published between 1990 and 2022 were obtained from the Scopus^® database and used in the bibliometric analysis. The steps implemented to retrieve the data on DNAPL transport in an aquifer and/or subsurface system are shown in

Table 1. Steps followed for data collection.

Search Steps	Steps Used in Search	Query on Scopus	Description
1	Search Keywords	TITLE-ABS-KEY	(((“Chlorinated Solvents” OR “DNAPL” OR “DNAPLS” OR (“Dense Non-Aqueous Phase” AND (“Liquid” OR “Liquids”))) AND (“Transport”) AND (“in” AND “the”) AND (“Groundwater” OR (“Ground” AND “Water”) OR “Aquifer” OR (“Saturated” AND “porous” AND “media”))))
2	Year	AND PUBYEAR > 1990	Published from 1990 to 2022 considered
3	Document Type	AND (LIMIT-TO (DOCTYPE, “ar”) OR LIMIT-TO (DOCTYPE, “re”) OR LIMIT-TO (DOCTYPE, “ch”))	Journal article, review paper, book chapter (journal articles, review papers, and book chapters searched, ignoring conference proceedings, etc.)
4	Document’s Language	AND (LIMIT-TO (LANGUAGE, “English”))	Language “English”
5	Subject Area	AND (LIMIT-TO (SUBJAREA, “ENVI”) OR LIMIT-TO (SUBJAREA, “EART”) OR LIMIT-TO (SUBJAREA, “CHEM”) OR LIMIT-TO (SUBJAREA, “ENGI”) OR LIMIT-TO (SUBJAREA, “AGRI”) OR LIMIT-TO (SUBJAREA, “CENG”) OR LIMIT-TO (SUBJAREA, “COMP”) OR LIMIT-TO (SUBJAREA, “MATE”) OR LIMIT-TO (SUBJAREA, “MATH”) OR LIMIT-TO (SUBJAREA, “MULT”))	Limited to Environmental Science, Earth and Planetary Sciences, Chemistry, Engineering, Agricultural, and Biological Sciences, Chemical Engineering, Computer Science, Material Science, Mathematics, and Multidisciplinary

.

2.1.2. Bibliometric Indicators and Tools Used

In this study, conventional bibliometric indicators, such as (i) trend of the number of publications, (ii) total citations per year, and (iii) total citations per article, were computed initially using Biblioshiny R-toolbox [33]. Then, the leading authors, most cited articles, most dominant journals, leading countries, and most productive affiliations of DNAPL transport in the subsurface were determined. Further, the correlations of the leading organizations and leading countries with the spatial distribution of DNAPL-contaminated sites were analyzed. A three-field plot (i.e., Sankey diagrams), co-occurrence network visualization of the author’s keywords and cited references, and thematic maps were obtained via Biblioshiny R-toolbox and VOSviewer [34].

3. Findings of Bibliometric Analysis

The preliminary information on the publications retrieved from Scopus on DNAPL transport in the subsurface is shown in

Table 2. Main description of DNAPL transport reviews from the Scopus^® database.

Description		Results
Timespan		1990:2022
Sources (Journals, Books, etc.)		107
Number of Research Documents		412
Average Years from Publication		13.60
Average Citations per Research Document		31.59
Average Citations per Year per Research Document		2.28
References		15,507
Document Types	Research Article	397
	Book Chapter	6
	Review	9
Document Contents	Keywords Plus (ID)	3625
	Author’s Keywords (DE)	882
Authors	No. of Authors	1060
	No. of Appearances of Authors	1531
	No. of Authors Who Published as Single-Authors	23
	No. of Authors Who Published Multi-Authored Documents	1037
Authors Collaboration	No. of Single-Authored Documents	25
	Documents per Author	0.389
	Authors per Document	2.57
	Co-Authors per Document	3.72
	Collaboration Index	2.68

. A total of 107 sources from the Scopus^® database were observed, which published 412 research documents from 1990 to 2022. Detailed information on several publication statistics is provided in Table 2. The year-wise variation in the publication statistics is shown in Figure 2. The number of publications increased significantly between 1997 and 1998, then increased steadily from 1998 to 2007, followed by a sharp decline in the number of articles (NoA) from 24 in 2007 to 7 in 2008 (Figure 2). A total of 21 and 24 articles on DNAPL transport were published in 2006 and 2007, while 14 articles were published in 2021. In addition, a fluctuating trend was observed in the mean total citations per article and mean total citations per year. Therefore, understanding the trends in the citations of published articles is important due to the varying trend in the number of publications with time. The highest value of the mean total number of citations (TCi) per article of 75.6 was observed in 1997; however, the mean TCi per article suddenly dropped to 22.87 in 1998. A fluctuating trend in the mean TCi per article was observed up to 2015, and then a decreasing trend in the mean TCi per article was observed from 2016 to 2022. The highest value of 7.69 for the mean TCi per year was observed in 2013, which dropped to 1.29 in the year 2014. From 2016 to 2022, an increasing trend in the mean TCi per year was observed. A nonuniform variation in the mean TCi per year has been observed over the past thirty years.

The variation in the number of publications and cumulative publications with time for the DNAPL transport in the subsurface system research field is shown in Figure 3. It is observed that an average of 5 articles were published from 1990 to 1997, which further increased to approximately 15 articles from 1998 to 2007 and dropped to 7 articles in 2008. In addition, it is observed that at least 14 articles related to DNAPL transport were published for most of the years from 2009 to 2021. A mild slope for the cumulative publications with time was observed from 1990 to 1997; however, an increasing slope was observed after 1998, showing the increased attention of researchers and funding organizations in DNAPL transport in the subsurface research field post-1997. Overall, the trends of the mean TCi per year and cumulative publications with time show the importance of the DNAPL transport research field and highlight the need to conduct such studies in the near future.

Contribution of Countries and Spatial Distribution of DNAPL-Contaminated Sites

The contribution of countries towards the DNAPL transport research field was quantified based on bibliometric survey-based parameters, such as the number of articles (NoA) published, total number of citations (TCi), and average article citations (AACi), as shown in Figure 4. In the USA, 173 articles were published, and it was observed as the leading country in the DNAPL transport research field (Figure 4a). The USA received a total of 6856 citations, as shown in Figure 4b, and ranked in the first position on the basis of TCi and NoA (Table S1 in the Supplementary Materials). However, the USA ranked fourth based on the AACi criterion, with a value of 39.63 (Figure 4c). For Canada, 42 articles were published, it received 2130 citations, and it was observed in the first position among the top ten leading countries on the basis of the AACi criterion. The United Kingdom, with 506 total citations and 24 articles, ranked in third place based on the TCi and NoA criteria; however, based on AACi, it ranked in the eighth position on the list. The USA received more than three times as many citations as Canada, showing that the DNAPL transport field from 1990 to present is led by the USA. Australia, with 247 TCi and 6 NoA pertaining to DNAPL transport, ranked ninth and tenth on the list, but based on an AACi value of 41.17, it ranked third (Table S1). Furthermore, the existence of any correlation between leading countries and the spatial distribution of DNAPL-contaminated sites was investigated. A highly DNAPL-contaminated site was reported at a Belfast industrial facility in Northern Ireland, with the highest value of $\mathrm{390,000}\mathsf{\mu }\mathrm{g}/\mathrm{L}$ of TCE [35], as shown in Figure 5. The maximum value of $\mathrm{250,000}\mathsf{\mu }\mathrm{g}/\mathrm{L}$ of TCE was also reported at a Canadian Forces base site in Borden, ON, Canada, while $\mathrm{170,000}\mathsf{\mu }\mathrm{g}/\mathrm{L}$ of dichloroethylene (DCE) was found in an industrial area of the Cape Canaveral Air Station, Cape Canaveral, FL, USA. [35]. Overall, it is observed that most of the DNAPL-contaminated sites were in the USA and Canada (Figure 5 and Table S2 in the Supplementary Materials), due to which most of the advancements in the DNAPL research field were led by organizations located in these countries. A detailed analysis of the findings of the bibliometric survey is presented in the Supplementary Materials.

Furthermore, the contributions of organizations towards the DNAPL research field were studied. The University of Waterloo ranked first based on the NoA criterion (39), followed by the University of Arizona with 28 articles, and the University of California ranked third with 14 articles (Table S3 of the Supplementary Material). It can be seen that 3 of the top 10 affiliations are from Canada, while 7 of the top 10 affiliations are from the USA, showing that the two countries dominate the DNAPL transport and modeling research field. Based on the analyzed bibliometric data, we found that the transport behavior of DNAPLs in aquifer systems dominated a large portion of the contaminant hydrogeology research field. Thus, the present study emphasizes the DNAPL transport behavior in heterogeneous aquifer systems.

4. What Controls the Dissolved-Phase Plume Evolution of DNAPLs in Heterogeneous Porous Systems?

DNAPLs, when spilled over a land surface, migrate downward dominantly through unsaturated regions and spread in longitudinal and transverse directions in the saturated zone below the water table, as shown in Figure 6. The spread and formation of DNAPL pools take place in the saturated zone if there exist any low-permeability regions (LPPMs), such as clay/silt lenses [39,40,41]. Thus, in this section, the potential physical, chemical, and biogeochemical governing factors controlling the evolution of dissolved-phase DNAPL plumes are discussed. Further critical observations are drawn from studies emphasizing DNAPL transport in heterogeneous porous systems.

4.1. DNAPL Fate and Transport under Dynamic Groundwater Flow

Groundwater flow gradients were observed to be as equally important as that of hydraulic conductivity for understanding the transport behavior under highly dynamic flow conditions at the Hanford Integrated Field Research Challenge (IFRC) field site [42]. It is observed that the position, architecture, and geometry of TCE can be changed even at the lowest flow velocity of $0.05\mathrm{m}/\mathrm{d}\mathrm{a}\mathrm{y}$ from laboratory-scale ($1.0\mathrm{m}\times 0.12\mathrm{m}\times 0.60\mathrm{m}$) experiments [14]. The displacement of a DNAPL pool ranging from $0\mathrm{m}$ at velocity at $v_w=0.0\mathrm{m}/\mathrm{d}\mathrm{a}\mathrm{y}$ to $1.45\mathrm{m}$ at $v_w=40.0\mathrm{m}/\mathrm{d}\mathrm{a}\mathrm{y}$ was observed, depicting the impact of high groundwater velocities. However, the study was limited to laboratory-scale conditions [14]. Similarly, a hydraulic gradient was stated as a predominant factor for the infiltration of DNAPL and promoted the downward and down-gradient migration, as observed from a 2D laboratory tank ($1\mathrm{m}\times 0.12\mathrm{m}\times 0.7\mathrm{m}$) study [43]. It is reported that an increase in the hydraulic gradient from $0.00$ to $0.08$ enhanced the infiltration rate and decreased the amount of remaining residual saturation [43]. Yang et al. [44] observed hydraulic gradients from 0.014 to 0.025 near the main source area of a TCE spill, indicating a strong seasonal recharge, which further resulted in the leaching of entrapped TCE from the unsaturated region. In contrast, low hydraulic gradients from 0.006 to 0.008 were observed in the industrial complex of Wonju city, Korea, depicting a limited recharge. A strong correlation between TCE concentrations and the change in hydraulic head ($∆h$) was observed near the source transect. An increase of 2 m of the hydraulic head caused a yield of 20,000 µg/L of TCE concentration; however, the influence was not significant for locations away from recharge area [44].

In a study carried out by Zheng et al. [15], heterogeneity was revealed as a dominant factor affecting the transport behavior of low-viscosity DNAPL (i.e., PCE) when the simultaneous influence of enhanced flow velocity and heterogeneity (${\sigma }_{lnk}$) was considered. An increase in the standard deviation of the log permeability (${\sigma }_{lnk}$) from 0.20 to 1.0 caused an increase by 207 and 633% in the variability of the spatial variances ${\sigma }_{xx}$ and ${\sigma }_{zz}$, respectively [15]. Spreading in a horizontal direction increased, and the infiltration rate decreased with the inclusion of layered, fine sand lenses in coarse sand, porous media, as the layered lenses promoted the lateral spread above the lenses [15]. These studies were limited to homogeneous systems such as sand layers with grain sizes $0.5–1.0\mathrm{m}\mathrm{m}$ and neglected the impact of LPPM (clay/silt lenses) on the transport behavior of DNAPL. Thus, it can be stated that there is a need for studies emphasizing heterogeneity and fluctuating groundwater table conditions to understand the fate and transport of DNAPLs.

4.2. Mass Transfer Processes between Aquifer and Low-Permeability, Porous Media

DNAPL plumes can persist in the heterogeneous soil-water system depending upon forward and backward diffusion through mobile and immobile regions [45,46]. While implementing the physical practices like flushing or after source removal, immobile contaminants captured in dissolved or sorbed phase in the LPPM either diffuse out due to the concentration gradients [47,48,49] or causing elongated plume for larger duration due to the fact of desorption [45,50]. In the study by Mobile et al. [51], depth dependence of mass transfer coefficients ($0.082-2.0\mathrm{d}\mathrm{a}{\mathrm{y}}^{-1}$) was observed for the DNAPL laboratory-scale experiments ($5.5\mathrm{m}\times 4.5\mathrm{m}\times 2.2\mathrm{m}$) conducted in sandy aquifer conditions. The differences in the maximum peak concentration and mass recovery were observed even for nonreactive tracer in saturated column condition, showing the effect of molecular diffusion coefficient [52]. It was observed that the nonideal transport behavior for a large heterogeneous porous system of Superfund site was mainly associated with mass transfer coefficient constraints from stagnant zones [49,53].

Nonideal mass-removal was observed even for homogenous soil-water systems, attributed solely to the arrangement of wells and mass transfer constraint [49] and back diffusion from captured masses in the heterogeneous system [54]. Tatti et al. [12] conducted laboratory-scale box model tracer (fluorescein) experiments and observed that remediation timeframe of contaminated aquifer region ($permeability\approx {10}^{-3}\mathrm{m}/\mathrm{s}$) consisting of low-permeability lenses ($permeability\approx {{10}^{-10}-10}^{-3}\mathrm{m}/\mathrm{s}$) decreased slowly and was independent of flow velocity. On contrast to observation given by Tatti et al. [12], orthogonal flow to the contaminated LPPM lenses was suggested as a measure for a faster diffusion rate at the interface between LPPM and aquifer region [55]. It is stated that the dissolved-phase by-products of DNAPL transfer between mobile and immobile regions of the soil-water system depends upon internal hydraulic gradients between contrasting permeability materials [56].

The temporal variation in the interphase mass transfer coefficient was observed from TCE flushing tank ($L80\mathrm{c}\mathrm{m},H40\mathrm{c}\mathrm{m},\mathrm{a}\mathrm{n}\mathrm{d}W5\mathrm{c}\mathrm{m}$) experiments [57]. Using diffusion coefficients measured in the lab, it was found that the mass storage of TCE ($35.0–89.4\mathrm{g}$) was $25–65\%$ of the mass ($137.7\mathrm{g}$), calculated using diffusion coefficients observed in the in situ field conditions [58]. The enhanced mass storage observed in the field was primarily caused by DNAPL that was stored in the cracks and diffused out through diffusion [58]. In addition, the Pearson correlation between the thickness of the LPPM (i.e., aquitard layer) and the effective diffusion coefficient was determined based on data compiled from various studies (Table S4 in the Supplementary Materials). The aquitard thickness varied from 0.015 m for a flow chamber experiment to 220 m for a field-scale study. Effective diffusion coefficient ($D^{*}$) was observed as −0.3603 for aquitard thickness which was negatively correlated.

4.3. Dissolution

Dissolution is considered a time-bound mass exchange reaction process [7]. It depends on several factors, such as the solubility of the DNAPLs, actual and solubility concentrations, contact area between DNAPL and groundwater, and flow conditions [1]. The aqueous phase solubility level of the DNAPL and dissolution-related mass transfer were identified as critical parameters in the absence of sorption processes based on a study by Clement et al. [59], which investigated entrapped DNAPL behavior in a laboratory aquifer model. Long-term persistence ($>50\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}\mathrm{s}$) of dense chlorinated solvents (TCE) was observed with a low dissolution rate due to the fact of a pool carrying $~6000\mathrm{k}\mathrm{g}$ of solvent [60]. The dissolution rate of the contaminant mass from fingers was suspected to occur faster than for pools, due to which aqueous phase concentrations were governed predominantly by the dissolution of DNAPL from pools [61]. Therefore, dissolution can be recognized as a dominant factor for the DNAPL plume evolution based on studies [60,61]. Further, the combined influence of the dissolution of DNAPL and matrix diffusion on the source zone longevity was investigated [17]. It was observed that the dissolution and matrix diffusion contributed equally to the source zone longevity for lower solubility chlorinated solvents (tetrachloroethene and PCE), while 97% of the longevity was found to be contributed by matrix diffusion only for higher solubility solvents such as dichloromethane [17].

DNAPL can dissolve in heterogeneous soil-water system under the influence of injection flow rate, distance to observation wells, and hydrogeological conditions, found to be governing the travel time of dissolved by-products [62]. Dissolution caused a large amount of mass discharged from the source zone during the loading period, which further resulted in a higher amount of back-diffusion flux for a shorter duration from an aquitard (low-permeability region) to an aquifer and finally prolonged the risk [45]. It was observed that an integrated scenario of dissolution and back diffusion represented the monitored PCE and TCE concentrations from matured site (Dover Air Force Base in Delaware, a contaminated site in Connecticut, and Naval Air Station Jacksonville in Florida) wells. On the other hand, back diffusion solely described the observed data for a site where a contaminant source was isolated [48]. However, it was observed that the degradation reactions and diffusion in the aquitard layer along the transverse direction were ignored. Thus, there is a need for studies in which degradation in the LPPM is considered.

4.4. Sorption–Desorption in the Aquifer (Mobile) and Aquitard (Immobile Type) Region

The long-term risks due to the fact of back diffusion from highly sorbing DNAPLs (retardation factor = 11) were much higher compared to conservative contaminants [10]. An increase in the contact time resulted in an increase in the elution tailing of trichloroethene (TCE) transport, as observed in the $7\mathrm{c}\mathrm{m}$ long column experimental set-up by Brusseau et al. [50]. This tailing was found to be primarily attributed to rate-limited desorption [50,63], while a minor contribution was observed by nonlinear sorption [63]. The depth of penetration of the contaminant in the aquitard was found to be inversely dependent on sorption; however, the long-term risks due to the fact of back diffusion were observed to be higher for sites with highly sorbing DNAPLs [45]. An increase in the transfer of 4-Cl-Nb mass into an immobile region was observed in comparison to a conservative contaminant, which depicted the influence of sorption on the clay matrix [64]. The retardation factor representing the sorbed mass in the LPPM region (silt/clay layer), geometry of the LPPM region (thickness of silt/clay lenses), and groundwater velocity in the adjoining higher permeability layer were found to be important factors influencing back-diffusive fluxes and the remediation timeframe [12,65]. It is stated that the plume persistence near the DNAPL source location was affected by nonideal sorption following its depletion [16]. Large variations in the plume were observed for different sorption models [16], indicating that the choice of the sorption model is crucial for an precise forecasting of DNAPL plume longevity. However, chemical and biological transformations were not considered in the modeling study, as daughter products could contribute to plume persistence [16].

4.5. Degradation Processes

Degradation is the process that breaks down DNAPL into less harmful substances and depends on factors such as the type and concentration of DNAPL, properties of the porous media, and environmental conditions. Thus, in this subsection, studies related to the degradation of DNAPLs are discussed. The influence of biodegradation (aerobic and anaerobic) reactions along with diffusive mass transfer on the plume persistence of methyl tert-butyl ether (MTBE)/tert-butyl alcohol (TBA) was studied at Vandenberg Air Force Base (VAFB) site [6]. The results revealed that degradation reactions decreased the long-term persistence of MTBE and TBA in sandy aquifers by lowering the mass available for back diffusion [6]. Compound-specific isotope analysis (CSIA) profiles depicted the feasibility of the degradation reactions of PCE, TCE, and chloroform in the clayey aquitard layer, irrespective of small pore sizes [41]. Wanner et al. [13] observed that the TCE back diffusion periods were shorter for aquitard-degradation scenarios in comparison to a nondegradation scenario. Long-term tailings in the breakthrough curves of daughter species (cis-DCE) were observed for scenarios representing high degradation rates ($half-life:30–80\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$), while for slow degradation rates ($half-life:<200\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$), aquifer contamination due to the parent (TCE), as well as daughter (cis-DCE), species was observed [13]. The aquitard depth was found to be a dominating factor affecting plume persistence compared to the concentration of daughter species [13]. Thus, it can be summarized that a nonuniform variation in the degradation rate (e.g., decreasing degradation with aquitard depth) can mimic the contaminant transport (chlorinated solvent) behavior under actual field conditions efficiently compared to conventional constant degradation rate scenarios [41].

4.6. Transverse Dispersion and Mixing

Transverse dispersion is a critical factor for plume spreading, and mixing [66,67,68,69,70] is still not emphasized compared to longitudinal dispersion in several existing numerical and analytical studies. Back diffusion, horizontal advection, and vertical transverse dispersion were revealed as the major reasons for trichloroethylene (TCE) plume tailing in the aquifer region, with concentrations ranging from $200$ to $2000\mathsf{\mu }\mathrm{g}/\mathrm{L}$ ($>MCL=5\mathsf{\mu }\mathrm{g}/\mathrm{L}$), after DNAPL source isolation [4]. The transverse vertical dispersivity and length of the silt/clay layer were found to be influencing operational variables for pump and treat operations based on simulations carried out in the developed in situ Remediation-MT3DMS model for field-scale case studies [65]. Transverse mixing of the conservative contaminant (e.g., bromide and uranine) was observed to be governed majorly by contrasts in hydraulic conductivity between fine matrix ($0.20-0.30\mathrm{m}\mathrm{m}$ grain size) and high permeability lenses (silica beads with a grain diameter of $1.0-1.5\mathrm{m}\mathrm{m}$) in comparison to actual pore-scale dispersivity, as revealed from 2D tank ($77.30\times 14.0\times 1.1\mathrm{c}\mathrm{m}$) experiments [71]. A study by Muniruzzaman and Rolle [29] concluded that microscale physical processes such as diffusion and electrochemical processes (e.g., Coulombic interactions) significantly affected the transport both in the homogeneous, as well as heterogeneous, porous media at the macroscale.

A comparison between ${\alpha }_{longitudinal},{\alpha }_{tranvserse},$ and ${\alpha }_{lateral}$ was also conducted on the basis of data collected from published studies to analyze the variation of the dispersivity (Table S5 in the Supplementary Materials). The average value of ${\alpha }_{longitudinal}=1.705\mathrm{m}$ was observed, which was ~13 times higher than the average ${\alpha }_{tranvserse}$ (${\alpha }_{tranvserse}=0.128\mathrm{m}$); furthermore, an average value of $0.032\mathrm{m}$ was observed for ${\alpha }_{lateral}$ (Figure 7). A noticeable gap between dispersivity values along different directions was revealed, indicating the significance of high contrast in the dispersivity values used in the modeling. It can be stated that there is a need to emphasize transverse dispersion/dispersivity in the simulation of DNAPL transport through a heterogeneous aquifer system.

5. Overview of Studies Pertaining to DNAPL Plume Evolution in Heterogeneous Porous Systems

Studies emphasizing DNAPL transport in subsurface systems can be categorized into (i) studies that focus on the source zone formation composed of DNAPL release and migration of different phases, (ii) studies that emphasize mass transfer processes between DNAPL and the aqueous phase, and (iii) studies related to DNAPL transport down-gradient of the source zone involving retardation, degradation, and mass transfer between contrasting permeability zones [20]. For example, studies analyzing the transport dynamics in the porous systems comprising low-permeability, porous media (LPPM), and adjoining aquifer regions come under the third category. In the present study, emphasis is given to studies in the third category in which the fate and transport behavior of DNAPLs were investigated in the region away from the source zone.

5.1. Classification of Studies Based on the Scale

5.1.1. Laboratory-Scale Studies

Luciano et al. [43] assessed the impacts of the hydraulic gradient and layered porous structure on the infiltration and redistribution of DNAPLs using image analysis for 2D laboratory-scale tank experimental conditions. Mass transfer coefficients were quantified by conducting multicomponent NAPL dissolution experiments by mimicking the field-scale heterogeneous source zone architecture [51]. The effect of the sorption contact time on the desorption and tailing of the concentration breakthrough curve was investigated using batch sorption and column experiments [50]. Trichloroethene (TCE) was used as a model organic compound in the experimental study, and the input pulse duration was varied as 42 min, 1.4 h, 2.8 h, 7 h, 35 h, and 15 days to study the contaminant transport behavior for different pulse durations [50]. Maruo et al. [76] evaluated the effects of well set-up and distances between injection and production wells on tracer concentrations using vertical 2D laboratory-scale experiments. Similarly, 2D bench-scale experiments were conducted to investigate the influence of heterogeneities on the transverse mixing of a conservative contaminant in a porous system [71]. The joint effect of the enhanced flow velocity and heterogeneity on the transport behavior of low-viscosity DNAPL was also analyzed [15]. The use of an ethanol–water mixture as a cosolvent for the enhancement and mobilization of DNAPL (TCE) in a saturated porous system was studied in a vertical flushing tank for different cosolvent contents, flushing velocities, and injection pattern [57].

The influence of groundwater velocities on contaminant back diffusion from LPPM was studied using a laboratory box model and subsequently via image analysis [12]. A groundwater circulation well system (GCW) was examined as remediation technology compared to a traditional pumping system for contaminated porous systems having LPPM via laboratory tests and numerical modeling [55]. An aquifer model was reconstructed using a 2-D tank set up to study the impact of DNAPL saturations, hydraulic conditions, and source configurations on DNAPL dissolution rate [77]. The effective diffusion coefficient of TCE was measured in laboratory experiments for three types of LPPM, viz., silt, silt–clay confined, and silt–clay expanded [78]. The mass storage of TCE with and without cracks in a 5 m thick low-permeability clay layer was assessed after 30 years of diffusion [58]. The performance of natural zeolite as a physical permeable reactive barrier was assessed in a bench-scale sand tank model for removing the BTEX compounds from groundwater [79]. Furthermore, the details related to the dimensions of the experimental set-up are shown in Table S6. It is observed that the length of the tank in several studies was chosen in a range from 0.60 to 5.5 m, while the width/thickness was kept from ~0.01 m to 4.5 m, and the height of tank was kept in the range from 0.14 to 2.2 m. Most studies emphasized the longitudinal and vertical migration of NAPL dissolved phase contaminant; however, they neglected the influence of the transverse direction in a horizontal plane. In most of the laboratory scale studies, the length was kept at ~1 m, while the height of the tank was kept at ~0.5 m (Table S6 in the Supplementary Materials). It should be noted that for the case of an ideal 2D modeling scenario, one of the dimensions of the tank can be 1–2 orders of magnitude less than the other dimension, as in studies by Ballarini et al. [71] and Zheng et al. [15], to actually mimic the experimental results in the 2D mathematical model. Thus, it can be inferred that most lab-scale studies were limited to a few meters and/or O (m), which reveals the need to analyze the fate and transport of contaminants at a large scale of O (Km), mimicking field conditions. Therefore, in the next subsection, a review of field-scale studies is presented.

5.1.2. Field-Scale Studies

A summary of field studies that examined DNAPL transport in the subsurface system is provided in this subsection. In early 2000, the back diffusion of DNAPL from an aquitard (clayey silt layer) to an adjoining aquifer (sand layer) region was observed for six years from an industrial site in CT, USA [4]. Temporal variations in the volatile chlorinated solvent (TCE, PCE, cis-DCE, and TCA) concentrations at 23 sites were assessed by calculating the percentage concentration change and by determining the trend using the Mann–Kendall method [80]. The impact of the seasonal rainfall variations and spatial variations in the geological conditions on the spatial and temporal variations in the TCE plume was analyzed for the Wonju industrial site in Korea [44]. Wanner et al. [41] effectively quantified the degradation of chlorinated hydrocarbons in LPPM via high-resolution compound-specific isotope analysis (CSIA). Groundwater flow scenarios at field-scale conditions were simulated using MODFLOW, and the persistence of dissolved contaminant plume in a layered system was investigated by coupling MT3DMS contaminant transport code with the random-walk-based method RWheat [49]. Similarly Guo and Brusseau [53] examined the mass removal processes from heterogeneous large groundwater systems using a 3D numerical model. The simulation capabilities of a multi-rate mass transfer (MRMT) model coupled with random walk particle tracking (RWPT) code were evaluated for field-scale conditions under steady and transient flow conditions by considering the heterogeneous domain using the geostatistical method-based program T-PROGS [56]. In the next section, the details of various developments in mathematical models and applications to simulate the DNAPL transport behavior are presented.

5.2. Classification of Studies Based on the Type of Modeling Framework Used

5.2.1. Analytical Modeling Studies

An analytical modeling framework includes mathematical models, usually referred to as closed-form models, which provide solutions in concise notation using well-known mathematical functions/formulas. An analytical modeling-based study utilized contaminant source mass-release models to study the fate and transport of DNAPLs in porous systems. Details of the several source mass release models and critical findings observed from these studies are presented in Table S7. Analytical studies implement integral, differential, and combinations of transformation techniques to solve contaminant transport equations, which eventually represent various physical and biogeochemical processes in porous systems [81,82]. For example, mass transfer processes between mobile and immobile water, residual NAPL and water, and aquifer solid and water were considered in the formulation of the simple analytical model [9]. The influence of the pools and fingers of chlorinated solvents on the properties of down-gradient contaminant plumes was investigated [61]. The combined effects of aquifer hydrodynamic heterogeneity, nonequilibrium dissolution, and NAPL source zone heterogeneity were represented well in a single parameter, namely, reactive travel time variability [83].

An aquitard diffusion model, based on dissolution and diffusion from DNAPL source, was developed to study the effects of source architecture and back-diffusive fluxes from LPPM regions [45]. Forward and backward diffusion through an aquitard were considered while deriving analytical solutions for aquitard–aquifer systems [10]. Furthermore, mobile and immobile water regions in the saturated column experiments were quantified using a single fissure diffusion model (SFDM)-based analytical approach [52]. The fate of reactive 4-Cl-nitrobenzene (4-Cl-Nb) in a well-controlled double-region porous system was quantified using SFDM-based modeling and a multitracer approach [64]. Step-change, linear, and exponential contaminant source depletion models were considered utilizing an analytical approach to examine the dynamic aquitard–aquifer diffusive mass transfer for nonreactive contaminants [48,84]. Aquitard concentration profiles were also evaluated by Yang et al. [48] using a coupled analytical–numerical model assuming three scenarios, viz., dissolution only, dissolution and back diffusion, and back diffusion only.

Analytical solutions were derived to analyze the diffusive flux exchange between aquifer and a single aquitard layer with the thickness varying from centimeters to hundreds of meters [85]. One-dimensional analytical solutions were derived to analyze the impacts of an exponential depleting source and back diffusion from LPPM on long tailings of BTCs [11]. However, the abovementioned analytical and semi-analytical solutions could only be applied to 1D saturated porous systems [11,48,84,85]. In addition, several semi-analytical/analytical models were coupled to numerical transport simulators to eliminate the limitations of analytical approaches. A spreadsheet-based modeling tool centered on diffusion, vertical advection, and degradation processes was developed for reconstructing the contaminant source history of chlorinated ethenes-contaminated sites [86]. The soil concentration versus depth data of TCE and PCE from LPPM or low-permeability zones was used in the analytical model [86]. The results revealed the negligible variation of the source concentration for sites where DNAPLs were present; however, declining source histories for sites where contaminant sources were removed, or where natural attenuation had occurred [86].

5.2.2. Semi-Analytical Modeling Studies

A semi-analytical modeling framework provides a solution that cannot be described with widely used mathematical functions/formulas and requires some form of computational technique such as numerical integration to obtain a solution. For example, a semi-analytical modeling framework, namely, FT-MIKSS, was developed by Atteia et al. [87] to approximate the dissolved-phase plumes of benzene, toluene, ethylbenzene, and xylenes (BTEX) in heterogeneous aquifers. The results obtained from FT-MIKSS were nearly identical to those predicted by numerical models (RT3D, PHT3D, and PHAST) for the nonreactive contaminant. However, FT-MIKSS demonstrated a computational speed that was 100 to 1000 times faster than the tested numerical models (RT3D, PHT3D, and PHAST) [87]. Additionally, using the FT-MIKSS model combined with sorption term, the reductive dechlorination and oxidation of chloroethene in the heterogeneous aquifers were simulated [88]. The approximations from FT-MIKSS matched well with other models (RT3D, PHT3D, and PHAST) for simple cases; however, for complex field cases, the results obtained from FT-MIKSS varied compared to the PHT3D, RT3D, and PHAST simulations [88]. The FT-MIKSS approach was proved as a valuable preliminary visualization tool for the development of complex numerical models and for faster fitting of transport parameters for field studies [88].

Numerical and semi-analytical fate and transport models were combined to simulate matrix diffusion [89,90]. The developed method worked well during the contaminant loading period, with the run periods varying from fractions of seconds to half a minute because a matrix-diffusive flux term was added as a concentration-dependent term to the discretized equations [89]. A comprehensive toolkit, known as ICET³ (Integrated Contaminant Elution and Tracer Test), was developed to enhance the characterization of mass-removal behavior and accurately depict the mass transfer processes that influence the transport of DNAPLs in complex soil-water systems [91]. It is observed that information related to site-specific processes and values of parameters from ICET³ would enhance the efficiency of a mathematical model, support the risk assessment framework, and promote cost-effective remediation and management procedures [91]. In most studies [11,48,84,85], only diffusion was assumed for contaminant transport in the aquitard or LPPM region. Therefore, it can be stated that a hybrid analytical-numerical modeling approach should be implemented whenever simulating DNAPLs under field-scale conditions.

5.2.3. Numerical Modeling Studies

The numerical modeling framework implements numerical methods such as the finite difference method, finite element method, and finite volume method to solve the flow and contaminant transport equation and can be applied to complex geometry and boundary conditions. The transport behavior of DNAPLs in an aquifer system has been the subject of numerous numerical modeling studies [12,13,54]. Key observations and details of models adopted in the numerical modeling-based studies are shown in Table S8 of the Supplementary Materials. In the early 2000s, a general-purpose reaction transport code (RT3D) was utilized to investigate the combined effect of rate-limited sorption, DNAPL dissolution, and kinetic-limited biological reactions [7]. Chapman and Parker [4] investigated the persistence of DNAPL plumes in an aquitard–aquifer system using 2D finite element software, namely, HydroGeosphere. The significance of considering intra-borehole vertical flow in fully-screened wells in the contaminant transport modeling was assessed using the Multi-Node Well (MNW) package within the MODFLOW-MT3DMS framework [42]. MODFLOW and a 2D reactive transport model (RT3D v2.5) were applied to investigate the vertical distribution of MTBE/TBA plumes in the saturated zone under field-scale conditions [6]. A 2D multiphase model (TMVOC and Petrasim) was used to analyze the impact of varying groundwater flow velocities on TCE infiltration and spreading behavior in a fully saturated sand layer [14].

The study examined the suitability of numerical models, specifically HydroGeoSphere, MODFLOW/MT3DMS, and FEFLOW, for simulating the diffusion into and out of LPPM in a sandy aquifer system [47]. It was found that a finite-difference-based numerical model (MODFLOW/MT3DMS) with suitable mesh size and time step closely matched the results of finite-element-based models such as Hydro Geosphere and FEFLOW [47]. However, their study considered single lumped values of dispersivity. In recent years, several modified versions of MT3DMS and MODFLOW have been created. For example, diffusion into/out of LPPM was included in the 3D MT3DMS modeling environment via a localized 1D analytical model [65]. A 3D stochastic numerical model, utilizing a lithofacies-based approach, was employed to accurately simulate the transport behavior of DNAPL at the Macro-dispersion Experiment (MADE) site in MS, USA [92]. The hydraulically heterogeneous domain was created using a Markov-chain-based stochastic technique. Furthermore, a 3D model was used to assess the long-term effectiveness of the pump and treat systems for TCE-contaminated sites [54]. Site-specific infiltration parameters were integrated into a fate and transport model for risk assessment of benzene at a Korean site and found an influence of the depth of the saturated zone on contaminant concentrations at different observation locations [93]. In another study by Guo et al. [75], MT3DMS and RWheat were compared while simulating pump and treat operations for mass removal at a contaminated site. Both models simulated measured field data with acceptable accuracy, although RWheat was better than MT3DMS on the basis of computational time [75].

Migration of the trichloroethene (TCE) plume originating from DNAPL source dissolution for an aquifer-aquitard system was investigated using a 2D numerical model for various degradation scenarios along with back diffusion [13]. A modified version of the MT3DMS was applied to understand the interplay between back diffusion, sorption, and DNAPL dissolution [16]. A comprehensive review of these studies reveals that the transverse dispersion coefficient was assumed as a lumped value, while the fitting parameters were not estimated via coupling of the transport model with an inverse optimization algorithm [49,54]. It can be stated that there is a need for numerical modeling-based field-scale studies in which unknown flow and transport parameters, such as dispersion and mass transfer coefficient, are computed by coupling the contaminant transport model with an inverse optimization algorithm.

6. Research Gaps and Future Perspectives

There are some potential gaps and research aspects that should be considered in future studies of DNAPL contamination in heterogeneous aquifer systems. A summary of the identified gaps and possible future directions are shown in Figure 8. These gaps are identified based on the review and bibliometric analysis of studies emphasizing the transport of DNAPLs in heterogeneous porous systems.

• Model calibration at the field-scale enhanced significantly when the effects of the dynamic flow model boundary conditions and intra-borehole flow were considered in the numerical models [42]. In addition, linear and exponential source depletion models effectively mimicked the back diffusion and its associated risks for gradual and noninstantaneous source depletion scenarios [84]. Thus, release models along with dynamic boundary and initial conditions should be involved in the mathematical modeling of DNAPL transport, especially under field-scale conditions;
• When modeling contaminant transport dynamics in highly heterogeneous aquifers, raw geological data with expert knowledge and an appropriate geostatistical technique can be taken into consideration;
• Microscale flow and transport processes through stagnant zones and low-permeability regions should be emphasized in greater amounts in the modeling framework when simulating the anomalous transport behavior of DNAPL under field-scale conditions;
• The space- and/or time-dependency of model parameters can be taken into account in 2D and 3D modeling frameworks for the simulation of DNAPLs under field-scale conditions;
• Spatial and temporal dependence of output metrics towards model input parameters can be quantified by applying spatial- and time-varying global sensitivity analyses.

7. Conclusions and Recommendations

The effective predictions of human health and ecological risks due to the presence of DNAPLs in subsurface systems are dependent upon contaminant transport model-based predictions. The prediction efficiency depends upon the concentration and mass fluxes at environmentally sensitive locations, which ultimately depends on a robust mathematical model. Thus, this article outlines a bibliometric survey and provides a concise review of DNAPL transport in porous systems to understand the past and present developments of modeling frameworks. Furthermore, this article provides an essential understanding of the factors and mechanisms governing dissolved-phase plume evolution of DNAPLs in heterogeneous porous systems. According to the bibliometric investigation of the top affiliations, the USA has seven affiliations, while Canada has three, indicating that these two nations have led the DNAPL transport research field over the past three decades.

Several studies related to DNAPL transport have been conducted to date; however, few articles have addressed the impact of low-permeability regions on transport dynamics. Even invariable functions of dispersion, mass transfer, and degradation coefficients were considered in modeling of DNAPLs at field-scale conditions. Thus, in the future, the variant flow and transport parameters dependent upon DNAPL concentration and/or distance and/or time can be considered in the modeling framework to understand their complex transport behavior in the heterogeneous porous system. The back diffusion and sorption mechanisms in low-permeability regions need to be studied extensively in DNAPL transport-related studies to elucidate the effects of nonlinear biogeochemical processes on the concentration and its associated risk metrics for the assessment of contaminated sites. Uncertainty and sensitivity analysis can be applied to analyze the effect of various model parameters on output indices, and thereafter a sensitivity-driven calibration approach can be implemented to determine the fitting parameters. The outcomes of the present study can provide a wide range of references, emerging topics, a list of dominant journals, and future perspectives for researchers to emphasize less focused on topics in the DNAPL transport research field.