Estimation of Free-Product PCE Distribution in Thick Multilayered Aquifers as Possible Long-Term Pollution Sources for Shallow and Deep Groundwaters, Using High-Precision Numerical Simulations

Abstract

Chlorinated organic compounds are Persistent Organic Pollutants (POPs) with high environmental continuity. These chemicals possess the ability to permeate into the environment across both unsaturated and saturated zones. This study examines the potential impact of perchloroethylene (PCE) releases in aquifer systems consisting of layers with varying permeability. The numerical simulations utilized the CactusHydro numerical code that employs a high-resolution shock-capturing flux conservative method to solve the non-linear partial differential equations of a three-phase immiscible fluid flow and study the migration of PCE into variably saturated zones to predict the spatial and temporal distribution of free products across multilayered aquifers. Two scenarios were simulated to predict how a first low-permeable horizon would affect the downward migration of PCE in its pure phase. The numerical simulations show that the multilayered aquifer system is vulnerable to an in-depth study of the PCE migration paths in the pure phase. However, very long times (about a year) are required for the DNAPL to be able to cross the semi-permeable horizons. The results in both scenarios show the quantity of mass (in kg and percentage concerning the initial one) that is left along the multilayered aquifer during the migration and after several years.

1. Introduction

Chlorinated organic compounds widely contaminate aquifer systems worldwide [1,2]. DNAPLs, which stand for dense, non-aqueous phase liquids, are substances that are denser than water and have very low water solubility. They can migrate under pressure and gravity forces through unsaturated and saturated media until reaching a bottom aquiclude (e.g., [3]). DNAPLs are typically found in areas impacted by human activity and are commonly used in dry cleaning, metal degreasing, and chemical production. According to the Stockholm Convention, several chlorinated organic compounds are classified as Persistent Organic Pollutants (POPs).

They persist in the environment and negatively impact human health (e.g., [4,5]). From 1960 onwards, numerous epidemiological studies have been conducted to assess the potential correlation between exposure to chlorinated compounds and the development of tumors, including non-Hodgkin lymphoma, bladder cancer, and multiple myeloma (e.g., [6]). In experimental studies, an increase in the incidence of mononuclear cell leukemia in rats was observed [7]. Additionally, an increase in kidney tumors was documented [8]. Male rats also exhibited higher rates of brain glioma and testicular interstitial cell tumors [6].

The behavior of chlorinated organic compounds in the subsurface has been under investigation since the 1980s (e.g., [9,10,11,12,13,14,15,16] for a review in modeling groundwater and surface water interactions). In detail, some studies were focused on relationships between DNAPLs and low-permeability media, such as clay aquitards (e.g., [17]), at sites where the groundwater quality will be significantly influenced by the integrity (sensu [18]) of the aquitards. As a matter of fact, low-permeability rocks are linked to high capillary resistance to DNAPLs, which causes the formation of the so-called DNAPL pools on top of very low-permeability aquitards [3]. Some studies were carried out at a laboratory scale to analyze possible factors that decrease the whole integrity of low-permeability media. For example [19,20], studying glaciolacustrine clay and clay till, respectively, suggested that fractures can serve as rapid conduits for DNAPL migration and significantly influence the whole DNAPL dynamics. DNAPL migration in fractured, low-permeability rocks is typically characterized by matrix diffusion, and the dual-porosity attribute can provide significant retention capacity. This is because the mass storage capacity of the matrix is generally much higher than that of the fractures [21].

Just as Sun studied the mechanical behavior of hard particles in rock fractures in [22], the movement of PCE in aquifer fractures is also affected by factors such as fracture angle, width, and roughness. Sun [22] points out that the trajectory shape, surface roughness, and opening angle of fractures have important influences on the mechanical action of particles. This also applies to our study of PCE migration in aquifers. Because the heterogeneity and complexity of aquifers will lead to the uncertainty of PCE migration, we need to consider these factors to more accurately predict the distribution of PCE.

To more accurately simulate the geological materials in the aquifer system, there is the joint particle model and rotation calculation model proposed by Sun and Huang [23] to solve the problems existing in traditional simulations, such as simply regarding soil particles as circular or elliptical. The research results on the influence of soil particle gradation and porosity on material properties in [23] have important implications for our understanding of the migration and distribution of PCE in aquifers because the properties of aquifers will directly affect the migration path and pollution range of PCE.

