Temporal moment-based approach to understand the dissolved-phase LNAPL recovery and associated characteristics in the porous system under dynamic groundwater table conditions

Pollution of soil–water resources due to the release of hydrocarbon contaminants is a serious threat to human and ecological health (Garg et al. 2015, 2017). Light non-aqueous phase liquids (LNAPLs) have major concerns in the subsurface due to their high sensitivity to varying subsurface conditions, which enable them to spread widely (Kim & Corapcioglu 2003; Oostrom et al. 2006). The LNAPLs are not freely miscible in water, but a small quantity of dissolved phase in the soil–water system is sufficient to seriously degrade soil–water quality (Essaid et al. 2015; Gupta et al. 2019; Gupta & Yadav 2020). Typically, the main sources of these dissolved hydrocarbons are residual mass, either lighter or denser-NAPLs, which may have been spilled in the past and generally occur near fuel production sites and petrol stations (U.S. EPA 2006). The hydrocarbon compounds may start dissolving from residual NAPL with soil moisture in the capillary zone and create a mobile plume, referred to as ‘plume’ in the manuscript, which subsequently moves to the downgradient locations due to advective flow (Chrysikopoulos 1995; Narayanan et al. 1998a, 1998b; Fagerlund et al. 2008; Yadav et al. 2012). In laboratory-scale conditions, the flow pattern changed from linear to non-linear, even with a gradual increase in the hydraulic gradient (Basack et al. 2022). Thus, it is important to understand the fate of plumes thoroughly under varying groundwater flow regimes in order to frame a suitable engineered remediation strategy. A series of two-dimensional (2D) sand tank experiments were performed by Gupta et al. (2019) to understand the fate of plumes under four groundwater table conditions. They observed a large plume with varying concentrations (15–200 mg/L) of toluene, a representative hydrocarbon, under rapid groundwater table fluctuation (GWTF), referred to as a rapid case in this manuscript. Their results are also in-line with Gupta & Yadav (2020), wherein they highlighted the dependency of biodegradation on the initial substrate concentrations. They observed degradation increases with an increase in initial concentration of toluene up to 50 ppm and remained maximum till 100 ppm before starting to decrease with an increment in substrate concentration. A low degradation rate at a high substrate concentration shows the toxic effects of pollutants on toluene degraders (Gupta et al. 2020). In order to improve treatability, it is necessary to estimate the recovery of dissolved hydrocarbon mass from the high-concentration zone first.

Several mathematical modeling techniques, such as the finite-difference method (FDM) and finite-element method (FEM), were applied to understand the groundwater flow (Beegum et al. 2018; Omar et al. 2019, 2020, 2021) and contaminant transport in the subsurface system (Šimůnek et al. 2016; Guo et al. 2019; Ekeleme et al. 2021). For example, the analytic element method (AEM) and FDM were used to study the groundwater flow dynamics in the Lower Ganga river basin from 2004 to 2017 (Omar et al. 2019). Their study found that the AEM performed better than the FDM for steady-state solutions. However, the FDM was found to be better than AEM in the case of transient modeling and heterogeneous aquifer system. Similarly, the AEM model was integrated with particle swarm optimization (PSO) and compared with the FDM-PSO model in the identification of flow parameters of unknown wells (Gaur et al. 2023). Their study found the AEM-PSO model as more effective than the FDM-PSO model. In another study, a loosely coupled SWAT–MODFLOW model was used to study the effects of recharge and abstraction on the groundwater levels of the Nairobi aquifer system (Nyakundi et al. 2022). A hybrid SWAT–MODFLOW model was found to be an efficient and realistic approach as it incorporated factors such as land use-land cover, topography, climatic factors, and soil type (Nyakundi et al. 2022). Furthermore, moment analysis based on experimental and model-simulated concentration data was found to be an effective approach to studying the flow and contaminant transport behavior in the subsurface system (Govindaraju & Das 2007; Blackmore et al. 2018; Guleria et al. 2020; Singh et al. 2021, 2022).

