Next Article in Journal
Data Modeling of Sewage Treatment Plant Based on Long Short-Term Memory with Multilayer Perceptron Network
Next Article in Special Issue
How to Minimize the Environmental Contamination Caused by Hydrocarbon Releases by Onshore Pipelines: The Key Role of a Three-Dimensional Three-Phase Fluid Flow Numerical Model
Previous Article in Journal
Tuning the Optical Properties of ZnO by Co and Gd Doping for Water Pollutant Elimination
Previous Article in Special Issue
Development of a Distributed Control System for the Hydrodynamic Processes of Aquifers, Taking into Account Stochastic Disturbing Factors
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Migration of DNAPL in Saturated Porous Media: Validation of High-Resolution Shock-Capturing Numerical Simulations through a Sandbox Experiment

1
Department of Chemistry, Life Sciences and Environmental Sustainability, University of Parma, 43124 Parma, Italy
2
Department of Engineering and Architecture, University of Parma, 43124 Parma, Italy
*
Author to whom correspondence should be addressed.
Water 2023, 15(8), 1471; https://doi.org/10.3390/w15081471
Submission received: 7 March 2023 / Revised: 4 April 2023 / Accepted: 6 April 2023 / Published: 10 April 2023

Abstract

:
This paper shows a comparison between experiments carried out in a laboratory-scale sandbox where the migration of a dense nonaqueous phase liquid (DNAPL), hydrofluoroether (HFE-7100), in a saturated porous medium was investigated, and validation was performed using high-resolution shock-capturing numerical simulations to resolve the nonlinear governing coupled partial differential equations of a three-phase immiscible fluid flow. The contaminant was released using a colored fluid as a tracer for a fixed time and pressures different from the atmospheric one into the saturated zone, first by using a column laboratory experiment, and then a sandbox-scale example with a hydraulic gradient. A digital image analysis procedure was used to determine the saturation distribution of the contaminant during its migration. These results are compared with the values determined for a DNAPL migration in a similar porous media through a numerical simulation. They show good agreement with the experimental results and also show that CactusHydro can follow the migration of a plume evolution very precisely and can also be used to evaluate the effects and environmental impacts deriving from leaks of DNAPL in saturated zones.

1. Introduction

Groundwater contamination due to the release of a dense non-aqueous phase liquid (DNAPL) characterized by a density higher than water that does not dissolve in or easily mix with water, such as oil, gasoline, and petroleum products, represents a significant environmental pollution problem due to the extensive use of these liquids in industrial or commercial processes. For these reasons, it is important to understand the migration fate of these contaminants in controlled laboratory conditions, such as porous medium in sandboxes with fixed boundary conditions and a type of dense contaminant.
When introduced into the environment, a DNAPL migrates downward, due to the gravitational effects, as a different liquid, first in the unsaturated zone, then in the saturated one. Besides this vertical movement, there is a lateral spreading due to capillary effects and hydraulic gradient. The distribution of the contaminant may vary depending on the medium characteristics, and entrapped DNAPL dissolves slowly into the groundwater flow and acts as a long-term source of contamination. Some examples have been proposed to investigate the migration of a DNAPL (e.g., [1,2,3,4,5,6,7]).
Two main measurement approaches for measuring contaminant dispersion in sandboxes aquifer are reported in the literature. One of them uses in-situ probes to measure the tracer concentration. The second one is based on an image analysis that allows for obtaining time lapse images of the tracer. Since this procedure is not invasive, it does not influence the fate of the flow and the concentration field [8,9]. In the hydraulic laboratory of the University of Parma, we set up an aquifer model sandbox (and a column) and a next-generation imaging system to investigate DNAPL migration using digital photography (a color camera with high resolution) to determine the concentration of the DNAPL in terms of reading and spatial details. In a previous work [5], a hydrofluoroether (HFE-7100), a colorless DNAPL, was used. Since the HFE-7100 is not miscible with water and the saturation concentration of fluorescein sodium salt is too low to be observed through imaging analysis, instead of colouring the HFE-7100, as was performed in a previous work [5], the experiments were carried out considering water with fluorescein sodium salt as the background fluid and HFE-7100 as the tracer. The contaminant can then be detected in the image analysis where the water is not colored.
The determination of the spilled DNAPL migration/spill in a porous medium (in particular, in a column/sandbox laboratory-scale model) can be described using numerical simulations of the governing equations of immiscible phase fluid flow in a porous medium. For example [10,11], laboratory studies applied imaging techniques [4], NAPL mass dissolution [4,5,12], and the relevance of heterogeneities [13,14]. In Citarella et al. [6], the main assumptions of the model are that the flow and transport phenomena are uncoupled and that the flow can be studied in a vertical plane; the porous medium was considered homogeneous and isotropic. The groundwater flow was reproduced using MODFLOW 2000 [15], and the transport process was reproduced by using MT3DMS with the total variation diminishing (TVD) method as the advection solver package [16].
The difference between the numerical approach used in [4] and this paper is that we used a numerical code, CactusHydro, recently introduced in [17,18], based on the high-resolution shock-capturing (HRSC) flux conservative method [19,20,21] that follows sharp discontinuities accurately and temporal dynamics of three-phase immiscible fluid flow in a porous medium; this shows the absence of spurious oscillations in the solution and converges to the ‘weak’ solution as the grid is refined. The time evolution was performed using a forward-in-time explicit (forward) method (instead of the most commonly used implicit ‘backward’ in time evolution method). CactusHydro is based on the Cactus computational toolkit [22,23,24], an open-source software framework for developing parallel high-performance computing (HPC) simulation codes, and the data are evolved on a cartesian mesh using Carpet [25,26]. The migration of the spilled DNAPL contaminant was considered immiscible. The effect of volatilization, biodegradation, or dissolution was not considered. The vertical and horizontal movement of the contaminant are coupled and are numerically resolved as a unique zone (and not separating the vertical movement from the horizontal one since the flow equations include both zones).
The purpose of this work is to validate the numerical simulations of a DNAPL leak in a saturated porous medium through the comparison between numerical and experimental results. In more detail, the HRSC flux conservative method and the CactusHydro code [16,17] were utilized to simulate the DNAPL leaking with a pressure different from the atmospheric one for a fixed time, and a sandbox experiment (after being calibrated, the system used a column experiment) was carried out to verify the effectiveness of the numerical simulations. At this stage, a comparison was made in order to verify the ability of this numerical approach to make reliable predictions about the DNAPL migration velocity, considering the whole DNAPL plume as a reference.