The following three-dimensional numerical study examines the potential effects of releasing perchloroethylene (PCE) in aquifer systems consisting of alternating layers of low-permeability (silt and clay) and high-permeability (gravel and sand) materials. The numerical code was introduced in [24,25] and validated through a sandbox laboratory experiment [26] and several applications that include free product extraction of DNAPLs in potential emergency scenarios [27,28,29,30]. It is based on the high-resolution shock-capturing (HRSC) flux conservative method [30,31,32], which accurately follows sharp discontinuities (shock fronts) and explicitly evolves the discretized (in space) variables. CactusHydro is written on the Cactus computational toolkit [33,34], an open-source problem-solving environment designed to enable parallel computing using simfactory, high-performance computing (HPC) simulation codes. The data evolved on a Cartesian mesh using Carpet [35,36]. The migration of a spilled DNAPL is considered immiscible, and the effects of volatilization, biodegradation, or dissolution are not considered. CactusHydro treats the vertical and horizontal movement of the contaminant in the variably saturated zone as a whole and it is numerically resolved as a unique zone (not separating the vertical movement from the horizontal one since the flow equation includes both zones).

The numerical code simulates the DNAPL migration through a wedge-shaped (from a few meters up to several hundred meters thick) multilayered aquifer starting from the top to a low-permeability layer and an alternation between low-permeability layers and more permeable ones. The simulation aims to predict the distribution of free products and the potential long-term pollution sources in both shallow and deep groundwaters within the aquifer. In detail, two main scenarios were simulated, with different characteristics of the different layers. The objective is to understand if and how some of the layer’s features (e.g., grain size, thickness, inclination, permeability) or the hydraulic gradient can influence the migration of DNAPLs in the subsurface. The test site is the alluvial aquifer of Parma Plain (Northern Italy), where perchloroethylene (PCE) was recently detected in groundwater [37,38], suggesting the existence of DNAPL sources. A previous study [39] presented a scenario with a smaller grid to investigate the influence of a low-permeability layer on the vertical migration of DNAPLs between two higher permeability horizons.

2. Hydrogeological Features and PCE Contamination at the Study Area

The study area is the alluvial plain located east of Parma town. From a geological point of view, the alluvial plain deposits are part of a well-known marine-regressive sedimentary succession [40,41]. This succession is composed of several synthems (sensu [42]). It starts with the Pliocene marine clay (FAA), which is discordantly eroded and onlapped by Pleistocene shallow marine, fan-deltas, and lagoon deposits (ATS and CMZ synthems) progressively evolving to continental depositional environments (AEI and AES synthems) since middle Pleistocene to Holocene (Figure 1).

In this work, the simulation of two scenarios is framed into a multilayer aquifer made of Pleistocene alluvial units (parts of AES synthem) that are characterized by alternating gravel-sand and silt silty-clay units (Figure 2). These alternations are controlled by the depositional dynamics of the Apennine paleo-streams and Po River, in combination with climatic–eustatic alternations related to the latest glacial and interglacial events. The aforementioned geological and, consequently, hydraulic characteristics of each individual layer depend on the locations of the single paleo-stream that deposited them. Therefore, the horizontal and vertical heterogeneity of the multilayer aquifer can be explained by the variability of the depositional environments and sub-environments. In addition, the inclination and the wedge shape of the sedimentary bodies and, in general, their geometries can also be influenced by tectonic buried thrust fronts, which create a northward thickening succession.

At a basin scale and in both simulated contexts, groundwater flows from southwest to northeast with depths that may be included between zero and several tens of meters, depending on the character of the aquifer layer and its location in the alluvial plain [44]. It is recharged by (i) local effective infiltration, which is maximum in the Apennine foothills area, and (ii) losing streams and artificial channels, which are hydraulically connected to the groundwater (e.g., [37,38,44,45,46]). During the AMIIGA CE32 research project (UE INTERREG program), some well tests have been performed at new monitoring wells to evaluate the hydraulic properties of the aquifer. Zanini et al. [37,38] estimated values of the hydraulic conductivity ranged from 1.2 × 10⁻⁵ to 4.9 × 10⁻⁵ m/s (mean 2.3 × 10⁻⁵ m/s; median 1.7 × 10⁻⁵ m/s) in gravel and sand horizons, and from 9.3 × 10⁻⁹ to 5.3 × 10⁻⁷ m/s (mean 1.6 × 10⁻⁷ m/s; median 9.7 × 10⁻⁸ m/s) in silt and clay layers.

Considering the available information, the calibration of the hydraulic conductivity process started from the knowledge of a few hydraulic tests available in the study area [38] and a well test performed on monitoring wells. The parameter estimation process was carried out through the software PEST (version 3.5) [37,38]. The hydraulic conductivity of the outcropping low-permeability horizon was set to 5 × 10⁻⁷ m/s [37]. The estimated hydraulic conductivity of the more permeable layer is (1.8–7.1 × 10⁻⁵ m/s) [38] and is consistent with the alluvial nature of the aquifer.