Temporal moment analysis (TMA) proved to be an effective method to understand the time-averaged response of contaminants like hydrocarbons, in the subsurface system as compared to the conventional way of analyzing breakthrough curves (BTCs) at downgradient location(s) or spatial concentration profile(s) (Pang et al. 2003; Govindaraju & Das 2007; Sharma et al. 2012; Goltz & Huang 2017). For example, flow and solute transport behavior at site-specific and laboratory-based column conditions was studied based on the temporal moments (Blackmore et al. 2014). Also, the relationship between the dispersion coefficient and seepage velocity for the transport of carbon disulfide vapor in the partially saturated 2 m long column was derived based on temporal moments (Kleinknecht et al. 2017). The experimental and theoretical data of tracer transport in soil columns under single- and two-phase flow conditions was analyzed based on temporal moments (Singh et al. 2021, 2022). On the basis of the temporal moments as output metrics, a comparison of the conventional model with the hybrid mathematical model to analyze the simulation capabilities was made (Guleria et al. 2020; Sharma et al. 2020). Thus, utilizing observed and physics-based simulated concentration data under various groundwater flow regimes via TMA can be a simple technique for estimating the dissolved hydrocarbon recovery and its associated parameters from plumes.

The experimental data from Gupta et al. (2019) were rare and unique, as the fate and transport of LNAPL were simulated under four groundwater table conditions. A question yet to answer is how much-dissolved hydrocarbon can be recovered from such a dynamic subsurface system. Furthermore, the role of groundwater table dynamics on the dissolved hydrocarbon recovery was unclear even after Gupta et al. (2019). Gupta et al. (2019) analyzed the transport dynamics of the dissolved-phase plume emanating from the LNAPL pool accumulated in a capillary fringe under stable and fluctuating groundwater table conditions based on the BTCs and spatial concentration profiles. In this research, it is estimated how much-dissolved hydrocarbon could be recovered from such a dynamic subsurface system which can help field practitioners and pollution remediation managers. It provides details on the role of groundwater table dynamics on dissolved hydrocarbon recovery and its associated parameters.

For modeling of toluene in the sand tank experimental setup conditions, the Galerkin FEM was used to solve the governing flow and transport equations in HYDRUS-2D software (Šimůnek et al. 2006). A 2D numerical domain with dimensions corresponding to those of the sand tank was created for the simulation of soil–water flow and LNAPL transport through the subsurface system. To solve the water flow equation, saturated moisture content was assumed as the initial moisture level, and the top side of the domain was taken as WT boundary condition. No-flux condition was taken for a bottom boundary, while a continuous flux-type boundary with pulse condition was considered at the left and right side of the domain for GWTF cases. In the GWTF-Stable case, influx and outflux boundaries without pulse conditions were assumed. Furthermore, to solve the LNAPL transport equation, zero toluene concentration was assumed as an initial condition for the whole numerical domain. An LNAPL source point was considered at the top-left side of the domain, mimicking the 2D sand tank conditions. To comprehend the LNAPL transport dynamics under various GWTF scenarios, simulations were run for 56 h.

The temporal moment of dissolved-phase contaminant concentrations can provide the time-averaged response of solute within the porous media (Naff 1992; Guleria et al. 2020). The first four moments (out of n possible moments), were commonly used to define the central distribution tendency (Govindaraju & Das 2007). For example, the zeroth moment relates to total mass density, and the first moment relates to the arithmetic mean (Pang et al. 2003; Guleria et al. 2019). In physical meaning, zeroth temporal moment (ZTM) reflects the contaminant mass recovered at a particular location or observation point. Furthermore, the first temporal moment describes the time when half of the injected mass has passed through the observation point or the average time spent by contaminant molecules in the porous system. The second temporal moment represents the variance of the BTC, which provides a measure of the degree of mixing of the mass within the porous system (Govindaraju & Das 2007; Guleria et al. 2020). Thus, temporal moments of toluene concentrations (mass recovered, mean residence time) were computed by calculating the absolute temporal moments at a location before determining the ratio of various absolute temporal moments (Figure 1). The zeroth, first, and second normalized temporal moments were computed using numerical integration in the MATLAB environment. After computing the various temporal moments, correlation analysis between observed and simulated moments was carried out to determine the discrepancies.

It was observed that the increasing rate of WT fluctuation leads to an increase in the dissolved toluene mass recovery and can be helpful in the enhancement of advection processes with increasing hydraulic gradient, as shown in Figure 2. Also, in the GWTF-Rapid case, a high rate of WT and pore velocity fluctuations might have advected the solute mass entrapped in the microscopic immobile regions due to which a large amount of mass seeps out of these regions, finally resulting in the highest solute mass recovery as compared to the GWTF-Stable case. The magnitude of toluene mass recovered was found to decrease at a mild slope for the GWTF-Stable case, which can be seen in Figure 2(a) from the isoline intervals. However, for the GWTF-Rapid case, toluene mass recovered was found to decrease at a high gradient while moving away from the LNAPL source location, as depicted in Figure 2(b). It is hypothesized that the location or observation point with highest solute mass recovered can be treated as a hot spot compared to other observation points, which need to be emphasized during remediation operations.