2. Materials and Methods

2.1. Governing Equations

The governing equations that describe three-phase fluid flow in a porous medium composed of nonaqueous ( n ), water ( w ), air ( a ), and a variably saturated zone were introduced in [16,17] and are given by:
t ρ n ϕ S n = x i ρ n k r n μ n k i j p a x j + ρ n g z x j x i ρ n k r n μ n k i j p c a n x j + q n ,
t ρ w ϕ S w = x i ρ w k r w μ w k i j p a x j + ρ w g z x j x i ρ w k r w μ w k i j p c a w x j + q w ,
t ρ a ϕ S a = x i ρ a k r a μ a k i j p a x j + ρ a g z x j + q a ,
where x i = x , y , z are the spatial cartesian coordinates; t is the time coordinate; ρ α   M L 3 , μ α   M L T , p α   M T 2 L , k r α , k i j   L 2 , g   L T 2 , z   L , q α   M L , and ϕ are the density for each phase, α = w , n , a , the dynamics viscosity of phase α , the pressure of phase α , the dimensionless relative permeability of phase α , the absolute permeability tensor, the gravity acceleration, the depth, the mass source/sink, and the porosity, respectively. A fourth equation accounts for S α , the dimensionless volumetric saturation of phase α , and satisfy the relation:
S w + S n + S a = 1 .
The system (1)–(4) is written in terms of the pressure p a (that is, the air pressure when S a is different from zero), and the saturations S w , S n , S a . It is also written in terms of the capillary pressure for the air–water phase, p c a w = p a p w , and the capillary pressure for the air–nonaqueous phase, p c a n = p a p n , where it is used, p w = p a p c a w , and p n = p a p c a n . The nonaqueous–water capillary pressure is given by p c n w = p n p w = p c a w p c a n (in contrast to Refs. [27,28], where the air gradient pressure is assumed to be negligible). Both the relative permeabilities k r w , k r n and k r a , and the capillary pressures are functions of the saturations, k r α = k r α S a , S n , S w , p c a n = p c a n S a , S n , S w and p c a w = p c a w S a , S n , S w .
For three-phases, they have been extended from the two-phase expressions [29] and can be found in [17,18]: k r w = S e w 1 / 2 1 1 S e w 1 / m m 2 , k r a = 1 S e t 1 / 2 ( 1 S e t 1 / m ) 2 m , and k r n = S e t S e w 1 / 2 1 S e w 1 / m m 1 S e t 1 / m m 2 , where the total effective liquid saturation, S e t , is defined in terms of the irreducible wetting phase saturation S w i r , i.e., S e t = S w + S n S w i r 1 S w i r . For the capillary pressure, instead, the van Genuchten model [30] is used where the effective saturation, S e , is given by: S e = 1 + α p c n 1 1 n , and α and n are model parameters. p c is the capillary pressure head: p c = p c 0 1 S e 1 / m 1 m , where m = 1 1 n , and p c 0 = α 1 is the capillary pressure at S e = 0 . Since p c a w = p c a n + p c n w , the capillary pressures are given by p c a n = p c a n 0 1 S e t 1 / m 1 m and p c a w = p c a n 0 1 S e t 1 / m 1 m p c n w 0 1 S e w 1 / m 1 m .
The strong nonlinearities in the partial differential equations PDEs system (1)–(4) are represented by the relationship between the permeability–saturation, which is responsible for creating sharp discontinuities (shocks), and the capillary pressure–saturation of each phase. Therefore, the dominant part of the multiphase PDEs flow that has to do with the water and nonaqueous phase is dominated by the hyperbolic part (the one proportional to gravity and the gradient of the pressure) rather than the one proportional to capillary pressures (elliptic part of the PDEs (1)–(4)). To numerically resolve the partial differential equations PDEs (1)–(4), we used the HRSC methods [19,20,21] introduced in [17,18], which treat the hyperbolic part of the PDEs and eliminate the oscillations/shocks and converge to the ‘weak’ solution of the system.

2.2. Experimental Equipment