From 2002 to 2019, PCE was detected in groundwater with concentrations up to 23 µg/L [38,39] using monitoring wells (from Parma Municipality) and piezometers drilled within the AMIIGA CE32 research project (UE INTERREG program). The PCE concentrations cyclically varied over time, with the highest detected in early recharge and the smaller concentrations in late recession [27]. For more than two decades, the presence of PCE has been detected. The cyclic variation in its concentrations over time and the consistent peak concentration observed in different hydrologic years indicate the presence of one or more free-product DNAPL pools upgradient of the observation wells mentioned in [38]. Therefore, taking into consideration (i) the layered hydrogeological setting of the whole system, as well as (ii) the existence of low-permeability layers interbedded within the coarse-grained ones, a 3D high-precision model has been implemented for PCE migration, forecasting its possible propagation towards the deep aquifer.

2.1. Mathematical Setup and Governing Equations

Here, it briefly reviews the method used to fix the notation that will be used in the following sections. The equation that describes a three-dimensional non-aqueous ($n$), water ($w$), air ($a$), and phase fluid flow and the unsaturated/saturated zone introduced in Refs. [24,25] is given as follows:

$$\frac{\partial }{\partial t}\left({\mathsf{\rho }}_n\mathsf{\varphi }{S}_n\right)=\frac{\partial }{\partial x^i}\left[{\mathsf{\rho }}_n\frac{k_{rn}}{{\mathsf{\mu }}_n}{k}^{ij}\left(\frac{\partial {\mathrm{p}}_{\mathrm{a}}}{\partial x^j}+{\mathsf{\rho }}_{n}g\frac{\partial z}{\partial x^j}\right)\right]-\frac{\partial }{\partial x^i}\left[{\mathsf{\rho }}_n\frac{k_{rn}}{{\mathsf{\mu }}_n}{k}^{ij}\left(\frac{\partial {\mathrm{p}}_{\mathrm{can}}}{\partial x^j}\right)\right]+q_n,$$

(1)

$$\frac{\partial }{\partial t}\left({\mathsf{\rho }}_w\mathsf{\varphi }{S}_w\right)=\frac{\partial }{\partial x^i}\left[{\mathsf{\rho }}_w\frac{k_{rw}}{{\mathsf{\mu }}_w}{k}^{ij}\left(\frac{\partial {\mathrm{p}}_{\mathrm{a}}}{\partial x^j}+{\mathsf{\rho }}_{w}g\frac{\partial z}{\partial x^j}\right)\right]-\frac{\partial }{\partial x^i}\left[{\mathsf{\rho }}_w\frac{k_{rw}}{{\mathsf{\mu }}_w}{k}^{ij}\left(\frac{\partial {\mathrm{p}}_{\mathrm{caw}}}{\partial x^j}\right)\right]+q_w,$$

(2)

$$\frac{\partial }{\partial t}\left({\mathsf{\rho }}_a\mathsf{\varphi }{S}_a\right)=\frac{\partial }{\partial x^i}\left[{\mathsf{\rho }}_a\frac{k_{ra}}{{\mathsf{\mu }}_a}{k}^{ij}\left( \frac{\partial {\mathrm{p}}_{\mathrm{a}}}{\partial x^j}+{\mathsf{\rho }}_{a}g\frac{\partial z}{\partial x^j}\right)\right]+q_a,$$

(3)

where $x^i=\left(x,y,z,t\right)$ are the spatial and time cartesian coordinates. A fourth equation for the saturations is given by

$$S_w+S_n+S_a=1.$$

(4)

The systems (1)–(4) are written in terms of the saturations, $S_{\mathsf{\alpha }},\mathsf{\alpha }=\left(w,n,a\right)$; the air pressure, $p_a$; the capillary pressures of air–non-aqueous phase, $p_{can}=\left(p_a–p_n\right)$ and air–water phase, $p_{caw}=\left(p_a-p_w\right)$; the absolute permeability tensor, $k^{ij}$; the gravity acceleration, $g$; depth, mass source/sink, and the porosity, $z,q_{\mathsf{\alpha }},\mathsf{\varphi },$ respectively. Moreover, it depends on the densities and dynamics viscosity of each phase, ${\mathsf{\rho }}_{\mathsf{\alpha }},{\mathsf{\mu }}_{\mathsf{\alpha }}$, respectively.

Let us define the product of porosity and saturation as ${\mathsf{\sigma }}_w\equiv \mathsf{\varphi }{S}_w,{\mathsf{\sigma }}_n\equiv \mathsf{\varphi }{S}_n,{\mathsf{\sigma }}_a\equiv \mathsf{\varphi }{S}_a$. The porosity depends on the pressure as $\mathsf{\varphi }={\mathsf{\varphi }}_0\left[1+c_R\left(p-p_0\right)\right],$ where $c_R$ is the rock compressibility, and ${\mathsf{\varphi }}_0$ is the porosity at the atmospheric value $p_0$.

Equation (4) can be written as ${\mathsf{\sigma }}_a+{\mathsf{\sigma }}_n+{\mathsf{\sigma }}_w={\mathsf{\varphi }}_0\left[1+c_R\left(p-p_0\right)\right]$. Its derivative becomes