The largest value of second temporal moment was observed for the GWTF-Rapid case in the top layer, depicting the strong variation of hydraulic conductivity field, which can be attributed to non-linear dependence of hydraulic conductivity and dispersion coefficient on moisture content level in the top layer. However, the lowest value of second temporal moment was observed for the GWTF-Stable case in the saturated region, which can be caused by anoxic conditions in the deeper saturated zone. Shorter mean residence time and larger variance of BTCs were observed in Figures 3 and 4 for a particular location, signifying the fast flow paths and dominance of advective transport. These observations are in-line with the findings of Blackmore et al. (2018) for contaminant transport experiments in waste rock piles. These similarities fortify the existence of microscale heterogeneity in the tank setup, though care was taken while packing the sand in the tank to create homogeneous conditions. It is also important to note that the mixing of end products via biodegradation of toluene may cause such variation in BTC, as this was a case of reactive transport.

Studies have reported the direct influence of WT fluctuations on the recovery rate of LNAPL (Wang et al. 2014; Teramoto et al. 2020). In the present study, the highest toluene mass recovery was observed under rapidly fluctuating WT conditions compared to slow (GWTF-Slow) and general (GWTF-General) cases. Similar findings of increased concentration of LNAPL dissolved-phase components due to WT fluctuations were reported for lab-scale (80 cm × 51 cm × 3.3 cm) aquifer conditions (Dobson et al. 2007). WT fluctuations caused an increase in biodegradation activity along with a dissolution of LNAPL components as compared to the model without fluctuation (Dobson et al. 2007). Not only the magnitude but the rate of WT fluctuation was observed as a governing factor controlling the spatiotemporal variation of LNAPL in the porous media (Oostrom et al. 2006). Similarly, in the present study, variation in the mass recovery and mean residence time was observed with several WT fluctuation scenarios, viz. GWTF-Rapid, GWTF-Slow.

In the study by Gupta et al. (2020), an increase in the microbial growth of hydrocarbon degraders was observed for fluctuating WT conditions in the column setup in comparison to static WT conditions. It is revealed that the WT fluctuation can act as a powerful tool for the remediation of NAPL-contaminated sites (Gupta et al. 2020). Similarly, the enhanced toluene mass recovery along with a decrease in the mean residence time was observed in the present study for the GWTF-Rapid case. Enhancement of the toluene mass recovery via manipulation of groundwater may help in reducing remediation costs and time. Also, an increase in the initial substrate concentration due to WT fluctuation will eventually help the remediation of contaminated sites governed by microbes. The quantification of the above findings via various temporal moments highlighted the significance of temporal moment-based approach and, thus, can be helpful in designing remediation strategies.

Hydrological interventions in the form of a dynamic groundwater act as a strong LNAPL removal tool. In this study, the ability of groundwater table manipulation has been tested to improve dissolved toluene mass from a plume zone. The study strongly supports the applicability of temporal moments to predict the time-averaged response of toluene plumes emanating from LNAPL pool in a soil–water system, which can be linked with management and remediation practices. The approach presented in the study would facilitate a monitoring campaign to determine sampling locations, operational parameters, etc., that can further support planning on the remediation strategies of various LNAPL-contaminated sites. However, there are a few limitations that could cause uncertainties in estimating temporal moments. The sensitivity analysis of the flow and transport parameters using TMA was not implemented. In the future, the impact of low permeability regions such as silt/clay lenses on the plume evolution in the saturated porous media can be investigated using TMA. The implementation of TMA at field-scale conditions could be the next step to check the applicability of temporal moment-based approach. In conclusion, we strongly recommend that field managers manipulate the groundwater table dynamics by implementing pumping operations to improve the hydrocarbon mass recovery.

The authors are thankful to the Department of Science and Technology (DST), Government of India, for funding this research under the scheme of Ramanujan fellowship. The first author would like to acknowledge the Indian Institute of Technology, Delhi, for providing a Doctorate Fellowship. The authors are also thankful to the University Grant Commission, New Delhi, for providing a JRF/SRF for this study.

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by A.G. and P.K.G. The first draft of the manuscript was written by A.G. and all authors commented on previous versions of the manuscript. The resource and mentorship were provided by S.C. and B.K.Y. All authors read and approved the final manuscript.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.