2.2.1. Tracer

The lab experiment’s objective was to evaluate the fate of a DNAPL in porous media. As in a previous work [5], hydrofluoroether, HFE-7100 (Sigma-Aldrich, Saint Louis, MO, USA. Product number: SHH0002, Formula: C5H3OF9), a non-toxic, non-flammable, colorless DNAPL, was used because of its similar properties to TCE. Considering that HFE-7100 is not miscible with water and the saturation concentration of fluorescein sodium salt is too low to be observed through an imaging analysis, instead of coloring the HFE-7100, as was performed in a previous work [5], the experiments were carried out considering water with fluorescein sodium salt as a background fluid and HFE-7100 as a tracer. Therefore, the HFE-7100 can be detected in the image analysis where the water is not colored.

2.2.2. Experimental Setup for DNAPL Migration in a Column

The preliminary tests were performed in a column in a steady state (see Figure 1). The column, with an internal diameter of 50 mm and a length of 50 cm, was filled with glass beads of 1 mm diameter as the porous media and saturated with water and fluorescein. Then, 50 mL of HFE-7100 was injected for 3 s at the top of the column, and the DNAPL movement was observed through photographic equipment every 0.5 s for 400 s.

2.2.3. Experimental Setup for DNAPL Migration in a Sandbox

The experiments were performed in a sandbox built with transparent plates; see Citarella et al., 2015 [6] for the details of the setup. The external dimensions of the sandbox were 1.20 m × 0.73 m × 0.14 m. Along the longest axis, the sandbox was made up of three parts (Figure 2): two tanks (upstream and downstream), which allow the regulation of the water level and, consequently, of the flux, and a central chamber (length L = 0.954 m, height H = 0.70 m and thickness TH = 0.10 m), which contains the porous medium. Two weirs controlled the water level in the upstream and downstream tanks to maintain the same boundary conditions without oscillations during the experiments. The porous medium consisted of glass beads with a diameter between 0.75 and 1 mm and a density of 1480   kg / m 3 . The material was packed in order to avoid non-uniformity of the media. The porosity of the medium was estimated at 37% [31], and the bulk hydraulic conductivity K was estimated to be K = 6 × 10 3   m / s [6].
The experimental device was placed in a darkroom to avoid external light contamination. Considering that the experimental apparatus was developed for using fluorescein sodium salt as a tracer, the darkroom was lighted by means of 8 monochromatic blue LED lights [6]. Then, 4.89 mL/s for 45 s of HFE-7100 was injected at the red dot (see Figure 2) and the DNAPL movement was observed through photographic equipment every 5 s for approximately 215 s. The spill rate was selected so as to simulate the transport resulting from instantaneous contaminant injection in the subsurface, at the laboratory scale.

2.2.4. Data Acquisition

The evolution of the DNAPL spread was evaluated through image analysis. The digital images were acquired with a Canon EOS 40D camera with a 16–35 mm zoom lens, placed inside the darkroom on a fixed tripod at a distance of 1.6 m from the sandbox. The camera was fully controlled by a computer, and the images were shot with a spatial resolution of 3888 × 2592 pixels in RAW format with 14-bit color depth for each channel. Ad hoc software developed in the Labview® environment, version 11.0.0.4029 (National Instruments, 2011) stored the injection and background rate and controlled the injection timing and the camera.

3. Results

This section shows the DNAPL migration three-dimensional HRSC numerical simulations results performed in both a column experiment and then a sandbox laboratory experiment. For the case of the column experiment, the DNAPL was released into a saturated zone for a fixed time equal to 3 s. For the case of the sandbox experiment, instead, the DNAPL was released into the saturated zone for a fixed time of 45 s. We compared these numerical results with the experimental ones coming from the DNAPL leakage from a column and the experimental results of the NAPL leakage in a sandbox experiment explained in Section 2.2.
From the numerical point of view, it is considered a three-phase immiscible fluid model composed of water, air, and nonaqueous phase liquid (DNAPL) defined in Section 2.1 [17,18] and includes, as an initial condition, a DNAPL contaminant that is released for a fixed time with a pressure higher than the atmospheric one. The temporal evolution of the migration of the immiscible contaminant was investigated following the saturation contour profiles of each phase along the saturated zone and at different times.

3.1. Numerical Setup of a DNAPL Migration

3.1.1. Numerical Setup of a DNAPL Migration in a Column