$$\frac{\partial {\mathsf{\sigma }}_a}{\partial t}+\frac{\partial {\mathsf{\sigma }}_n}{\partial t}+\frac{\partial {\mathsf{\sigma }}_w}{\partial t}={\mathsf{\varphi }}_0{c}_R\frac{\partial p}{\partial t},$$

(5)

for constant densities and viscosities. Using (1–3), we obtain

$$\frac{\partial {\mathsf{\sigma }}_{\left(\mathsf{\alpha }\right)}}{\partial t}=-\frac{\partial }{\partial x^i}\left[F_{\left(\mathsf{\alpha }\right)}^i\left(S_w,S_n,S_a,p\right)\right]+\frac{\partial }{\partial x^i}\left[Q_{\left(\mathsf{\alpha }\right)}^i\left(S_w,S_n,S_a,p\right)\right].$$

(6)

The coupled partial differential equation systems (5) and (6) as a function of $p,{\mathsf{\sigma }}_w,{\mathsf{\sigma }}_n,{\mathsf{\sigma }}_a$ must be resolved together with an explicit form of the relative permeabilities and capillary pressures. See [20,21] for details. In this paper, we introduced a new module for the initial conditions that reads multilayered surfaces (that will be explained in the following sections).

2.2. Hydrogeological Parameters of the Numerical Model

The geological settings of both the analyzed contexts implemented in the work simulations represent the physical basis for the simulation and have recreated interpolating punctual subsurface stratigraphic data mainly derived from water well boreholes (ISPRA and Regione Emilia-Romagna geo databases) using the Leapfrog Geo software (version 2021.1.3). This first step allowed us to obtain the three-dimensional top and bottom surfaces of the aquifers and aquitards object of the modeling, which will then be used to reconstruct the 3D grid for the numerical model using CactusHydro code [24,25]. In this paper, we consider one numerical model situated upstream of the alluvial aquifer of Parma with a hydraulic gradient of $1.47\times {10}^{-3}$ and another model situated downstream (see Figure 2) with a hydraulic gradient of $7.2\times {10}^{-4}$.

Table 1. Parameter values used in the numerical simulations of the migration of PCE in a variably saturated zone.

Parameter	Symbol	Value
Absolute permeability (gravel and sand)	$k_{gs}$	5.102 × 10−12 m2
Absolute permeability (silt and clay)	$k_{sc}$	5.102 × 10−14 m2
Rock compressibility	$c_R$	4.35 × 10−7 Pa−1
Porosity	${\mathsf{\varphi }}_0$	0.37
Water viscosity	${\mathsf{\mu }}_w$	10−3 kg/(ms)
Water density	${\mathsf{\rho }}_w$	103 kg/m3
DNAPL viscosity	${\mathsf{\mu }}_n$	0.844 × 10−3 kg/(ms)
DNAPL density	${\mathsf{\rho }}_n$	1643 kg/m3
Air viscosity	${\mathsf{\mu }}_a$	1.8 × 10−5 kg/(ms)
Air density	${\mathsf{\rho }}_a$	1.225 kg/m3
Van Genuchten	$\left(n,m\right)$	$\left(2.68,1-\frac{1}{2.68}\right)$
Irreducible wetting phase saturation	$S_{wir}$	0.045
Surface tension air–water	${\mathsf{\sigma }}_{aw}$	7.199 × 10−2 N/m
Interfacial tension non-aqueous–water	${\mathsf{\sigma }}_{nw}$	4.44 × 10−2 N/m
Capillary pressure air–water at zero saturation	$p_{caw0}$	676.55 Pa
Capillary pressure air–non-aqueous at zero saturation	$p_{can0}$	259.83 Pa

shows the parameter values used in the simulations of the different phases of fluid flow, including the PCE properties. For example, the PCE density is ${\mathsf{\rho }}_{\mathrm{n}}=1643 \mathrm{kg}/{\mathrm{m}}^3$, which is denser than water. We also considered two different values of the absolute permeability based on the values of the hydraulic conductivity previously mentioned. These values are $k_{gs}$ for a hydraulic conductivity equal to 5.0 × 10⁻⁵ m/s in a coarse-grained horizon (gravel and sand) and $k_{sc}$ for a hydraulic conductivity equal to 5.0 × 10⁻⁷ m/s in fine-grained layers (silt and clay). The relation between the absolute permeability and the hydraulic conductivity is given by $k=\frac{K{\mathsf{\mu }}_w}{{\mathsf{\rho }}_{w}g}$. Finally, the bottom of the aquifer is considered impermeable with a value of $k=5.102\times {10}^{-19}{\mathrm{m}}^2$ [47].

