A Screening Model to Predict Entrapped LNAPL Depletion

: Accidental leakage of hydrocarbons is a common subsurface contamination scenario. Once released, the hydrocarbons migrate until they reach the vicinity of the uppermost portion of the saturated zone, where it accumulates. Whenever the amplitude of the water table fluctuation is high, the light non ‐ aqueous phase liquid (LNAPL) may be completely entrapped in the saturated zone. The entrapped LNAPL, comprised of multicomponent products (e.g., gasoline, jet fuel, diesel), is responsible for the release of benzene, toluene, ethylbenzene, and xylenes (BTEX) into the water, thus generating the dissolved phase plumes of these compounds. In order to estimate the time required for source ‐ zone depletion, we developed an algorithm that calculates the mass loss of BTEX compounds in LNAPL over time. The simulations performed with our algorithm provided results akin to those observed in the field and demonstrated that the depletion rate will be more pronounced in regions with high LNAPL saturation. Further, the LNAPL depletion rate is mostly controlled by flow rate and is less sensible to the biodegradation rate in the aqueous phase.


Introduction
Accidental leakage of petroleum fuels from storage tanks or transfer pipelines is one of the main sources of subsurface contamination. The leaking hydrocarbons migrate within the unsaturated zone until they reach the vicinity of the saturated zone, where they accumulate as light non-aqueous phase liquid (LNAPL). When the LNAPL remains mostly as an immiscible fluid in the capillary fringe, the volatilization in the air/LNAPL interface represents the prevailing mechanism of LNAPL depletion with respect to benzene, ethylbenzene, toluene, and xylenes (BTEX) [1][2][3][4][5]; however, as the water level fluctuates in response to the cyclical alternation of the dry and rainy seasons, there is a continuous redistribution of the LNAPL because it migrates to lower aquifers as water levels decrease [6][7][8][9][10]. When the water level rises, LNAPL is retained by capillary force in the saturated zone, a phenomenon known as entrapment [4][5][6][7][8][9][10], which has been demonstrated in numerous previous studies (e.g., [11][12][13][14][15][16]). When entrapment occurs, LNAPL is entrapped by the capillarity in the saturated zone; in contrast, when the water level drops, LNAPL is released, and previously isolated ganglia clump together to gain mobility. Once entrapped, LNAPL releases soluble compounds into the groundwater, particularly BTEX, which represents the most soluble compounds [15,[17][18][19].
In tropical and subtropical climates, the amplitude of the water table fluctuation is typically high [15,20]. Under these conditions, the LNAPL may be completely entrapped below the water table, as described by Isler et al. [21]. The release of soluble compounds into the water from LNAPL comprises a multicomponent mixture (e.g., gasoline, jet fuel, diesel) and occurs because of the limited rate mass transfer phenomenon (non-equilibrium condition). In the case considered herein, this phenomenon is categorized as being diffusive. The total mass that is transferred from the LNAPL to the groundwater is calculated by multiplying the concentration gradient between the water and the effective compound solubility by the mass transfer coefficient and the phase-specific interfacial area, as demonstrated by Miller et al. [19]: where J = the mass flow of residual LNAPL into the water (ML −2 T −1 ), = the mass transfer coefficient (LT −1 ), = the maximum concentration of aqueous soluble compound (ML −3 ), C = the concentration of the water-soluble compound (ML −3 ), and = the specific interfacial area between the water and LNAPL (L 2 ).
Mass transfer will occur until the concentration equilibrium between the water and the entrapped LNAPL is reached. Because LNAPL represents a complex mixture of organic compounds with different solubility values, its dissolution becomes more complicated than that of a single compound. Further, the concentration of individual organic compounds is in equilibrium with the water that is in contact with LNAPL in the pore space. The equilibrium concentration depends on the composition of the mixture, as demonstrated by Banerjee [22] and Mackay et al. [23]. In water, the solubility of each compound is lower than the solubility of pure substances and can be determined in a porous medium with two immiscible liquids using Raoult's law: where , = the concentration of compound i in water (mol/L), , = the molar fraction of compound i in the mixture, , = the activity coefficient of compound i in the mixture, and = the maximum solubility of compound i in water (mol/L).
Because of the interphase mass transfer, there are strong variations in the concentration in the source area over time, wherein a decreasing trend is normally observed. The concentration decline is related to several factors, such as the LNAPL saturation, flow velocity, and biodegradation rates in the aqueous phase. Teramoto and Chang [15], Mackay et al. [23], Huntley and Beckett [24], and Thornton et al. [25], among others, demonstrated that the continuous loss of water-soluble compounds leads to continuous depletion of LNAPL in the source zone. Actually, there are a large number of analytical and numerical approaches, varying in scale and complexity, to simulate NAPL dissolution into the aqueous phase [18,[26][27][28][29][30][31][32][33][34][35][36][37][38][39]. Most of these approaches are based on empirical Sherwood-Gilland models that do not consider explicitly the NAPL/water interface area and were exclusively validated by lab-scale experiments. Based on the assumption that entrapped LNAPL continually loses mass to water, this study aims to present an alternative and simplistic numerical strategy for the long-term prediction of LNAPL depletion based on a phenomenological model to be used in practical purposes.

Numerical Model to Simulate the LNAPL Depletion
In order to predict the temporal change in the BTEX concentration in both the aqueous and nonaqueous phases within the source zone, in which the LNAPL is mostly entrapped, we developed a numerical model within the representative elementary volume (REV) that represents a discrete portion of the LNAPL source zone (Figure 1). The temporal change in the BTEX concentration is driven by interphase mass transfer from the LNAPL to groundwater and biodegradation in the aqueous phase.

Prediction of the BTEX Concentration in Groundwater
To predict the entrapped LNAPL depletion using the proposed model, we employed the algorithm proposed by Teramoto and Chang [15]. The model proposed by these authors, described by Equations (3)-(7), was developed to simulate LNAPL depletion by modifying the model proposed by Liu et al. [35] to include the first-order decay term to simulate biodegradation: , , , , , where , = the concentration of compound i in the aqueous phase, = the concentration, , , = the equilibrium concentration of compound i, , = the mass of compound i (M), = the coefficient of activity of compound i, = the water volume in a porous medium, = the mass transfer , , and , are marked by the subscript t and t − 1 to indicate the current and previous time steps, respectively.
The nonlinear system of Equations (3)-(6) is numerically solved using the iterative Gauss-Seidel method with central time weighting until convergence is reached. The model of Teramoto and Chang [15] was developed to reproduce mass transfer within an REV of the aquifer ( Figure 1). According to the formulation of Teramoto and Chang [15], the fluctuations in the water level may promote an increase or decrease in concentrations in the aqueous phase depending on the distribution of LNAPL within the smear zone. By using Equations (3)-(6), we understand that the concentration gradient is maintained by the Qin influx. The mass removal is mediated by the outflow of REV (Qout) and the firstorder decay kinetics λ.
The effective solubility (Equation (3)) is calculated by multiplying the activity coefficient ( ) by the molar fraction of the compound of interest in LNAPL and the solubility of the pure compound in water. The value of , according to several previous researchers [17,[40][41][42], is equal to or close to unity for aromatic hydrocarbons.
The main assumption of the proposed model is the complete entrapment of LNAPL below the water table (Figure 1). Whenever the mobile and residual phases in the unsaturated zone comprise an expressive fraction of LNAPL, the proposed model does not produce reliable results.

Prediction of the BTEX Concentration in LNAPL
The volume of water (Vw) in the REV is calculated as the saturated volume of the VER multiplied by porosity, discounting the proportion of pores filled by the LNAPL, according to Equation (8): where h represents the height of the saturated zone within the VER, n is the effective porosity, and Snw is the LNAPL saturation in the porous medium. The term , in Equation (3) represents the sum of all the constituents in the mixture (j = 1 to m) corresponding to the mass of the entrapped LNAPL in the porous medium. The volume of LNAPL in VER can be obtained from Equation (9): where and represent, respectively, the molecular weight and density of LNAPL. The concentration of the remaining monoaromatic compounds in LNAPL at each time step t can be obtained from Equation (10) as follows: where _ , represents the concentration of compound i in LNAPL at time step t.

Case Study
In order to demonstrate the reliability of the proposed algorithm, we applied it to a jet fuelcontaminated site located in the municipality of Paulínia, São Paulo, Brazil that has been previously described by other researchers [15,[42][43][44][45][46][47][48]. The volume of jet fuel is present in the subsurface with an estimated volume of 520 m 3 .
The studied shallow aquifer is composed of a heterogeneous Cenozoic unit, comprising medium to coarse sand, clayey sands, sandy clays, and clayey silts, deposited in a meandering river environment [42]. The geological heterogeneity explains the high variability of the hydraulic conductivity, ranging from 1.2 × 10 −7 m/s to 2.4 × 10 −4 m/s. The average hydraulic gradient is 0.0036 in the northeast region and increases to the southwest, toward the discharge zone, reaching values of 0.0176.
Teramoto and Chang [46] demonstrated that the water table fluctuations in the aquifer studied have high amplitudes, favoring the entrapment of LNAPL. Isler et al. [21] verified that most of the LNAPL is distributed as an entrapped phase nearly 2 m below the water table via the use of laserinduced fluorescence (Figure 2), as a result of strong water table fluctuations recorded in the field. As demonstrated by Teramoto and Chang [15], the concentrations of BTEX in the aqueous phase are directly related to the concentrations in the LNAPL. Teramoto and Chang [15] demonstrated that the decline of the BTEX concentration occurs slowly in the case of high LNAPL saturation (the center of the LNAPL source zone) within the pore system and occurs rapidly in the case of low LNAPL saturation (borders of the LNAPL source zone). Thus, the BTEX concentration in the LNAPL will reduce slowly in the center of the LNAPL source zone. Figure 3; Figure 4 exhibit the ethylbenzene and toluene concentrations, respectively, in the LNAPL. These figures demonstrate that the high concentration is observed in the center of the LNAPL source zone and declines towards the borders.

Results and Discussion
Under the entrapment conditions, mass transfer to water will be the most important mechanism for LNAPL depletion. Thus, to obtain more accurate estimates of LNAPL mass removal rates, the monitoring of BTEX compound concentrations in the aqueous phase must be emphasized. Based on field data, it is possible to make predictions about the long-term depletion of the source area and to make estimates of the time required for rehabilitation of the contaminated area. Appropriate strategies are required for the correct management of contaminated areas, quantification of the risk to human health, and assessment of the need for remediation.
The tests conducted using the algorithm we have developed will allow us to improve our understanding of the factors that control LNAPL depletion in the saturated zone. In the tests undertaken, the predominant factors for the drop in concentrations in the aqueous phase are controlled by several variables, such as the flow velocity, LNAPL saturation, LNAPL/water interfacial area, and biodegradation rate.

Prediction of LNAPL Depletion
The molar fraction of the BTEX compounds was obtained from the concentration of these compounds in fresh jet fuel. From the molar fraction of these compounds in jet fuel, the effective solubility of these compounds in LNAPL was calculated. The concentrations of the BTEX compounds in jet fuel, as well as the effective solubility values, are presented in Table 1.

Simulation of the BTEX Concentration Drop in LNAPL
The decline of the BTEX concentration in groundwater is directly related to the reduction in the molar fraction of these compounds in the entrapped LNAPL. Several parameters govern the rate of LNAPL depletion, mainly the flow rate, biodegradation kinetics, and LNAPL saturation. In order to evaluate the relative importance of these parameters, we performed several simulations, changing the input values of the model within the representative range commonly observed in the field. Table 2 presents the parameter values used in the algorithm we developed to simulate variations in the benzene concentration in an area contaminated by entrapped jet fuel, based on the initial molar fraction presented in Table 1. Because only monoaromatic compounds have effective solubility, and their mixing ratio is greatly reduced in jet fuel, it is possible to approximate the LNAPL mass ( , ) as being constant throughout the simulation.  Figure 4 shows the simulation results of the benzene concentration drop in both the aqueous and non-aqueous phases over 20 years. The effective depletion of benzene is associated with its reduced molar fraction within aviation jet fuel (7.36 × 10 −5 ); this depletion is associated with its high values of water solubility (1786 mg/L). Because of these characteristics, benzene concentrations tend to fall exponentially over time, describing first-order decay, which is the typical behavior of a Fickian phenomenon. By contrast, given its more expressive molar fraction (2.32 × 10 −3 ) and lower solubility (187 mg/L), ethylbenzene will be depleted less efficiently. The ratio of the removed BTEX relative to the remaining BTEX is significantly higher for a low saturation in the simulated example, resulting in a long-tail decreasing trend with respect to the BTEX concentration in the LNAPL ( Figure 5).

Depletion with a Change in the LNAPL Saturation
In order to evaluate the impact of Snw on the LNAPL depletion, we simulated three distinct scenarios, varying the Snw values. To reproduce the LNAPL depletion in the center of the LNAPL source zone, we assigned Snw the values of 0.22, 0.38 and 0.53. The model parameterization employed in the performed simulations is presented in Table 3. It is noted in Figure 6 that depletion will occur more intensely in regions with low LNAPL saturation and less intensely in regions with high LNAPL saturation. This result is in accordance with the results presented by Teramoto and Chang [15] and Ebehardt et al. [17]. These authors have demonstrated that a high LNAPL saturation can sustain high concentrations over a time scale that could span decades, while a low saturation would completely deplete the LNAPL within a few years. The difficulty in obtaining reliable values of may hamper the accurate prediction of LNAPL depletion using our approach; however, reliable estimates of this parameter may be derived from the floating phase thickness in monitoring wells, as demonstrated by Jeong and Chaberneau [8] and Lenhard et al. [12]. Alternatively, Gatsios et al. [16] and Teramoto et al. [49] demonstrated that may be estimated by the laser-induced fluorescence method via empirical models.
The specific interfacial area (A), a parameter intrinsically related to non-aqueous phase saturation, represents another parameter that is very challenging to estimate in the field. The parameter A may be estimated using , as previously demonstrated by Helder and Celia [50], Culligan et al. [51], and Brusseau et al. [52]; however, for practical purposes, this parameter may be optimized manually by trial-and-error during model calibration.

Depletion with a Change in Flow Rate
In addition to the LNAPL saturation, our analysis revealed that the parameter that controls the LNAPL source-zone depletion is the groundwater flow rate (Qin). To illustrate the impact of the groundwater flow rates on LNAPL depletion, we performed simulations varying the flow rate (Qin), while the remaining parameters were maintained at a constant value in three distinct scenarios. The Q or flow rates employed were 0.0122, 0.0243 and 0.0486 m/day, and the obtained results are presented in Figure 7. The results indicate that the LNAPL depletion will be more effective in aquifers with higher flow velocities (Q) and less intensive in aquifers with lower flow rates. As illustrated in Figure 4, increments in the flow rates promote sharp variations in the LNAPL depletion rate, prompting the concentration in the aqueous phase to also fall. The increase in the dissolution rate of LNAPL related to the increase in the flow rate is supported by several experimental studies, such as those of Seagren et al. [28], Saba et al. [18], Geller and Hunt [53], and Powers et al. [54].

Depletion with the Change in Biodegradation Rate
The algorithm proposed herein uses first-order decay kinetics to simulate hydrocarbon biodegradation in the aqueous phase. Because the metabolic action of native bacteria promotes the removal of BTEX from water, the concentration gradient between the aqueous phase and the effective solubility is increased, inducing an increase in the mass transferred. To evaluate the impact of biodegradation rates on LNAPL depletion in relation to benzene, three distinct scenarios were simulated. In the first scenario the biodegradation was not active, in the second scenario the λ was 0.0693 days −1 (10-day half-life), and in the third scenario the lambda was 0.0139 days −1 (half-life of 50 days). Figure 8 illustrates the benzene concentration decrease in the aqueous phase for the three simulated scenarios.
By comparing the various results obtained from simulations with varying Snw (Figure 6), Q (Figure 7), and λ ( Figure 8) values, we can conclude that the LNAPL depletion is less sensitive to the biodegradation rate of BTEX in the aqueous phase. On the other hand, the LNAPL depletion is highly sensitive to variations in the flow rate and to variations in Snw in particular; thus, the Snw and Q are key parameters to be obtained in the field to produce reliable predictions.

Applications in Real Contamination Cases
Because the BTEX concentration measurements in LNAPL are far less common than those measured in groundwater, the model calibration should be performed using a time-series of the BTEX concentration in groundwater. The prediction of the LNAPL source-zone depletion in real cases using data on temporal changes in the BTEX concentration in the aqueous phase was previously demonstrated by Teramoto and Chang [15]; however, it is necessary to emphasize that, because of the large number of parameters employed in the model, many of which are obtained from the model's calibration, there are problems related to the non-uniqueness of the performed simulations. In the case of a non-singular model, the different sets of values employed may offer similar results. To minimize the uncertainties associated with non-singularities, it is recommended that two or more compounds be adjusted simultaneously, because the set of REV and LNAPL parameters (Q, Vw, A, and Snw) should be the same for both of these compounds. Figure 9 illustrates the simultaneous adjustment for total xylenes and ethylbenzene in the aqueous phase. The presence of subsurface hydrocarbons represents one of the most common problems around the world. Given the risk to human health from many LNAPL constituents, particularly monoaromatic compounds, severe penalties for the mass removal or destruction of remediation systems are usually put in place by environmental agencies. Many studies have shown, however, that the processes associated with natural source zone depletion are efficient in the mass destruction of organic contaminants and in long-term risk reduction. In the scenario described, no intervention is required for mass destruction.
The proposed model represents an upscaled mass transfer approach that explicitly assumes the rate-limiting dissolution process. Seagren et al. [28] and Lekmine et al. [29] experimentally demonstrate that the local equilibrium assumption best reproduces the mass transfer from LNAPL to groundwater. On the other hand, Mobile et al. [27], Nambi and Powers [30], and Zhu and Sykes [34] indicate that the local equilibrium is unsuitable, and the rate-limiting assumption offers the most realistic simulations. Thus, despite the ability of our work to reproduce the LNAPL dissolution, the most suitable approach to reproduce the LNAPL dissolution remains imperfectly answered.

Conclusions
In tropical and subtropical environments, the high amplitude of water table fluctuations favors LNAPL entrapment. Based on the assumption that LNAPL depletion in the entrapped phase of the saturated zone occurs predominantly by solubilization in groundwater, a numerical algorithm was developed that reproduces the continuous drop of the BTEX molar fraction in LNAPL for more realistic predictions. The results point to the potential of the methodology proposed in this work to estimate the natural depletion rates of the source area and to predict the time required for the rehabilitation of the contaminated area. The algorithm proposed herein can thus be employed to assist in the management of contaminated areas and support decision-making.