We first consider a DNAPL leak, HFE-7100, of density ρ h f e 7100 = 1500   k g / m 3 [4], released for three seconds (at a spill rate of 0.0167   k g / s ) inside a saturated column zone, such as the one in Figure 1. The contaminant is released for 3 s with a pressure of 4118.0   P a . The saturated zone grid geometry is assumed to be a parallelepiped (see Figure 3) of 0.060   m long from x = 0.030 , + 0.030   m (left-hand side), 0.060   m wide from y = 0.030 , + 0.030   m (right-hand side), and 0.360   m depth, z = + 0.0 , 0.360 m , with a spatial resolution of d x = d y = d z = 0.01   m , and a time step size d t = 0.001   s . The DNAPL is released into the saturated zone on top of a parallelepiped grid placed at z = 0.040 , 0.035   m , x = 0.005 , + 0.005   m , and y = 0.005 , + 0.005   m , at t = 0   s , as shown in Figure 3. Initially, the saturation of the contaminant inside this box is one. The experimental porous medium is a saturated zone and is composed of spherical glass beads equivalent to sand and porosity equal to 0.37 in the numerical model. All the boundary conditions are no-flow except for the infiltration zone on top of the parallelepiped and an impermeable zone at the bottom. The legend at the right-hand side of Figure 3 indicates the saturation contour of the DNAPL, σ n = S n ϕ in color bars. The left-hand side of Figure 3 shows the z x plane, while the right-hand side shows the z y plane.
Table 1 gives the material properties and parameter details used in the numerical simulations performed with CactusHydro [17,18] together with the initial conditions shown in Figure 3. It shows the density of DNAPL HFE-7100, which is 1500   k g / m 3 . The porosity of the porous medium is fixed to be 0.37, which represents a sand geological structure [31]. The value of the hydraulic conductivity was measured to be K = 6.0 × 10 3   m / s . From this value, we calculated the absolute permeability, k = K μ w ρ w g , where we used the density and the dynamics viscosity values of the water, while g = 9.8   m / s 2 , and gives k = 6.122 × 10 10   m 2 (this is the value that has been used in the numerical model; see Table 1). We set this value for the entire grid with k x = k y = k z . The relative permeabilities and capillary pressure were obtained using the equations in Section 2.1. The van Genuchten α parameter is given by α = p c ρ w g 1 , where p c is the capillary pressure head equal to 676.55   P a .
Using the values of the surface tension and the interfacial tension of the different phases and, in particular, for the DNAPL HFE-7100 in Table 1 (see Ref. [4]), one can calculate the coefficients β a n = σ a w σ a n = 71.75 13.60 = 5.28 and β n w = σ a w σ n w = 71.75 35.59 = 2.02 , the capillary pressure nonaqueous-water, p c n w S w = p c a w β n w = 334.93   P a , and the capillary pressure air-nonaqueous, p c a n = p c a w p c n w = 341.62   P a .
The value of the dynamic viscosity for the DNAPL HFE-7100 was set to be 1.35 × 10 3   k g / m s [32]. The irreducible wetting phase saturation S w i r = 0.045 (in Table 1) was taken from Refs. [33,34] corresponding to a porous medium ‘Sand’.

3.1.2. Numerical Setup of a DNAPL Migration in a Sandbox

We now consider the same contaminant as before, the DNAPL leak, HFE-7100, whose density is ρ h f e 7100 = 1500   k g / m 3 , released for forty-five seconds (at a spill rate of 4.89 × 10 3   kg / s ) inside a saturated sandbox laboratory experiment (similar to Figure 2). The contaminant is released with a pressure of 3118.0   P a and a saturation of one. The direction of the release is in the direction of the x axis (left and right). The saturated zone grid geometry is assumed to be a parallelepiped of 1.16   m long from x = 0.88 , + 0.28   m (see Figure 4), 0.10 m wide from y = 0.05 , + 0.05   m , and 0.80   m depth, z = + 0.15 , 0.65   m , with a spatial resolution of d x = d y = d z = 0.01   m , and a time step d t = 0.001   s . The DNAPL is released in the saturated zone on top of a parallelepiped grid placed at z = 0.180 , 0.170   m , x = 0.025 , + 0.025   m , and y = 0.100 , + 0.100   m , at t = 0   s , as shown in Figure 2. The porous medium is only a saturated zone and is composed of spherical glass beads equivalent to sand [31]. All the boundary conditions are no-flow except for the infiltration zone on top of the parallelepiped and an impermeable zone at the bottom. The legend at the right-hand side of Figure 4 indicates the saturation contour values of DNAPL, σ n = S n ϕ , in color bars. A hydraulic gradient of 0.02083 originates a flow from the right to the left (Figure 4).
Table 1 gives the material properties and parameter details used in the numerical simulations performed with CactusHydro code [17,18], and the initial conditions are shown in Figure 4. Most parameters are similar to the column case, since we used this case to calibrate the sandbox model, except for the hydraulic conductivity, which changed slightly to be K = 6.92 × 10 3   m / s . From here, we obtain the absolute permeability, k = K μ w ρ w g , that is k = 7.061 × 10 10   m 2 , over the whole grid, with k x = k y = k z .

3.2. Comparison of DNAPL Migration between Experiment and Numerical Simulations in a Column

After being released for three seconds, using an injector placed in the saturated zone, the DNAPL migrates downward to the saturated zone under the influence of gravity (see Figure 5), which shows a comparison between the experimental results inside the column (left-hand side of each of the four panels) and the numerical results of the saturation σ n contours of the DNAPL migration in the z x plane (the right-hand side of each panel) for different times. The σ n = S n ϕ , is the product of the saturation of the DNAPL multiplied by the porosity of the porous medium, and the color bar indicates its value. The column is represented using a grid dimension of 0.060 m × 0.060 × 0.360   m at different times: 12 s, 24 s, 50 s, 90 s. The left-hand side of each panel time shows the experimental result, while the saturation contours of the DNAPL migration in the z x plane are placed on the right-hand side of each panel. Although we have numerical outputs for any particular time, we have chosen two equispaced time frame (12 s and 24 s), plus one when the contaminant reaches the medium part of the column (50 s), and the last one corresponding to the contaminant reaching the bottom part of the column (90 s).
The first panel shows the time at 12 s, the second at 24 s, the third at 50 s, and the fourth at 90 s, respectively. The contaminant is released for three seconds with a pressure of 4118.0   P a . Then, it moves downward (no hydraulic gradient is present), but lateral spreading also occurs due to the effect of the viscosity/capillary pressures (see Table 1 for details). Notice how the contaminant fills the entire column at each time.
The initial time in which the contaminant arrives at the bottom is at approximately 90 s (fourth panel). The final time it completely arrives at the bottom is after around 220 s. For each time, the transient saturation contours are viewed in the z x (see Figure 5). The agreement between the experimental and the numerical simulation results is very good. In Table 2, we show the values used in Figure 5.