In our formulation, we use $p_{caw}=\left(p_a-p_w\right)$ and $p_{can}=\left(p_a-p_n\right)$ as the two independent capillary pressures for a three-phase fluid flow. The third one is $p_{cnw}=\left(p_n-p_w\right)=\left(p_{caw}-p_{can}\right).$ They are an extension from the two-phase expressions [11] (see [24,25] for details).

Both, the capillary pressures and relative permeabilities depend on the saturations, $k_{r\mathsf{\alpha }}=k_{r\mathsf{\alpha }}\left(S_a,S_n,S_w\right),p_{can}=p_{can}\left(S_a,S_n,S_w\right)$, and $p_{caw}=p_{caw}\left(S_a,S_n,S_w\right)$. We use the van Genuchten model [48] for the capillary pressures (see [24,25]). The van Genuchten parameter [11,48] for “sand” is as follows:

$$\mathsf{\alpha }={\left(\frac{p_c}{{\mathsf{\rho }}_{w}g}\right)}^{-1}=0.145 \mathrm{c}{\mathrm{m}}^{-1}=14.5 {\mathrm{m}}^{-1}$$

The capillary pressure for air–water is $p_{caw}=\frac{{\mathsf{\rho }}_{w}g}{\mathsf{\alpha }}=\frac{{10}^3 \mathrm{kg}/{\mathrm{m}}^3 9.81 \mathrm{m}/{\mathrm{s}}^2}{14.5 {\mathrm{m}}^{-1}}=676.55$ Pa. Using the interfacial tension ${\mathsf{\sigma }}_{nw}=4.44\times {10}^{-2} \mathrm{N}/\mathrm{m}$ and ${\mathsf{\sigma }}_{aw}=7.199\times {10}^{-2}\mathrm{N}/\mathrm{m}$, we obtain ${\mathsf{\beta }}_{nw}=\frac{{\mathsf{\sigma }}_{aw}}{{\mathsf{\sigma }}_{nw}}=1.62$ and $p_{cnw}\left(S_w\right)=\frac{p_{caw}}{{\mathsf{\beta }}_{nw}}=417.62$ Pa. The capillary pressure for air–non-aqueous is $p_{can}=p_{caw}-p_{cnw}=259.83$ Pa.

3. Three-Dimensional Numerical Results

3.1. Geological and Geometric Layers Found Characters

Although the two simulated scenarios are located in the same geological context and within a short distance of each other (about 3.0 km), several differences have been found to influence DNAPL migration in a singular way.

Table 2. Main geological characters founded in the “upstream scenario”.

Upstream Area
Layer	Prevalent Lithologies	Thickness (Average Value in m)	Inclination (°)
A1	Silt and clay	8.8	0.34
G1	Gravel and sand	49.1	1.29
A2	Silt and clay	9.3	1.58
G2	Gravel and sand	14.3	1.5
A3	Silt and clay	21.0	1.65
G3	Gravel and sand	36.0	1.72

and

Table 3. Main geological characters founded in the “downstream scenario”.

Downstream Area
Layer	Prevalent Lithologies	Thickness (Average Value in m)	Inclination (°)
A1	Silt and clay	16.8	0.19
G1	Gravel and sand	31.9	0.09
A2	Silt and clay	4.3	0.2
G2	Gravel and sand	13.7	0.09
A3	Silt and clay	11.5	0.17
G3	Gravel and sand	25.1	0.57

allow us to capture the main character of the layer and compare the two areas.

Layer thickness variation is linked to sedimentary dynamics of the ancient alluvial systems that originated these bodies. Climatic–eustatic cycles, rambling riverbeds, and the size dimensions of alluvial systems are among the main factors to refer to. The highest inclination in the “upstream scenario” is related to closer proximity to the Apennine growing chain and its associated buried thrust fronts. The angle generally appears to be increasing along with depth; older (and deeper) sediments have been subjected to tectonic uplift for a longer time compared to the younger ones that are more parallel to the current ground level.

3.2. Numerical Results of the Upstream Scenario of the Parma Alluvial Aquifer

The 3D numerical model describes two scenarios of a free-product PCE released into the environment that migrates downward through the multilayered aquifer and forecasts its possible propagation toward the deep aquifer. It is built on reconstructing by interpolating points of several surfaces that indicate the presence of low-permeability layers interbedded within the coarse-grained ones. We investigated two types of scenarios. The first one is located upstream in the south Parma alluvial aquifer, and we will call it an upstream scenario. At a basin scale, the groundwater flows from southwest to northeast [45] with a hydraulic gradient equal to $1.47\times {10}^{-3}.$ The second scenario is situated downstream in the north of the Parma alluvial aquifer, and we will define it as a downstream scenario. The groundwater has a hydraulic gradient of $7.2\times {10}^{-4}.$

Figure 3 shows the upstream scenario of the Parma alluvial aquifer, which is composed of eight surfaces using geo-referenced coordinates (and a grid resolution of 0.5 m). The absolute permeability between the red surface ${\mathrm{A}}_1$ and the brown surface ${\mathrm{G}}_1$, the cyan surface ${\mathrm{A}}_2$ and the purple surface ${\mathrm{G}}_2$, and the yellow surface ${\mathrm{A}}_3$ and the gray surface ${\mathrm{G}}_3$ is ${\mathrm{k}}_{\mathrm{gs}}$. The rest has an absolute permeability equal to ${\mathrm{k}}_{\mathrm{sc}}$ (see Table 1 for details). The green surface corresponds to the surface level. The surface colored in blue is the groundwater. The base aquiclude is impermeable, with a permeability value of $5.102\times {10}^{-19} {\mathrm{m}}^2$.

Before importing these surface data in the CactusHydro code, a grid for the CactusHydro simulation was defined by interpolating these data and rotating the surfaces so that the hydraulic gradient moves along the x-coordinate. Figure 4 shows the different surfaces (shown in Figure 3) for the upstream scenario, where the center of the surfaces corresponds to the value $\left(x,y\right)=\left(0,0\right)$. The z-coordinate elevation appears in the legend.

In the numerical model, a fixed volume of PCE is released in the subsurface, which then moves downward through the unsaturated zone located upstream of the alluvial Parma aquifer. See Figure 5 for the cubic grid, which also contains the interpolating eight surfaces (Figure 4). The system is a parallelepiped of 360.0 m long, from $x=\left[-180.0,180.0\right] \mathrm{m}$ on the left-hand side, 360.0 m wide from $y=\left[-180.0,180.0\right] \mathrm{m}$ on the right-hand side of Figure 5, and 146.0 m in the z-coordinate from $z=\left[49.0,-97.0\right] \mathrm{m}$. The spatial resolution is $dx=dy=dz=2.0 \mathrm{m}$, with a time step resolution of $dt=0.2 \mathrm{s}.$ The grid contains $\left(360\times 360\times 146\right)/8=2365200$ cells. The PCE released has a volume of $8 {\mathrm{m}}^3$ and saturation one. It is positioned on top of the grid at $\left(x,y,z\right)=\left(0,0,43.0\right) \mathrm{m}$. The boundary conditions are impermeable in the bottom (the green zone below). The imposed boundary conditions are such that there is no flow in the direction orthogonal to the side of the wall of the numeric grid. Indeed, to model the physical setting, we rotated the simulation domain so that the flow goes parallel to the x-coordinate at the edges (see Figure 4).

The system consists of an unsaturated zone (air and PCE) in light green in Figure 5, a groundwater table surface (in blue in Figure 3), and a saturated zone (depicted in pastel). The green color on the upper side is the atmosphere, while the bottom of the aquifer (Pliocene clay) represents the base aquiclude. The unsaturated zone is dry. The two pink rectangles on the figure’s upper side magnify the region at the bottom. The legend on the right-hand side indicates the saturation contours times the porosity of the PCE $\left({\mathsf{\sigma }}_n=S_n\times \mathsf{\varphi }\right)$ at different times (starting from time zero). We show the saturation in the $z–x$ plane (the left-hand side) and the $z–y$ plane (the right-hand side).

Once the PCE is released at the atmospheric pressure, it migrates downward into the dry, unsaturated zone under the influence of gravity. This is the first layer with a low permeability composed of silt and clay (see Table 1 for details). When the PCE reaches the groundwater table, it continues its migration downward, being denser than the water. See Figure 6, which shows the saturation contours ${\mathsf{\sigma }}_n=S_n\mathsf{\varphi }$ in three dimensions. It takes 47 days and 9.8 h to reach the coarse-grained surface.

The PCE now passes through the more coarse-grained layer and, after 312 days and 21.3 h, arrives at the second low-permeability layer situated approximately at z = −14 m. See Figure 7 for details. Notice that part of the PCE remains through the porous medium, having an irreducible water saturation equal to 0.045 (see Table 1). Most contaminants remain trapped in the unsaturated zone, although part continues its migration downward. It can also be noticed that the PCE slightly tends to move to the right-hand side (see left-hand side of Figure 7) due to a hydraulic gradient, although it prefers to migrate in-depth, as this layer is not completely impermeable.

After 1327 days and 9.8 h, the PCE enters the second coarse-grained permeability layer. See Figure 8 for details. It slightly accumulates just above the second low-permeability layer. Generally, moving inside this low-permeability layer takes longer than in the coarse-grained one. The PCE takes approximately one thousand days to cross this second low-permeability layer.

Figure 9 shows the fate of the PCE migration in the upstream scenario after 2094 days and 17.8 h (almost six years). The PCE has already crossed the second coarse-grained layer and arrived at the third low-permeability surface (in red). Notice that the contaminant generally accumulates above the second low-permeability layer and tends to move slightly in the right direction (see the left-hand side) following the hydraulic gradient, even if it prefers a vertical migration toward the depth.

Calculation of the Distribution of PCE Mass in the Upstream Aquifer Scenario

Now, let us investigate the contaminant distribution throughout the porous medium. See

Table 4. Mass (in kg) of the contaminant after 2094 days and 17.8 h in the upstream aquifer scenario and the contaminant’s percentual in the different z levels. Notice that the sum of the last column is 100%.

Upstream Aquifer Scenario
z Altitude (m. a.s.l)	Type of Layer	Mass in kg of PCE After 2094 Days	Percentage of Contaminant Trapped (%)
46 < z < 52	-	0.0	0.00
36 < z < 46	Silt and clay	1749.6	35.88
−14 < z < 36	Gravel and sand	1519.6	31.25
−22 < z < −14	Silt and clay	1282.2	26.34
−50 < z < −22	Gravel and sand	302.4	6.33
−58 < z < −50	Silt and clay	9.53	0.20
−90 < z < −58	Gravel and sand	0.0	0.00
−102 < z < −90	-	0.0	0.00

for details on the upstream scenario. At zero time, the PCE is located on top of the grid in (x,y,z) = (0,0,38) (see Figure 5) with a ${\mathsf{\sigma }}_n=S_n\mathsf{\varphi }=0.37$ and $S_n=1$, which means $M=4863.3 \mathrm{kg}$.

After 2094 days, the distribution of the PCE mass and the percentual concerning the initial distribution of the mass at time zero is shown in Table 4. Approximately 36% of contaminant remains trapped in the unsaturated zone, as seen in Figure 9. A similar percentage is located in the first coarse-grained and the second low-permeability layers, with 31.25% and 26.36%, respectively. The rest of the mass remains deeper within the immiscible fluid, and its mass remains constant. The migration continues downward, but because most contaminant is stuck in the porous medium, it eventually stops moving down further, especially in the deeper layers.

3.3. Numerical Results of the Downstream Scenario of Parma Alluvial Aquifer

A similar approach was used for the downstream scenario of the Parma alluvial aquifer. Figure 10 details the eight surfaces, including the piezometric level, with geo-referenced coordinates (and a grid resolution of 0.5 m). The colors are like the previous case, although the surfaces and locations are different. Indeed, the green surface corresponds to the surface level, while the groundwater table surface is colored in blue. Between the red surface $A_1$ and the brown surface $G_1$, $A_2$ (cyan) and $G_2$ (purple), and $A_3$ (yellow) and $G_3$ (gray), the absolute permeability is $k_{gs}$. At the same time, the rest has an absolute permeability of $k_{sc}$ and, thus, a low permeability (see Table 1 for details). The base aquiclude (as in the previous upstream scenario) is impermeable, with a permeability value of $5.102\times {10}^{-19} {\mathrm{m}}^2$.

Also here, like in the previous case, before importing these surface data in the CactusHydro code, we have defined the grid for the CactusHydro simulation by interpolating these data and rotating the surfaces so that the hydraulic gradient is moving along the x-coordinate. Figure 11 shows the different surfaces for the downstream scenario where the center of the surfaces corresponds to the value $\left(x,y\right)=\left(0,0\right)$. The z-coordinate elevation is represented on the legend of each image.

Now consider the numerical model for the downstream scenario of the Parma alluvial aquifer. A fixed volume of PCE (density ${\mathsf{\rho }}_{\mathrm{n}}=1643 \mathrm{kg}/{\mathrm{m}}^3$) is released in the subsurface as in the previous case. The interpolating eight surfaces in a cubic grid geometry are shown in Figure 12. The variably saturated zone is 360.0 m long, with $x=\left[-180.0,180.0\right] \mathrm{m}$ (left-hand side), 360.0 m wide with $y=\left[-180.0,180.0\right] \mathrm{m}$ (right-hand side), and 110.0 m in the z-coordinate from $z=\left[44.0,-66.0\right] \mathrm{m}$. Spatial resolution of $dx=dy=dz=2.0 \mathrm{m}$ and time step resolution of $dt=0.2 \mathrm{s}$ are used. The grid contains $\left(360\times 360\times 110\right)/8=1782000$ cells. The PCE with a volume of $8 {\mathrm{m}}^3$ and a saturation one is positioned on top of the grid at $\left(x,y,z\right)=\left(0,0,38.0\right) \mathrm{m}$.

When the PCE migrates downward into the unsaturated zone, it encounters the first low-permeability layer composed of silt and clay (see Table 1 for details). Since this first low-permeability layer is thicker than the previous case, the contaminant takes longer to arrive at the more coarse-grained layer. Indeed, 322 days and 8.9 h (see Figure 13). Due to a $S_{\mathrm{wir}}$ different from zero, most contaminants remain in the first low-permeability layer.

After 549 days and 22.2 h, the PCE crosses the more coarse-grained layer and arrives at the second low-permeability layer. See Figure 14 for details. The immiscible PCE moves slightly in the right direction due to a hydraulic gradient of $7.2\times {10}^{-4}$, as can be seen on the Z–X plane (while it is symmetric in the Z–Y plane), but continues to migrate deeper.

After 1042 days and 23.1 h, the PCE crossed the second low-permeability layer and arrived at the bottom of the second, more coarse-grained layer. See Figure 15 for details. It can also be noticed that the PCE accumulates on top of the second low-permeability layer, which is relatively thin.

Figure 16 shows the fate of the PCE migration after 2453 days and 23.5 h (almost seven years). The PCE crossed the third low-permeability layer. See Figure 16 for details. It also accumulates on top of the second and third low-permeability layers while migrating downward.

Calculation of the Residual PCE Mass in the Downstream Aquifer Scenario

The distribution of the PCE all over the porous medium in the downstream scenario is investigated. See

Table 5. Mass in Kg of the contaminant after 2585 days in the downstream aquifer scenario and the percentual of contaminant in the different levels in terms of z-coordinate. Notice that the sum is 100%.

Downstream Aquifer Scenario
z Altitude (m a.s.l.)	Type of Layer	Mass in kg of PCE After 2585 Days	Percentage of Contaminant Trapped (%)
39 < z < 47	-	0.0	0.00
23 < z < 39	Silt and clay	2093.8	43.06
−9 < z < 23	Gravel and sand	739.3	15.20
−13 < z < −9	Silt and clay	534.2	10.98
−27 < z < −13	Gravel and sand	771.4	15.86
−31 < z < −27	Silt and clay	625.4	12.86
−49 < z < −31	Gravel and sand	99.2	2.04
−71 < z < −49	-	0.0	0.00

for details on the downstream scenario. At time zero (see Figure 12), the contaminant is released into the unsaturated zone, and its mass is equal to 4863.3 Kg distributed on a volume of $8{ \mathrm{m}}^3$.

After 2585 days, approximately 43.05% of the contaminant remains trapped in the first low-permeability layer (partially unsaturated zone). It is bigger, around 10%, than the previous case of the upstream scenario since the first layer here is thicker. Also, for this reason, the rest of the mass reaches the more profound layer but in less quantity concerning the upstream scenario.

4. Conclusions

In this study, 3D numerical modeling is presented to explore the migration of a free-product PCE (DNAPL) released in a multilayered aquifer system. This aquifer system is characterized by low-permeability layers (silt and clay) alongside coarse-grained permeable layers (gravel and sand). The study’s primary objective was to forecast the distribution of the free product in multilayered aquifers and to comprehend the potential long-term sources of pollution for both shallow and deep groundwaters.

The attention was focused on the test site located in the alluvial aquifer of Parma Plain (northern Italy), where PCE was recently detected in groundwater [37,38], suggesting the existence of PCE sources. Two different scenarios were investigated, one in the upstream and the other in the downstream zone. Although they are situated very near each other, they present differences in the characterization of the different surfaces, although both start from a first low-permeability layer. The values used for the hydraulic conductivities were taken from calibrations (see Refs. [37,38]).

The numerical code CactusHydro (introduced in [24,25]) was used to develop 3D numerical modeling. It is based on the conservative high-resolution shock-capturing (HRSC) conservative method applied to groundwater that follows sharp discontinuities accurately and temporal dynamics of the fluid flow. The migration of the spilled DNAPL is considered immiscible, and the effects of volatilization, biodegradation, or dissolution are not considered. CactusHydro treats the vertical and horizontal movement of the contaminant in the variably saturated zone as a whole and is numerically resolved as a unique zone (not separating the vertical movement from the horizontal one since the flow equations include both zones).

The results indicate that CactusHydro can simulate the PCE contamination migration in space and time (see Table 4 and Table 5) and predict, among others, the free-product distribution in multilayered aquifers (up to several hundred meters thick) and the percentage distribution of possible pollution sources in the long-term scenario. The case study demonstrated that a potential release of PCE can migrate downward through low-permeability layers, even though a very low velocity was estimated and a high percentage of the mass is trapped within the migration pathway. These results suggest that the maximum depth at which the PCE can migrate through a thick multilayer aquifer can be greatly influenced by the aquifer features and the progressive contaminant trapping. At the same time, the DNAPL significantly trapped into the different layers (both lower- and higher-permeability ones) of a thick multilayer aquifer represents a widespread contamination source for groundwater within shallower and deeper parts of a single aquifer system.

The following steps of this research will examine, in multilayered systems, the role of the slope of low-permeability horizons in the contamination migration using the CactusHydro numerical code.