3.3. Comparison of DNAPL Migration between Experiment and Numerical Simulations in a Sandbox

After being released for forty-five seconds, using an injector placed in the saturated zone (see Figure 4) that released the DNAPL hde7100 at a pressure of 3118.0   P a , it migrates downward into the saturated zone under the influence of gravity. In Figure 6, a comparison between the experimental results of the sandbox (the left-hand side of each of the four panels) and the numerical results of the saturation σ n contours of the DNAPL migration in the z x plane (the right-hand side of each panel) for four different times are shown. Although we have numerical outputs for any particular time, we have chosen two equispaced time frames (5 s and 25 s), plus one when the contaminant reaches the medium part of the column (50 s), and the last one corresponding to the contaminant reaching the bottom part of the column (75 s).
The first panel shows the time at 5 s, the second at 25 s, the third at 50 s, and the fourth at 75 s, respectively. The injection of the DNAPL into the saturated zone is performed in both x directions (see the first pane, the right-hand side). There is a hydraulic gradient and an underline motion with Q = 9.4   m L / s in the left direction. A green rectangle has been highlighted for comparison with the real experimental image.
After 75 s, the contaminant arrives at the bottom of the green rectangle (last panel, the right-hand side), similar to the experimental counterpart on the left-hand side. The comparison between both the numerical and experimental results is very good, even though the numerical model simulates a slightly faster migration than what was observed in the sandbox (Table 3). At the same time, the difference between the observed and the simulated transport velocity progressively decreases during the migration of the contaminant in the subsurface. In Table 3, we show the values used in Figure 6.

4. Conclusions

This paper investigates the migration of a DNAPL HFE-7100 through a column and a sandbox laboratory scale porous media and compares their results with three-dimensional numerical simulations using a CactusHydro [17,18] conservative HRSC method that precisely follows the advective part of the fluid flow and resolves the hyperbolic part of the nonlinear governing coupled PDFs of a three-phase immiscible fluid flow. The contaminant is released using a colored fluid as a tracer injected into a saturated zone, first using a column laboratory experiment, then using a sandbox-scale system with a hydraulic gradient and a fixed time.
We investigated the temporal evolution of the migration of the immiscible DNAPL in a porous medium following the saturation contour profiles of the three-phases fluids flow (air is zero in the saturated zone) numerical simulations. The comparison with the experimental results shows very good agreement, even though the numerical model simulates a slightly faster migration than that observed in the sandbox. At the same time, the difference between the observed and the simulated transport velocity progressively decreases during the migration of the contaminant in the subsurface.
The next step required having quantitative DNAPL saturation data from the experimental images and comparing them with the numerical data for these homogeneous conditions. This means quantifying the images in terms of DNAPL saturations. Future developments could include the application of experimental tests carried out in a sandbox filled with heterogeneous porous media; they can lead to a detailed analysis of the fluid flow phenomena and their relationship with the capillary pressures.
In a wider context, the experimental validation of this numerical approach further confirms the possibility of applying it to (i) reconstruct reliable environmental scenarios caused by DNAPL releases in the subsurface (with emphasis on porous media) and (ii) design (also through numerical simulations) the best hydraulic barrier system to be constructed to optimize the free-product DNAPL extraction in potential emergency scenarios.
The limitation of the work is that, at this stage, the experimental validation was made at the laboratory scale in homogeneous porous media. Therefore, further validations will be carried out by (i) quantifying the images in terms of DNAPL saturations in homogeneous and heterogeneous conditions at the laboratory scale and (ii) at a site scale by analyzing (and simulating) the effect of real DNAPLs releases in the subsurface.

Author Contributions

Conceptualization, A.F., F.C. and A.Z.; methodology, A.F., F.C. and A.Z.; experiment, A.Z.; numerical, A.F.; validation, A.F., F.C. and A.Z.; writing and reviewing the draft, A.F., F.C. and A.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data is contained within the article.

Acknowledgments

This work used resources at the University of Parma at (https://www.hpc.unipr.it). This research benefited from the equipment and framework of the COMP-R Initiative, funded by the ‘Departments of Excellence’ program of the Italian Ministry for University and Research (MUR, 2023–2027). A.F. and F.C. acknowledge financial support from the PNRR MUR project ECS_00000033_ECOSISTER. We also thank two anonymous reviewers for their interesting comments/questions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Goswami, R.R.; Clement, T.P. Laboratory-scale investigation of saltwater intrusion dynamics. Water Resour. 2007, 43, W04418. [Google Scholar] [CrossRef]
  2. Silliman, S.E.; Zheng, L. Comparison of observations from a laboratory model with stochastic theory: Initial analysis of hydraulic and tracer experiments. Transp. Porous Media 2001, 42, 85–107. [Google Scholar] [CrossRef]
  3. Rathfelder, K.M.; Abriola, L.M.; Taylor, T.P.; Pennell, K.D. Surfactant enhanced recovery of tetrachloroethylene from a porous medium containing low permeability lenses 2. Numerical simulation. J. Contam. Hydrol. 2001, 48, 351–374. [Google Scholar] [CrossRef]
  4. Luciano, A.; Viotti, P.; Papini, M.P. Laboratory investigation of DNAPL migration in porous media. J. Hazard. Mater. 2010, 176, 1006–1017. [Google Scholar] [CrossRef] [PubMed]
  5. Luciano, A.; Mancini, G.; Torretta, V.; Viotti, P. An empirical model for the evaluation of the dissolution rate from a DNAPL-contaminated area. Environ. Sci. Pollut. Res. 2018, 25, 33992–34004. [Google Scholar] [CrossRef]
  6. Citarella, D.; Cupola, F.; Tanda, M.G.; Zanini, A. Evaluation of dispersivity coefficients by means of a laboratory image analysis. J. Contam. Hydrol. 2015, 172, 10–23. [Google Scholar] [CrossRef] [PubMed]
  7. Seyedpour, S.M.; Valizadeh, I.; Kirmizakis, P.; Doherty, R.; Ricken, T. Optimization of the Groundwater Remediation Process Using a Coupled Genetic Algorithm-Finite Difference Method. Water 2021, 13, 383. [Google Scholar] [CrossRef]
  8. Aksoy, A.; Guney, M.S. Experimental determination of three-dimensional dispersivities in homogeneous porous medium. Environ. Earth Sci. 2010, 60, 383–393. [Google Scholar] [CrossRef]
  9. Castro-Alcala, E.; Fernàndez-Garcia, D.; Carrera, J.; Bolster, D. Visualization of mixing processes in a heterogeneous sand box aquifer. Environ. Sci. Technol. 2012, 46, 3228–3235. [Google Scholar] [CrossRef]
  10. Praseeja, A.V.; Sajikumar, N. A review on the study of immiscible fluid flow in unsaturated porous media: Modeling and remediation. J. Porous Media 2019, 22, 889–922. [Google Scholar] [CrossRef]
  11. Soga, K.; Page, J.W.E.; Illangasekare, T.H. A review of NAPL source zone remediation efficiency and the mass flux approach. J. Hazard. Mater. 2004, 110, 13–27. [Google Scholar] [CrossRef]
  12. Engelmann, C.; Lari, K.S.; Schmidt, L.; Werth, C.J.; Walther, M. Towards predicting DNAPL source zone formation to improve plume assessment: Using robust laboratory and numerical experiments to evaluate the relevance of retention curve characteristics. J. Hazard. Mater. 2021, 407, 124741. [Google Scholar] [CrossRef] [PubMed]
  13. Zheng, F.; Gao, Y.; Sun, Y.; Shi, X.; Xu, H.; Wu, J. Influence of flow velocity and spatial heterogeneity on DNAPL migration in porous media: Insights from laboratory experiments and numerical modelling. Hydrogeol. J. 2015, 23, 1703–1718. [Google Scholar] [CrossRef]
  14. Parker, B.L.; Cherry, J.A.; Chapman, S.W. Field study of TCE diffusion profiles below DNAPL to assess aquitard integrity. J. Contam. Hydrol. 2004, 74, 197–230. [Google Scholar] [CrossRef]
  15. Harbaugh, A.W.; Banta, E.W.; Hill, M.C.; McDonald, M.G. MODFLOW-2000, the U.S. Geological Survey Modular Ground-Water Model—User Guide to Modularization Concepts and the Ground-Water Flow Process; Open File Report 00-92; United States Geological Survey: Reston, VA, USA, 2000.
  16. Zheng, C.; Wang, P.P. MT3DMS: A Modular Three-Dimensional Multispecies Transport Model for Simulation of Advection, Dispersion, and Chemical Reactions of Contaminants in Groundwater Systems; Documentation and User’s Guide No. SERDP-99-1; U.S. Army Engineer Research and Development Center: Vicksburg, MS, USA, 1999.
  17. Feo, A.; Celico, F. High-resolution shock-capturing numerical simulations of three-phase immiscible fluids from the unsaturated to the saturated zone. Sci. Rep. 2021, 11, 5212. [Google Scholar] [CrossRef]
  18. Feo, A.; Celico, F. Investigating the migration of immiscible contaminant fluid flow in homogeneous and heterogeneous aquifers with high-precision numerical simulations. PLoS ONE 2022, 17, e0266486. [Google Scholar] [CrossRef] [PubMed]
  19. Kurganov, A.; Tadmor, E. New high-resolution central scheme for non-linear conservation laws and convection-diffusion equations. J. Comput. Phys. 2000, 160, 241–282. [Google Scholar] [CrossRef] [Green Version]
  20. Lax, P.; Wendroff, B. Systems of conservation laws. Commun. Pure Appl. Math. 1960, 3, 217–237. [Google Scholar] [CrossRef] [Green Version]
  21. Hou, T.Y.; LeFloch, P.G. Why nonconservative schemes converge to wrong solutions: Error analysis. Math. Comp. 1994, 62, 497–530. [Google Scholar] [CrossRef]
  22. Allen, G.; Goodale, T.; Lanfermann, G.; Radke, T.; Rideout, D.; Thornburg, J. Cactus Users’ Guide. 2011. Available online: http://www.cactuscode.org/documentation/UsersGuide.pdf (accessed on 1 January 2023).
  23. Cactus Developers. Cactus Computational Toolkit. Available online: http://www.cactuscode.org/ (accessed on 1 January 2023).
  24. Goodale, T.; Allen, G.; Lanfermann, G.; Massó, J.; Radke, T.; Seidel, E.; Shalf, J. The Cactus Framework and Toolkit: Design and Applications. In Vector and Parallel Processing—VECPAR’2002, Proceedings of the 5th International Converence, Porto, Portugal, 26–28 June 2002; Lecture Notes in Computer Science; Springer: Berlin/Heidelberg, Germany, 2003. [Google Scholar]
  25. Schnetter, E.; Hawley, S.H.; Hawke, I. Evolutions in 3D numerical relativity using fixed mesh refinement. Class. Quantum Gravity 2004, 21, 1465–1488. [Google Scholar] [CrossRef] [Green Version]
  26. Schnetter, E.; Diener, P.; Dorband, E.N.; Tiglio, M. A multi-block infrastructure for three-dimensional time-dependent numerical relativity. Class. Quantum Gravity 2006, 23, S553. [Google Scholar] [CrossRef] [Green Version]
  27. Faust, C.R. Transport of immiscible fluids within and below the unsaturated zone: A numerical model. Water Resour. Res. 1985, 21, 587–596. [Google Scholar] [CrossRef]
  28. Faust, C.R.; Guswa, J.H.; Mercer, J.W. Simulations of three-dimensional flow of immiscible fluids within and below the unsaturated zone. Water Resour. Res. 1989, 25, 2449–2464. [Google Scholar] [CrossRef]
  29. Parker, J.C.; Lenhard, R.J.; Kuppusamy, T. A parametric model for constitutive properties governing multi-phase flow in porous media. Water Resour. Res. 1987, 23, 618–624. [Google Scholar] [CrossRef]
  30. van Genuchten, M.T. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 1980, 44, 892–898. [Google Scholar] [CrossRef] [Green Version]
  31. Chiapponi, L. Water retention curves of multicomponent mixtures of spherical particles. Powder Technol. 2017, 320, 646–655. [Google Scholar] [CrossRef]
  32. Hu, X.; Meng, X.; Wei, K.; Li, W.; Wu, J. Compressed Liquid Viscosity Measurements of HFE-7000, HFE-7100, HFE-7200, and HFE-7500 at temperatures from (253 to 373) K and pressures up to 30 Mpa. J. Chem Eng. Data 2015, 60, 3562–3570. [Google Scholar] [CrossRef]
  33. Freeze, R.A.; Cherry, J.A. Groundwater Book; Prentice-Hall, Inc.: Englewood Cliffs, NJ, USA, 1979. [Google Scholar]
  34. Carsel, R.F.; Parrish, R.S. Developing Joint Probability Distributions of Soil Water Retention Characteristics. Water Resour. Res. 1988, 24, 755–769. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Experimental equipment for the column filled with glass beads of 1 mm diameter as the porous media and saturated with water and fluorescein.
Figure 1. Experimental equipment for the column filled with glass beads of 1 mm diameter as the porous media and saturated with water and fluorescein.
Water 15 01471 g001
Figure 2. Sketch of the experimental device. Constant head boundaries upstream and downstream; the porous media was laterally confined by iron plates, and at the top, it is involved in the capillary fringe. The red dot is the source location. Dimensions are in mm. The blue arrow indicates the flow direction.
Figure 2. Sketch of the experimental device. Constant head boundaries upstream and downstream; the porous media was laterally confined by iron plates, and at the top, it is involved in the capillary fringe. The red dot is the source location. Dimensions are in mm. The blue arrow indicates the flow direction.
Water 15 01471 g002
Figure 3. The grid geometry used in the numerical simulation of an immiscible DNAPL inside a saturated column and a grid dimension of 0.060 m × 0.060 × 0.360 m, at the initial time t = 0 s. The red box is the immiscible DNAPL at the top of the parallelepiped in the z x plane (left-hand side) and the z y plane (right-hand side), respectively.
Figure 3. The grid geometry used in the numerical simulation of an immiscible DNAPL inside a saturated column and a grid dimension of 0.060 m × 0.060 × 0.360 m, at the initial time t = 0 s. The red box is the immiscible DNAPL at the top of the parallelepiped in the z x plane (left-hand side) and the z y plane (right-hand side), respectively.
Water 15 01471 g003
Figure 4. The grid geometry used in the numerical simulation of an immiscible DNAPL inside a saturated sandbox and a grid dimension of 1.16 m × 0.10 × 0.80 m, at the initial time t = 0 s. The blue box is the immiscible DNAPL at the top of the parallelepiped in the z x plane.
Figure 4. The grid geometry used in the numerical simulation of an immiscible DNAPL inside a saturated sandbox and a grid dimension of 1.16 m × 0.10 × 0.80 m, at the initial time t = 0 s. The blue box is the immiscible DNAPL at the top of the parallelepiped in the z x plane.
Water 15 01471 g004
Figure 5. Comparison between the experimental results on the DNAPL migration in a column and three-dimensional numerical results on the saturation contours σ n = S n ϕ of the DNAPL in a column using a grid dimension of 0.060   m × 0.060 × 0.360   m , at different times: 12 s, 24 s, 50 s, 90 s. The left-hand side of each panel time shows the experimental result, while the saturation contours of the DNAPL migration in the z x plane are placed on the right-hand side of each panel.
Figure 5. Comparison between the experimental results on the DNAPL migration in a column and three-dimensional numerical results on the saturation contours σ n = S n ϕ of the DNAPL in a column using a grid dimension of 0.060   m × 0.060 × 0.360   m , at different times: 12 s, 24 s, 50 s, 90 s. The left-hand side of each panel time shows the experimental result, while the saturation contours of the DNAPL migration in the z x plane are placed on the right-hand side of each panel.
Water 15 01471 g005
Figure 6. Comparison between the experimental results on the DNAPL migration in a sandbox and three-dimensional numerical results on the saturation contours σ n = S n ϕ of the DNAPL in a sandbox using a grid dimension of 0.060   m × 0.060 × 0.360   m at different times: 5 s, 25 s, 50 s, 75 s, respectively. The left-hand side of each panel time shows the experimental result in the sandbox, while the saturation contours of the DNAPL migration in the z x plane are placed on the right-hand side of each panel.
Figure 6. Comparison between the experimental results on the DNAPL migration in a sandbox and three-dimensional numerical results on the saturation contours σ n = S n ϕ of the DNAPL in a sandbox using a grid dimension of 0.060   m × 0.060 × 0.360   m at different times: 5 s, 25 s, 50 s, 75 s, respectively. The left-hand side of each panel time shows the experimental result in the sandbox, while the saturation contours of the DNAPL migration in the z x plane are placed on the right-hand side of each panel.
Water 15 01471 g006
Table 1. List of parameters used for the three-dimensional numerical simulations of a DNAPL leak in a column.
Table 1. List of parameters used for the three-dimensional numerical simulations of a DNAPL leak in a column.
ParameterSymbolValue
Absolute permeability, m 2 k 6.122 × 10 10
Rock compressibility, P a 1 c R 4.35 × 10 7
Porosity ϕ 0 0.37
Water viscosity, k g / m s μ w 10 3
Water density, k g / m 3 ρ w 10 3
DNAPL HFE-7100 dynamic viscosity, k g / m s μ n 1.35 × 10 3
DNAPL HFE-7100 density, k g / m 3 ρ n 1500
Air viscosity, k g / m s μ a 1.8 × 10 5
Air density, k g / m 3 ρ a 1.225
Van Genuchten parameter n , m 2.68,0.627
Irreducible wetting phase saturation S w i r 0.045
Surface tension DNAPL, N / m σ n a 13.60 × 10 3
Interfacial tension DNAPL, N / m σ n w 35.59 × 10 3
Surface tension water, N / m σ a w 71.75 × 10 3
Capillary pressure air-water at zero saturation, P a p c a w 0 676.55
Capillary pressure DNAPL-water at zero saturation, P a p c n w 0 334.93
Capillary pressure air-nonaqueous at zero saturation, P a p c a n 0 341.62
Resolution, m Δ x = Δ y = Δ z 0.01
Table 2. Values of the position of the contaminant migration in the column experiment, and in the numerical model output, for different times.
Table 2. Values of the position of the contaminant migration in the column experiment, and in the numerical model output, for different times.
Time (s)Position of the Contaminant in the Experimental Column (m)Position of the Contaminant in the Numerical Model (m)
12 0.125 0.135
24 0.170 0.175
50 0.240 0.245
90 0.320 0.350
Table 3. Values of the position of the contaminant migration in the sandbox experiment, and in the numerical model output, for different times.
Table 3. Values of the position of the contaminant migration in the sandbox experiment, and in the numerical model output, for different times.
Time (s)Position of the Contaminant in the Experimental Sandbox (m)Position of the Contaminant in the Numerical Model (m)
5 0.200 0.225
25 0.280 0.335
50 0.380 0.435
75 0.435 0.480
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Feo, A.; Celico, F.; Zanini, A. Migration of DNAPL in Saturated Porous Media: Validation of High-Resolution Shock-Capturing Numerical Simulations through a Sandbox Experiment. Water 2023, 15, 1471. https://doi.org/10.3390/w15081471

AMA Style

Feo A, Celico F, Zanini A. Migration of DNAPL in Saturated Porous Media: Validation of High-Resolution Shock-Capturing Numerical Simulations through a Sandbox Experiment. Water. 2023; 15(8):1471. https://doi.org/10.3390/w15081471

Chicago/Turabian Style

Feo, Alessandra, Fulvio Celico, and Andrea Zanini. 2023. "Migration of DNAPL in Saturated Porous Media: Validation of High-Resolution Shock-Capturing Numerical Simulations through a Sandbox Experiment" Water 15, no. 8: 1471. https://doi.org/10.3390/w15081471

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop