Next Article in Journal
Hybrid Optimization Algorithms of Firefly with GA and PSO for the Optimal Design of Water Distribution Networks
Next Article in Special Issue
Modeling of Distributed Control System for Network of Mineral Water Wells
Previous Article in Journal
Sludge Management in the Textile Industries of Bangladesh: An Industrial Survey of the Impact of the 2015 Standards and Guidelines
Previous Article in Special Issue
Migration of DNAPL in Saturated Porous Media: Validation of High-Resolution Shock-Capturing Numerical Simulations through a Sandbox Experiment
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Department of Chemistry, Life Sciences and Environmental Sustainability, University of Parma, 43124 Parma, Italy
*
Author to whom correspondence should be addressed.
Water 2023, 15(10), 1900; https://doi.org/10.3390/w15101900
Submission received: 17 March 2023 / Revised: 10 May 2023 / Accepted: 12 May 2023 / Published: 17 May 2023

Abstract

:
The contamination impact and the migration of the contaminant into the surrounding environment due to the presence of a spilled oil pipeline will cause significant damage to the natural ecosystem. For this reason, developing a rapid response strategy that might include accurate predictions of oil migration trajectories from numerical simulation modeling is decisive. This paper uses a three-dimensional model based on a high-resolution shock-capturing conservative method to resolve the nonlinear governing partial differential equations of the migration of a spilled light nonaqueous liquid oil contaminant in a variably saturated zone employed to investigate the migration of the oil pipeline leakage with great accuracy. The effects of the oil type density, gasoline, and diesel oil, the unsaturated zone depth, its saturation, the hydraulic gradient, and the pressure oil pipeline are investigated through the temporal evolution of the contaminant migration following the saturation profiles of the three-phase fluid flow in the variably saturated zone. The calculation results indicate that the leaking oil’s pressure is the parameter that significantly affects the contaminants’ arrival time at the groundwater table. Additionally, the water saturation of the unsaturated zone influences the arrival time, as the water saturation increases at a fixed depth. The unsaturated zone depth significantly influences the contaminant migration in the unsaturated zone. At the same time, the oil density and the hydraulic gradient have limited effects on the contaminant migration in the variably saturated zone.

1. Introduction

Oil and its derivates are usually transferred using pipelines that link oil fields to refineries or cargo oil tankers to oil storage, for example. Pipelines are considered the safest and most economical solution to transfer oil products over long distances. Still, problems of leakage have been frequently observed (e.g., [1,2]), with a negative impact on the environment (e.g., [3,4]). Construction defects, third-party damage, damages caused by landslides or earthquakes (e.g., [5,6]), as well as corrosion by the transmission fluid (e.g., [7,8]), can cause pipeline breakage.
Based on statistical analyses, the failure probability curve for oil and gas pipelines is divided into three stages [9,10]; (i) the initial stage of pipeline production, (ii) the stable operation stage of pipelines, and (iii) the aging stage of pipelines. During the first stage (which usually lasts from 0.5 to 2 years), there is a high probability of breakage, mainly due to problems in pipe material, construction, etc., with approximately five failures per 1000 km of pipeline. The process linked to the second stage usually lasts 15 to 20 years, mainly due to corrosion and third-party damages. In this case, the number of failures is approximately twice per 1000 km of pipeline. The probability of pipeline failures increases significantly during the third stage, mainly due to increased wear and corrosion.
In these scenarios, several detection techniques for pipeline leaks have been developed to reduce the loss of transported products in pipelines, even though these techniques are mainly designed for water pipelines. Overall, the different approaches can be categorized as follows (according to [10]): (i) automatic, semiautomatic, or manual methods; (e.g., [11]); (ii) direct or indirect methods (e.g., [12]); (iii) optical and nonoptical methods (e.g., [13]); (iv) hardware- or software-based approaches (e.g., [14]). Nevertheless, despite the efforts to minimize the adverse effects of pipeline failures, the protection of aquifer systems impacted by oil loss can only be guaranteed by optimizing the remediation actions in terms of effectiveness and rapidity through the effective simulation of fluid migration in the subsurface, from the pipeline failure to the groundwater.
The contamination impact and the migration of contaminants into the surrounding environment vary based on the quantity and quality of spilled oil products. Among the possible consequences, the one related to the infiltration of nonaqueous phase liquids lighter than water (LNAPL) is of utmost importance. LNAPLs, such as crude oil and several refined products, are hydrophobic liquid organic compounds infiltrating vertically through the unsaturated zone (after spills from a pipeline). A fraction of the same product is trapped (discontinuous phase) within the aquifer medium’s interconnected pores (or fractures). In contrast, another fraction (continuous phase) flows in the aquifer on top of the groundwater table. The continuous phase migrates (or can be mobilized) according to the hydraulic gradient. Some constituents (e.g., the so-called BTEX) will dissolve in groundwater to be transported downgradient regarding the pipeline breakage. However, LNAPL migration is difficult to predict (and simulate) due to its features and the hydraulic heterogeneities of the aquifer media.
The determination of spilled oil contaminant migration in the continuous phase can be described using the governing equations of immiscible phase fluid flow in a porous medium. Predicting and simulating a multiphase flow of hydrocarbon fluids in groundwater is challenging due to the nonlinear nature and coupled properties of its governing equations. Indeed, there are still many scientific challenges to accurately describe the spilled oil’s subsurface multiphase flow process. Predicting the contaminated range due to spilled oil requires computationally massive numerical solutions of multiphase flow equations written in terms of soil parameters, which in most cases are not precisely known due to spatial heterogeneities, parameter estimation errors, etc. [10]. The nonlinearities arise since the capillary pressure and permeability of each fluid phase are a function of the saturation and, due to the gravity and pressure gradients, are responsible for creating sharp (shock) fronts and rarefaction, which may introduce significant errors in the numerical simulations.
Many numerical models have been developed in the last decades to determine the migration of the oil saturation distribution in aquifers of various geometries that are subjected to different boundary water flow conditions. See, for example, a numerical model for subsurface oil pollutant migration (from a leaking underground tank or pipeline) modeled by a nonlinear convection-diffusion equation with power law nonlinearities [15], a two-dimensional oil spill numerical model for two-phase immiscible fluids in a variably saturated zone using a finite-element method [16], modeling of an immiscible-miscible nonaqueous phase of organic liquid compounds flows in the subsurface [17], and a numerical study using a modified USGS solute transport numerical model to predict high-density hydrocarbon migration [18].
Immiscible multiphase flow in porous media can be investigated using laboratory experiments and numerical investigation results of immiscible multiphase fluid flow [19] using an iterative compositional model for subsurface multiphase fluid flow [20]. Three-dimensional numerical modeling of multiphase flow and transport has been investigated in porous media [21], and in fractured media [22].
Different numerical methods have been used to resolve the partial differential equations that describe the multiphase fluid flow. A comparison between numerical formulations, in particular, the Picard and Newton–Raphson linearization approach for two-phase fluid flow in porous media, can be found in [23,24]. In [25], a three-dimensional element-free Galerkin (EFG) code for simulating two-phase fluid flow in porous materials is developed. Underwater spreading of oil leakage from damaged submarine pipelines was investigated numerically in [26], using a coupled Eulerian–Lagrangian CFD model [27] or modeling spills in ice-infested water [28]. Transport and biodegradation modeling of gasoline spills in soil-aquifer systems using a central upwind scheme for a compressible two-phase fluid flow model [29,30] or using buried steel pipelines was also investigated [31]. More recently, a high-resolution shock-capturing numerical method has been introduced to simulate three-dimensional immiscible three-phase fluid flow in variably saturated zones [32,33].
Despite the efficiency of the recent numerical models for modeling subsurface multiphase flow, analytical solutions to very simplified cases might be interesting for gaining some physical insight into the flow processes. See, for example, an oil plume distribution study in an anisotropic porous medium [34], in an aquifer near an impermeable barrier [35], and the study of the migration of a fuel pollutant [36]. They describe the nonlinear diffusion convection with a power law nonlinearity (using an approximation valid for low saturations) [34,35,36], where the solution is characterized by a shock wave moving forwards or with the inclusion of the gravity term [37]. Or the problem of an oil-water injection into a semi-infinite reservoir is described by a power-law nonlinear convection-diffusion equation [38], which includes the capillary pressure effect. It assumes that the position of the moving front of the water saturation profile is proportional to the square root of time.
Xu et al., 2015 [29] presented a gasoline spill model of a soil-aquifer system similar to the one performed in this paper. The numerical simulations for the transport and biodegradation of the gasoline spill were built up by combining the hydrocarbon spill screening model (HSSM) [39,40] for the unsaturated zone and a modified MT3DMS (modular three-dimensional multispecies transport model) [41,42] in the saturated one. The numerical model is divided into three steps: (1) the vertical migration of gasoline in the vadose zone, (2) the spread and dissolution processes on the groundwater table, and (3) the transport and biodegradation of dissolved benzene from gasoline in the saturated zone. The first two steps, (1) and (2), are simulated using HSSM, and the last step, (3), is simulated using a modified MODFLOW/MT3DMS. Although the actual flow in the vadose zone is three-dimensional, HSSM treats migration through the vadose zone as being one-dimensional downward.
The main difference between [29], which investigates a similar pipeline leak, and this paper is that (1) we use the flow equations for immiscible fluid flow applied to both vertical and horizontal zones. The variably saturated zone is treated and resolved as a whole. The system is not divided into zones and treated differently as in [29]; (2) we consider the migration of the three phases through the variably saturated zone as immiscible. The effects of volatilization/biodegradation or dissolution are not considered (while Ref. [29] considers dissolution of fluids since it uses the transport equations to determine concentrations); (3) we use the CactusHydro code that uses an HRSC flux conservative method (introduced in [32,33] to groundwater) that treats the advective part of the flow equation for each phase better than other numerical methods when there is a discontinuity since it is a conservative flux method. In our approach, the air density is constant (volume conservation corresponds to mass conservation). This simplification is correct as far as the temperature of the air phase is almost constant, and the pressure of the air is always of the order of the atmospheric pressure (no high-pressure air bubbles in the system). The computation assumptions are constant density and constant viscosity in each phase.
The new contribution of the paper is that, in the investigation of an oil spill pipeline, using the HRSC method, we can follow the saturation contours of the contaminant very precisely and investigate the behavior of the oil spill (with a pressure different from the atmospheric one and for a fixed time) as a function of several parameters.
The purpose of this work is to investigate the effects of an oil pipeline failure in a range of hydrogeological conditions by varying the unsaturated depth zone, the water saturation of the unsaturated zone, the oil contaminant density, two different values of the hydraulic gradient, and an oil pipeline pressure using the high-resolution shock-capturing flux (HRSC) conservative method and the CactusHydro code recently introduced in [32,33]. This code allows us to accurately simulate three-phase immiscible fluid flow (in this case, gasoline, diesel, or any LNAPL) in a porous medium, minimizing the time necessary to implement the subsequent remediation actions and optimizing the type, number, and technical features of the actions. The uniqueness of this study is that it reveals the oil pipeline leakage in a variably unsaturated zone and the migration characteristics under different influencing factors, some of which have not been studied before. The results of this investigation can be used as a reference point for dealing with emergencies involving oil leaks from pipelines.

2. Materials and Methods

2.1. Study Area

The study area falls in the eastern portion of central Italy, between the Apennine chain and foothills to the southwest and the Adriatic Sea coast to the northeast. From a geological point of view, it is characterized by the typical Plio-Pleistocene sedimentary succession found all over the periadriatic Marche-Abruzzo coastal zone. Several studies have already described and discussed this stratigraphic succession in portions of the eastern Apennines [43,44,45,46,47] and the Adriatic foredeep [48,49,50]. Specifically, this succession is dominated by a regressive marine depositional sequence that occurs through shallow marine to foreshore units to evolve into a coastal lagoon with alluvial channeling deposits. Compressional tectonics of the growing Apennines dominated the area during the placement of this succession.
Linking to the official cartographic nomenclature (ISPRA-CARG), the area is characterized by the Mutignano Formation (FMT), which is composed of shallow marine clayey lithologies (FMTa), transitioning upward to beach sands and sandstones with intercalation of imbricated gravels (FMTd). At the top of the succession, outcropping at ground level, the presence of the Ripa Teatina Clays and Conglomerates (RPT), deposited in a marine to the continental transition area, covers the underlying units.
Due to their lithological differences, the geological units mentioned above constitute the main hydrogeological complexes in the area. FMTd and RPT units, mainly the sandy-gravel sediments located therein, generally represent the main aquifers. Conversely, the FMTa unit is the bottom aquiclude (Figure 1).

2.2. Hydrogeological and Geological Data

The presence of technical and stratigraphic reports related to the construction of water wells in the study area is included in ISPRA’s public database (https://www.isprambiente.gov.it/it/banche-dati/banche-dati-folder/suolo-e-territorio/dati-geognostici-e-geofisici, 1 January 2023) provided a helpful dataset for the purposes of this paper. In particular, obtaining a series of stratigraphic logs derived from perforation cuttings enriched by a series of hydrogeological data, including hydraulic head measures, was possible.
Existing geological cartography at different scales also furnished valuable support in finalizing the following work.

2.3. Mathematical Model

The numerical model developed in this paper is based on a conceptual model describing the risk of leakage of an immiscible LNAPL contaminant (gasoline or diesel oil) released to the environment from an onshore oil pipeline within a heterogeneous aquifer system.
The oil pipeline is one meter from the ground surface in the unsaturated zone. The possible spill would migrate toward the saturated zone under the effect of the gravity force, and it would be released for 3600 s before being noticed and subsequently stopped by the owner of the pipe. Due to a hydraulic gradient, the immiscible contaminant would move together with the flow while it moves downward.
A three-dimensional numerical model is set up using the numerical code CactusHydro, introduced in [32,33]. CactusHydro resolves the governing equations that describe the migration of a three-phase fluid flow in a porous medium composed of nonaqueous (n), water (w), air (a), and a variably saturated zone. It is based on a high-resolution shock-capturing (HRSC) flux conservative method [52,53,54] that follows sharp discontinuities accurately and temporal dynamics of three-phase immiscible fluid in a porous medium, shows the absence of spurious oscillations in the solution and converges to the ‘weak’ solution as the grid is refined. The time evolution is performed using a forward-in-time explicit method (instead of the commonly used implicit ‘backward’ in time evolution method). The method used in HRSC flux conservative belongs to the class of Monotonic Upstream-Centered Scheme for Conservation Law (MUSCL) suggested by van Leer in 1973, and the Kurganov–Tadmor (KT) scheme [52] that is up-to-second-order accurate in space using analytical information on the form of the flux. However, as in this work, using the first-order upwind formula for the fluxes, and the minmod flux limiter, only the point values of the flux are required. That is a great advantage since it can use tabulated values for permeabilities. Consequently, the integration of time is performed in first order (for consistency) using the Euler method.
CactusHydro is based on the Cactus computational toolkit [55,56,57], an open-source software framework for developing parallel high-performance-computing (HPC) simulation codes, and data is evolved on a Cartesian mesh using Carpet [58,59]. The migration of the spilled LNAPL contaminant is considered immiscible. The effect of volatilization, biodegradation, or dissolution is 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 equation includes both zones).
In this paper, we introduce an initial condition that mimics a fixed-time oil spill from a broken pipeline to predict oil spill trajectories and precisely follow the migration of the immiscible contaminant, varying the distance of the oil pipeline to the groundwater table, the saturation of the unsaturated zone, oil contaminant density, and analyzing two different values of the hydraulic gradient (0.04 and 0.004). To simplify the three-dimensional numerical simulations, we assume the variably saturated zone’s grid geometry to be a parallelepiped of 160 m long with an impermeable oil pipeline along this direction situated at one-meter depth from the land surface, 80 m wide, a variably unsaturated zone ranging between 1 and 20 m, and a saturated zone of 16 m. The leakage is represented by a cube (into the pipeline) containing a volume of immiscible contaminant. It comes out with high pressure at the atmospheric pressure of the unsaturated zone, approximately 2.0 × 10 6   Pa for 3600 s. The range of the unsaturated zone thickness was selected to simulate the possible impact on both (i) the shallower perched temporary groundwater and (ii) the deeper perennial groundwater.

2.4. Hydrogeological and Hydrocarbon Phase Parameters

First, it investigated two types of LNAPL contaminants, gasoline and diesel oil. Table 1 shows the hydrocarbon phase properties, such as density, viscosity, and parameter details used in the numerical simulations for the gasoline leak. In particular, the density of the gasoline is 750   kg / m 3 , which corresponds to an LNAPL being less dense than water.
For a three-phase fluid flow, we need two capillary pressures. In our formulation, we used 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 . The nonaqueous-water capillary pressure is then given by, p c n w = p n p w = p c a w p c a n . 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 [32,33,60],
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 ,
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, it uses the van Genuchten model [61] 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
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 rock compressibility c R is taken from [62], considering a porous medium to be sand. The porous medium’s effective porosity (at a pressure of one atmosphere) is 0.43, corresponding to clean sandy/silty sandy [63]. We consider a constant and isotropic hydraulic conductivity, since that was the value measured in the study area to be K = 6.8 × 10 5   m / s . From here, we determined the value of the absolute permeability ( k = K   μ w ρ w   g ), see Table 1, which guarantees a contaminant’s leak of approximately 2.78 × 10 4   kg in 3600 s at a pressure of approximately 2.0 × 10 6   Pa . Indeed, the absolute permeability is k = 2.059 × 10 11   m 2 all over the grid with k x = k y = k z . The numerical simulations were implemented using the van Genuchten model and an irreducible wetting phase saturation of 0.045, which is compatible with porous sand [64]. The van Genuchten α parameter is defined as:
α = p c ρ w g 1 ,
where p c is the capillary pressure head. For a clean sand/silty sand, the value is given by α = 0.145   cm 1 = 14.5   m 1 . From Equation (6), we can determine the capillary pressure at zero saturation between air and water, which gives:
p c a w = ρ w g α = 10 3   kg / m 3 9.81   m / s 2 14.5   m 1 = 676.55   Pa
Using the value of the superficial tension σ n w = 2.6 × 10 2   N / m and σ a w = 6.5 × 10 2   N / m we have that β n w = σ a w σ n w = 6.5 × 10 2   N / m 2.6 × 10 2   N / m = 2.5 , and p c n w S w = p c a w β n w = 270.62   Pa . Then, p c a w = p c a n + p c n w . From here, we get the capillary pressure at zero saturation, p c a n = p c a w p c n w = 676.55 270.62 = 405.93   Pa .
Table 2 shows the hydrocarbon phase properties, such as density, viscosity, and parameter details used in the numerical simulations for the diesel oil leakage. In particular, the density of the diesel oil is 830   kg / m 3 , which also corresponds to an LNAPL being less dense than water.
Table 1 and Table 2 differ in the value of the immiscible contaminant released. As for the capillary pressure, following the same procedure used for the gasoline case and using the value of the superficial tension of the diesel oil, σ n w = 3.0 × 10 2   N / m and σ a w = 6.5 × 10 2   N / m . Indeed, we have that β n w = σ a w σ n w = 6.5 × 10 2   N / m 3.0 × 10 2   N / m = 2.17 , and p c n w S w = p c a w β n w = 311.77   Pa . Since p c a w = p c a n + p c n w . From here, we get the value of the capillary pressure at zero saturation, p c a n = p c a w p c n w = 676.55 311.77 = 374.68   Pa .

3. Results

3.1. Hydrogeological Conceptual Model

In all stratigraphic logs, continuous clays typical of a marginal marine environment (FMTa) are noted at the base of the succession, transitioning to shoreline sands and gravels (FMTd). Gravels and clays associated with lagoon areas with fluvial channel incisions (RPT) are superimposed on FMTd with a discordant contact. Geometric evidence of geologic contacts, as well as local geomorphological features, suggested the presence of two faults within the studied aquifer system (Figure 2). The aquifer bottom was found to be the contact between the sands and gravels of FMTd and the underlying clays of FMTa. Inside the FMTd and RPT units, several clay strata were found to generate temporary perched aquifers (Figure 2), in agreement with findings in other heterogeneous systems worldwide (e.g., [64]). The western fault, associated with the local rise of the clayey bottom, defines a local groundwater divide. The deeper groundwater generally flows with a few units of hydraulic gradient while abruptly increasing to approximately 10% within the eastern fault (Figure 2). This local increase in hydraulic gradient is clearly associated with a low-permeability fault core [65,66,67,68,69,70,71,72,73,74] that acts as a barrier to groundwater flow, therefore causing the aquifer medium to be a sort of basin-in-series system (sensu [66]).

3.2. Three-Dimensional Numerical Simulations, Results, and Discussions

This section presents twelve results out of twenty-eight (see Supplementary Materials) different three-dimensional numerical simulations of immiscible leaks of LNAPL (gasoline or diesel oil) that migrate in a variably saturated zone from an oil pipeline. The migration is investigated by varying: the unsaturated depth zone, the hydraulic gradient, and the water saturation of the unsaturated zone starting from dry soil porous medium. It is considered a three-phase fluid model and follows the temporal evolution of the contaminant migration following the saturation profiles of each phase along the variably saturated zone. We also determine the time for the LNAPL contaminant to arrive at the groundwater table from the oil pipeline and the leakage location. Different situations are analyzed in which the distance from the contaminant to the groundwater table varies from 1 to 2, 5, 10, and 20 m, with two different values of the hydraulic gradient (0.04 and 0.004) and different values of water saturation in the unsaturated zone: 0 (dry soil), 0.20, and 0.50. In particular, the results for unsaturated depth zones equal to 2, 5, and 20 m can be found in the Supplementary Materials.

3.2.1. Numerical Simulations of a Gasoline Leak from an Oil Pipeline in a Dry Zone

First, it is considered a leak of gasoline released into the environment from an oil pipeline situated in the unsaturated zone at a one-meter depth from the land surface for a fixed time of 3600 s (a time in which the pipeline owner realizes an emergency before deciding to close the pipeline). The immiscible contaminant moves inside the oil pipeline with an experimental value of the pressure that may vary between 2.0 2.4 × 10 6   Pa , while the rest of the unsaturated zone is at the atmospheric pressure. After the contaminant leakage, it migrates from the unsaturated zone toward the saturated one under the effect of the gravity force. For the transient numerical simulations, we consider a three-phase fluid model composed of water, air, and a fixed quantity of gasoline. We follow the contaminants released through the variably saturated zone. Figure 3 shows one of the grid geometries used in this paper. We assume the variably saturated zone grid geometry to be a parallelepiped of 160   m long from x = 100 , + 60   m (left-hand side), 80   m wide from y = 40 , + 40   m (right-hand side), and 18   m depth from z = + 2 , 16   m . We consider (otherwise mentioned) a spatial resolution of dx = dy   =   dz   = 0.50   m , and a time step resolution of d t = 0.025   s . The leakage of contaminants is situated at x , y , z = 0 , 0 , 1   m . The porous medium is composed of an unsaturated dry zone (air-NAPL) and a saturated one (where initially there is only water), separated by the groundwater table. The groundwater table crosses the point x , y , z = 0 , 0 , 0   m . All boundary conditions are no-flow, except in the infiltration zone on top of the parallelepiped. The bottom part of the parallelepiped is impermeable, with an absolute permeability of 2.059 × 10 15   m 2 (ten thousand times smaller than the value of the rest of the domain, see Table 1 for details). The dimension of the grid has been chosen to cover the minimum possible part of the aquifer. However, it is big enough to avoid the boundary effect (finite grid size) on the dynamics of the contaminant.
After being released, from the oil pipeline, into the environment, the LNAPL migrates downward into the unsaturated zone under the influence of gravity. The contaminant is released in all directions in the Cartesian coordinate system, assuming no preferential direction, with a pressure equal to 2.0 × 10 6   Pa . Figure 4 shows the three-dimensional numerical simulation results of the saturation contours, σ n = S n ϕ (do not confuse with the interfacial tension sigma) at different times, for the three-phase immiscible fluid flow (water, gasoline and air) in a dry unsaturated zone, using a spatial grid resolution of 0.50   m and a grid dimension of 160 × 80 × 18   m . The saturation contours are viewed on the left-hand side on the z x plane. On the right-hand side, on the z y plane. The green rectangle zone in each panel is amplified in the upper area at each time indicated.
The contaminant generally moves predominantly downward in the unsaturated zone due to the gravity force. Lateral spreading may also appear due to the presence of capillary pressures and the high pressure of the contaminant itself. Figure 4 shows the initial time, t = 0   s (Figure 4a), and the position of the contaminant. At 204 s (Figure 4b), the contaminant has arrived at the groundwater table. Although a consistent quantity remains in the capillary fringe, some material enters the saturated zone due to the elevated contaminant pressure. The third row (Figure 4c) shows the contaminant at one day and 4.4 h. Notice that the leak spill from the pipeline has already stopped. The fourth row (Figure 4d) shows the fate of the immiscible contaminant, which has already arrived and remained, in the capillary fringe of an LNAPL, moving into the left (negative) direction of the x-axis together with the water (see the left-hand side) due to a hydraulic gradient of 0.04. (See the blue arrows directed to the left-hand side). In contrast, the right-hand side remains symmetric around the plane z     y . After 24 days and 2.0 h, the contaminant reaches position x = 62   m .
Figure 5 shows three-dimensional numerical simulation results of the depth as a function of the water saturation S w (blue points), gasoline saturation S n (red points), and air saturation S a (green points) at various times for a gasoline leak shown in Figure 4. Initially, at t = 0   s (left-hand side), there is a sharp front of immiscible contaminant located at z = 1.0   m (red points) with a saturation S n = 1 . Below z = 1.0   m , this saturation goes abruptly to zero since initially the contaminant is situated only on top of the parallelepiped. See also Figure 4, the first row. Instead, the water saturation (blue points) is zero in the unsaturated zone (being completely dry). It is equal to one in the saturated zone, initially filled up with water. The air saturation (green points) equals one in the unsaturated region with no contaminants. Notice that the sum of the three saturations is always one at any depth value, as stated in Equation (6). After 5.7 h, the immiscible gasoline has already reached the saturated zone. Being an LNAPL, it remains around the capillary fringe. Since it reaches the aquifer with high pressure, it enters a bit into the saturated zone (center, see the red point that reaches almost a depth of 5 m. Compared with Figure 4, the third row). After two days and 8.9 h, the contaminant tends to move in the left direction due to a hydraulic gradient. Additionally, the contaminant that initially enters the saturated zone returns to the capillary fringe and stays in this zone (right-hand side).
Figure 6 shows the three-dimensional numerical simulation results of the saturation contours, σ n = S n ϕ , at different times, for the three-phase immiscible fluid flow (water, gasoline, and air) using a grid resolution of 0.50   m and a grid dimension of 160   m long from x = 100 , + 60   m (left-hand side), 80   m wide from y = 40 , + 40   m (right-hand side), and 27   m depth from z = + 11 , 16   m . The leakage contaminant is initially situated at x , y , z = 0 , 0 , 10   m in the oil pipeline. After one hour (at 3686 s), the immiscible contaminant still has not reached the groundwater table (Figure 6a). At 6553 s, the contaminant arrives at the groundwater table and enters the saturated zone (see (Figure 6b) at 7168 s) but remains in the capillary fringe zone, being an LNAPL. Then it moves in the left direction due to a hydraulic gradient equal to 0.04 (Figure 6c). After 12 days and 1.2 h (Figure 6d), the contaminant reaches position x = 41.0   m (Figure 6d). Notice that the contaminant saturation on top of the grid is decreasing since the leakage in the pipeline is closed after one hour. The number of cells used for this numerical simulation is 160 × 80 × 27 × 8 = 2 , 764 , 800 .
Figure 7 shows similar three-dimensional numerical simulation results to Figure 4, but the hydraulic gradient is 0.004 instead of 0.04 (see Figure S1 in the Supplementary Materials for details). Figure 7 shows this migration on the y   x plane. It can be noticed that, after 21 days and 10.3 h, the contaminant has reached position x = 20   m (on the left-hand side of the panel, the right-hand side is its amplification). The hydraulic gradient has a key role in the movement of the contaminant in the direction of the groundwater flow and has the consequence that the migration in the flow direction is lower when the hydraulic gradient is smaller.

3.2.2. Numerical Simulations of a Gasoline Leak from an Oil Pipeline in an Unsaturated Zone

The effects of the water saturation in the unsaturated zone are investigated using the same situation as before, but with a porous soil unsaturated zone that is not completely dry but has nonzero water saturation. Consider two different cases, S w = 0.20 and S w = 0.50 . This situation can be imagined as continuous rain, partially filling the unsaturated zone with water. The grid used here is like the one shown in Figure 3 (for one meter). The boundary conditions are similar to the previous cases. Except that, at each point of the grid in the unsaturated zone, S w is different from zero, and at later times, on top of the grid, this value is kept constant. In this way, the water saturation of the unsaturated zone remains constant.
Figure 8 shows the numerical simulation results similar to Figure 4, except that the unsaturated zone now has S w = 0.20 instead of S w = 0.0 (dry soil). The first difference from Figure 4 is that there is a vertical water flow in the unsaturated zone (see the blue lines) due to the gravity force (besides the one going in the left direction in the saturated zone). Initially, the contaminant is released on top of the grid, and at 204 s already arrived at the groundwater table (the second row) in Figure 8. In this case, the presence of water in the unsaturated zone does not substantially change the situation in the dry case (see the third row at one day and 4.4 h, in both figures).
Figure 9 shows a similar case as in Figure 8, but the water saturation in the saturated zone is increased to S w = 0.50 . In this case, the contaminant reaches the groundwater table before 204 s (Figure 9b), which means there is no change from the previous case due to the high pressure of the contaminant leakage. After one day and 4.4 h, the groundwater table level increases due to the filling of the voids in the porous medium. See (Figure 9b) and (Figure 9c), respectively.
Similar numerical simulations were performed, varying the depth of the unsaturated zone. Figure 10 shows the saturation contours σ n = S n ϕ , at different times. The source of the contaminant is initially situated at x , y , z = 0 , 0 , 10   m in the oil pipeline, and the saturation of the water in the unsaturated zone is S w = 0.20 . Figure 11, instead, shows a similar scenario to Figure 10, but the saturation of the water in the unsaturated zone is increased to S w = 0.50 . Notice the presence of vertical water flow (see the blue lines) compared to Figure 6, which corresponds to a dry unsaturated zone case.
The presence of water saturation in the unsaturated zone accelerates the arrival time of the contaminant into the aquifer. In the case of Figure 10, the contaminant arrives at the groundwater table at 4915.0 s (see (Figure 10b)), while in Figure 11, it arrives at 4710.4 s (see (Figure 11b)). The more water content in the unsaturated zone, the faster the contaminant reaches the saturated zone. Additionally, water in the unsaturated zone inhibits the spill of contaminants from the pipeline compared to the dry soil case (Figure 6). Indeed, the quantity of contaminant spilled after one hour is 213,904 kg in Figure 6, 197,943 kg in Figure 10, and 158,647 kg in Figure 11. The saturation contour values of σ n are lower in Figure 11 than in the case of Figure 10.

3.2.3. Numerical Simulations of a Diesel Oil Leak from an Oil Pipeline in a Dry Zone

Consider now a diesel oil leak (denser than gasoline) released into the environment from the oil pipeline in the unsaturated zone, just as in the previous cases. The parameters of the numerical simulations are listed in Table 2. Figure 12 shows similar numerical simulation results as in Figure 4 (dry soil) but using diesel oil instead of gasoline for the saturation contours, σ n = S n ϕ , at different times, for the three-phase immiscible fluid flow (water, diesel oil, and air) using a grid resolution of 0.50   m and a grid dimension of 160   m long from x = 100 , + 60   m (left-hand side), 80   m wide from y = 40 , + 40   m (right-hand side), and 18   m depth from z = + 2 , 16 m . The leakage contaminant is initially located at x , y , z = 0 , 0 , 1   m in the oil pipeline (see (Figure 12a)). After 204 s (Figure 12b), the contaminant arrived at the groundwater table. Once in the groundwater table, it moves to the left-hand side due to a hydraulic gradient of 0.04 (Figure 12c). After nine days and 0.8 h, the diesel oil contaminant reaches position x = 12.5   m (see left-hand side of (Figure 12d)), in contrast with Figure 4 in which, after nine days and 11.6 h, the gasoline contaminant reaches position x = 39.5   m .
Figure 13 shows numerical simulation results similar to those in Figure 6 (dry soil) but uses diesel oil instead of gasoline. In Figure 13, the contaminant arrives at the groundwater table long after the gasoline case of Figure 6, say, 3.081 days. In general, it takes longer for diesel oil to arrive at the groundwater table with respect to gasoline and has less mobility to move together with the groundwater flow generated by a hydraulic gradient equal to 0.04.
See the Supplementary Materials for a similar case with a hydraulic gradient equal to 0.004 instead of 0.04.

3.2.4. Numerical Simulations of a Diesel Oil Leak from an Oil Pipeline in an Unsaturated Zone

The effects of water saturation in the unsaturated zone are also investigated in the case of diesel oil. Consider the two different cases, S w = 0.20 and S w = 0.50 , in the unsaturated zone. Figure 14 shows the numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.20 . After 204 s, the contaminant has already arrived at the groundwater table, as in the case of Figure 12. Notice how the contaminant enters the saturated zone (due to elevated pressure from the pipeline). Then the contaminant starts to float and moves along with the groundwater flow (see the third row). After five days and 14 h, the diesel oil reaches position x = 8.5   m to be compared with the x = 9.5   m in the case of a diesel oil spill in a dry unsaturated zone (see Figure 12) after five days and 22 h.
Figure 15 shows the same situation as Figure 14, but with S w = 0.50 . As can be seen, the contaminant arrived at the groundwater table at 204 s (Figure 15a) since the effects on the oil pipeline pressure are much stronger than the effects on the saturation of the unsaturated zone. After five days and 14.0 h (Figure 15d), the diesel oil contaminant does not move much compared to the previous case, with water saturation equal to 0.20 and dry soil. Notice that the level of the groundwater table increases (Figure 15d) due to the unsaturated zone being filled with water.
Figure 16 and Figure 17 show the same situation as Figure 14 and Figure 15, respectively, except that the diesel oil spill is now released at 10 m instead of one meter from the oil pipeline. The diesel oil contaminant in the case of S w = 0.20 , Figure 16, arrives at the groundwater table at 1.754 days, while the water saturation is S w = 0.50 ,  Figure 17, the contaminant arrives at the groundwater table at 0.853 days. Both results are to be compared with a similar case but in a dry unsaturated zone (Figure 6) in which the diesel oil arrives at the groundwater table at 3.319 days. See the following subsections for details.
Further numerical simulation results on the saturation contours of a multiphase immiscible fluid flow using gasoline or diesel oil for 2 m, 5 m, 20 m, and a hydraulic gradient of 0.04 and 0.004 are reported in the Supplementary Materials.

3.3. Effects on the Density of the Contaminant

Two species of oil (gasoline and diesel oil) with densities of 750   kg / m 3 and 850   kg / m 3 , respectively, are considered for analyzing the influence of oil density on oil spreading in the case of an oil pipeline leaking and dry soil. Table 3 reported the time employed for the LNAPL contaminant in a numerical simulation to arrive at the groundwater table and the position in the x -coordinate at one day and 4.4 h.
Figure 18 shows the depth of the oil pipeline as a function of the arrival time at the groundwater table in unsaturated dry soil for two different species of oil (gasoline, red line, and diesel oil, green line) and a hydraulic gradient of 0.04. As can be seen from this figure, there are no substantial differences between both species for depths below 2.0 m since the oil spills have a very high pressure with respect to the atmospheric pressure of the unsaturated zone. For this reason, other effects are not noticeable. For depths greater than 5 m, it takes diesel oil to arrive at the groundwater table much longer than gasoline, although it is much denser. This is because diesel oil’s viscosity exceeds gasoline’s (see Table 1 and Table 2). In general, diesel oil spreads less than gasoline once it arrives at the groundwater table. Table 3 also shows that the arrival time to the groundwater table does not depend on the value of the hydraulic gradient in the saturated zone, as it should. But the hydraulic gradient affects the displacement in the x direction. See, for example, the last column of Table 3, which clearly shows that the displacement is smaller when the hydraulic gradient equals 0.004.

3.4. Effects on the Water Saturation of the Unsaturated Zone

The effects on the water saturation of the unsaturated zone are reported in Table 4. The time employed for the LNAPL contaminant in a numerical simulation to arrive at the groundwater table and the position in the x coordinate at one day and 4.4 h are investigated as a function of the depth of the unsaturated zone and the two different values of the water saturation (hydraulic gradient of 0.04), as shown in Figure 8, Figure 9, Figure 10 and Figure 11 and Figure 14, Figure 15, Figure 16 and Figure 17. Figure 19 shows the water saturation of the unsaturated zone as a function of the arrival time at the groundwater table for both gasoline (left-hand side) and diesel oil (right-hand side), as reported in Table 4. The data was taken from the numerical simulations. The arrival time at the groundwater table (for a fixed depth) increases as the water saturation decreases for LNAPL densities, although the arrival time for the diesel oil case is much longer than the gasoline case since the viscosity of the diesel oil is an order of magnitude higher than the gasoline one (the right-hand side). The water saturation of the unsaturated zone influences the arrival time. As the water saturation increases, the arrival time decreases (for a fixed depth). Table 3 and Table 4 show the effects on the thickness of the unsaturated zone (the oil pipeline is situated at one-meter depth in the unsaturated zone, and the unsaturated zones vary from 1.0 up to 20.0 m. The more profound the unsaturated zone is, the more time the contaminant takes to arrive at the groundwater table.

3.5. Effects on Pressure in the Oil Pipeline

All the simulations previously shown consider an oil pipeline pressure of 2.0265 × 10 6   Pa during the time of the leakage, which was fixed to be one hour, and an atmospheric pressure value in the unsaturated zone as the initial condition. The high pressure inside the oil pipeline has the most critical effect on the arrival time of the contaminant at the groundwater table when the distance between the leakage and the groundwater table is smaller than two meters, as was previously shown. Consider now the case where the pressure inside the oil pipeline is also atmospheric. Table 5 shows the numerical simulation data results for the arrival time for the gasoline spill with atmospheric pressure, different values for the water saturation in the unsaturated zone, and a hydraulic gradient of 0.04.
Figure 20 shows the arrival time of the gasoline spill at the groundwater table as a function of the water saturation in the unsaturated zone at a one-meter depth (see Table 5) and the case in which the contaminant spill has an atmospheric pressure (green line) and another case in which the spill pressure is 2.0 × 10 6   Pa   (red line, previous cases). As can be seen, in the case of high pressure, the water saturation of the unsaturated zone does not affect the arrival time very much since the depth is not much longer (red line), and the effect of the high pressure is stronger. In the case of a gasoline spill with atmospheric pressure, the arrival time to the groundwater table decreases when the water saturation of the unsaturated zone increases.

3.6. Validation of the Results

In a previous publication [32,33], the numerical method (HRSC) and the code were validated using several toy models and a real experiment, such as a sandbox. Since here, it was not possible to validate the numerical results to experimental ones, we had to rely on classical convergence tests running the same code at different resolutions. Figure 21 shows comparison results on the saturation contours for a gasoline contaminant spill from an oil pipeline at two different times using a grid resolution of 0.25 m and a time step size of 0.00625 s (first and third row) and a grid resolution of 0.50 m and a time step size of 0.0025 s (second and fourth row, see Figure 4). Both results are similar and show that using two different resolutions, the results are similar within a 5% error.

4. Conclusions

In this work, we presented a three-dimensional numerical modeling investigation on the effects of an oil pipeline failure in a range of hydrogeological conditions, using the high-resolution shock-capturing (HRSC) flux conservative method, and the CactusHydro code, recently introduced in [32,33]. We investigated two different density types (gasoline and diesel oil) in a variably saturated zone and the effects on the variation of the depth of the unsaturated zone, three different saturations of the unsaturated zone that include the dry soil case, the hydraulic gradient and the pressure of the oil pipeline leakage. A model limitation is that at this stage, we consider these fluids immiscible, and the processes of volatilization and/or dissolution is not considered. We investigate the contour saturation of the contaminant and its migration from the oil pipeline leakage with great accuracy. The results indicate that the leaking oil’s pressure is the parameter that most significantly affects the contaminant’s arrival at the groundwater table. Additionally, the water saturation of the unsaturated zone influences the arrival time at the groundwater table. Indeed, the arrival time decreases (for a fixed depth) when the water saturation of the unsaturated zone increases. This is because the contaminant moves faster (easiest downward) when the saturation increases. This behavior was also observed experimentally [75]. The unsaturated depth zone significantly influences the contaminant migration in the saturated zone/capillary fringe, although it results in no significant change in the vertical direction. The unsaturated zone depth significantly influences the contaminant migration in the unsaturated zone. At the same time, the oil density and the hydraulic gradient have limited effects on the contaminant migration in the variably saturated zone. In particular, the gasoline case moves faster than the diesel oil case (for a fixed depth) since the viscosity of the diesel oil is an order of magnitude higher than the gasoline one.
In the application field, CactusHydro Code represents an advanced tool to be applied in the environmental risk assessment caused by hydrocarbon releases by onshore pipelines, in accordance with the provisions of the Legislative Decree 26 June 2015, n. 105 “Implementation of Directive 2012/18/EU relating to the control of the danger of major accidents connected with dangerous substances” (the so-called Seveso III Directive). The assessment of the aquifer’s vulnerability and the analysis of the environmental consequences of a major accident are specific components of the safety reports that must be prepared by the plant managers as required by the aforementioned legislation. With this in mind, the results of applying the CactusHydro code can be used as a support of utmost importance for identifying the most effective actions aimed at preventing and/or reducing the probability and extent of “pollution” and damage to environmental receptors in the event of a major accident. In conclusion, it should be emphasized that the CactusHydro code can also be used at 360 degrees to evaluate the effects and environmental impacts deriving from leaks of hydrocarbons and toxic substances, i.e., dense nonaqueous phase liquid from underground or above-ground tanks.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/w15101900/s1.

Author Contributions

Conceptualization, A.F., R.P., E.S. and F.C.; methodology, A.F., R.P., E.S. and F.C.; software, A.F.; validation, A.F.; writing—review and editing, A.F., R.P., E.S. and F.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data is contained within the article and the Supplementary Materials.

Acknowledgments

This work used resources from the University of Parma at (https://www.hpc.unipr.it, 1 January 2023). 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 three anonymous reviewers for their interesting comments/questions.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Li, X.; Chen, G.; Zhu, H. Quantitative risk analysis on leakage failure of submarine oil and gas pipelines using Bayesian network. Process Saf. Environ. Prot. 2016, 103, 163–173. [Google Scholar] [CrossRef]
  2. Xinhong, L.; Guoming, C.; Renren, Z.; Hongwei, Z.; Jianmin, F. Simulation and assessment of underwater gas release and dispersion from subsea gas pipelines leak. Process Saf. Environ. Prot. 2018, 119, 46–57. [Google Scholar] [CrossRef]
  3. Abbas, M.; Jardani, A.; Ahmed, A.S.; Revil, A.; Brigaud, L.; Bégassat, P.; Dupont, J.P. Redox potential distribution of an organic-rich contaminated site obtained by the inversion of self-potential data. J. Hydrol. 2017, 554, 111–127. [Google Scholar] [CrossRef]
  4. Gainer, A.; Cousins, M.; Hogan, N.; Siciliano, S.D. Petroleum hydrocarbon mixture toxicity and a trait-based approach to soil invertebrate species for site-specific risk assessments. Environ. Toxicol. Chem. 2018, 37, 2222–2234. [Google Scholar] [CrossRef]
  5. Yan, Y.; Xiong, G.; Zhou, J.; Wang, R.; Huang, W.; Yang, M.; Wang, R.; Geng, D. A whole process risk management system for the monitoring and early warning of slope hazards affecting gas and oil pipelines. Front. Earth Sci. 2022, 9, 812527. [Google Scholar] [CrossRef]
  6. Ning, P.; Jiang, Y.-J.; Tang, J.-J.; Xie, Q.-J. Research on the Early Warning Model for Pipelines Due to Landslide Geohazards under Multiple Influencing Factors. Water 2023, 15, 693. [Google Scholar] [CrossRef]
  7. Nesic, S. Key issues related to modelling of internal corrosion of oil and gas pipelines—A review. Corros. Sci. 2007, 49, 4308–4338. [Google Scholar] [CrossRef]
  8. Papavinasam, S. Corrosion Control in the Oil and Gas Industry; Elsevier: Amsterdam, The Netherlands, 2013; ISBN 97801239-73061. [Google Scholar]
  9. Farshad, M. Plastic Pipe Systems: Failure Investigation and Diagnosis; Elsevier: Oxford, UK, 2011. [Google Scholar]
  10. Liu, S.; Wang, Y.; Liang, Y. Environmental consequences analysis of oil spills from onshore pipelines with parametric uncertainty. Process Saf. Environ. Prot. 2020, 14, 123–134. [Google Scholar] [CrossRef]
  11. Murvay, P.S.; Silea, I. A survey on gas leak detection and localization techniques. J. Loss Prev. Process Ind. 2012, 25, 966–973. [Google Scholar] [CrossRef]
  12. Folga, S.M. Natural Gas Pipeline Technology Overview; Argonne National Lab.: Argonne, IL, USA, 2007. [Google Scholar]
  13. Batzias, F.; Siontorou, C.; Spanidis, P.M. Designing a reliable leak bio-detection system for natural gas pipelines. J. Hazard. Mater. 2011, 186, 35–58. [Google Scholar] [CrossRef]
  14. Qin, B.; Yunping, Z.; Min, F.; Xiaojian, S. Leakage detection technology of oil and gas transmission pipelines and its development trend. Petrol. Eng. Construct. 2007, 33, 19–23. [Google Scholar]
  15. Hayek, M. A model for Subsurface oil pollutant migration. Transp. Porous Med. 2017, 120, 373–393. [Google Scholar] [CrossRef]
  16. Hochmuth, D.P.; Sunada, D.K. Ground-Water Model of Two-Phase Immiscible Flow in Coarse Material. Groundwater 1985, 23, 617–626. [Google Scholar] [CrossRef]
  17. Pinder, G.F.; Abriola, L.M. On the simulation of nonaqueous phase organic compounds in the subsurface. Water Resour. Res. 1986, 22, 109S–119S. [Google Scholar] [CrossRef]
  18. Hossain, M.A.; Corapeilglu, M.Y. Modifying the USGS Solute Transport Computer Model to Predict High-Density Hydrocarbon Migration. Groundwater 1988, 6, 717–723. [Google Scholar] [CrossRef]
  19. Høst-Madsen, J.; Høgh-Jensen, K. Laboratory and numerical investigations of immiscible multiphase flow in soil. J. Hydrol. 1992, 135, 13–52. [Google Scholar] [CrossRef]
  20. Reeves, H.W.; Abriola, L.M. An iterative compositional model for subsurface multiphase flow. J. Contam. Hydrol. 1994, 15, 249–276. [Google Scholar] [CrossRef]
  21. Lagendijk, V.; Forkel, C.; Köngeter, J.; Braxein, A. Three-dimensional numerical modeling of multiphase flow and transport. J. Environ. Sci. Health Part A 2001, 36, 1473–1489. [Google Scholar] [CrossRef] [PubMed]
  22. Hoteit, H.; Firoozabadi, A. An efficient numerical model for incompressible two-phase flow in fractured media. Adv. Water Resour. 2008, 31, 891–905. [Google Scholar] [CrossRef]
  23. Ataie-Ashtiani, B.; Raeesi-Ardekani, D. Comparison of Numerical Formulations for Two-phase Flow in Porous Media. Geotech. Geol. Eng. 2010, 28, 373–389. [Google Scholar] [CrossRef]
  24. Park, C.-H.; Böttcher, N.; Wang, W.; Kolditz, O. Are upwind techniques in multi-phase flow models necessary? J. Comput. Phys. 2011, 230, 8304–8312. [Google Scholar] [CrossRef]
  25. Samimi, S.; Pak, A. A novel three-dimensional element free Galerkin (EFG) code for simulating two-phase fluid flow in porous materials. Eng. Anal. Bound. Elem. 2014, 39, 53–63. [Google Scholar] [CrossRef]
  26. Sun, Y.; Cao, X.; Liang, F. Investigation on underwater spreading characteristics and migration law of oil leakage from damaged submarine pipelines. Process Saf. Environ. Prot. 2019, 127, 329–347. [Google Scholar] [CrossRef]
  27. Sun, Y.; Cao, X.; Liang, F.; Bian, J. Investigation on underwater gas leakage and dispersion behaviors based on coupled Eulerian-Lagrangian CFD model. Process Saf. Environ. Prot. 2020, 136, 268–279. [Google Scholar] [CrossRef]
  28. Yang, M.; Khan, F.; Garaniya, V.; Chai, S. Multimedia fate modeling of oil spills in ice-infested waters: An exploration of the feasibility of fugacity-based approach. Process Saf. Environ. Prot. 2015, 93, 206–217. [Google Scholar] [CrossRef]
  29. Xu, Z.; Chai, J.; Wu, Y.; Qin, R. Transport and biodegradation modeling of gasoline spills in soil-aquifer system. Environ. Earth Sci. 2015, 74, 2871–2882. [Google Scholar] [CrossRef]
  30. Ahmed, M.; Saleem, M.R.; Zia, S.; Qamar, S. Central Upwind Scheme for a Compressible Two-Phase Flow Model. PLoS ONE 2015, 10, e0126273. [Google Scholar] [CrossRef]
  31. Zhang, J.; Liang, Z.; Han, C.J. Numerical Modeling of Mechanical Behavior for Buried Steel Pipelines Crossing Subsidence Strata. PLoS ONE 2015, 10, e0130459. [Google Scholar] [CrossRef]
  32. 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] [PubMed]
  33. 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]
  34. Pistiner, A. Oil plume distribution in an anisotropic porous medium. Transp. Porous Med. 2007, 70, 293–304. [Google Scholar] [CrossRef]
  35. Pistiner, A. Oil plume distribution in an aquifer near an impermeable barrier. Transp. Porous Med. 2009, 76, 67–75. [Google Scholar] [CrossRef]
  36. Pistiner, A.; Shapiro, M.; Rubin, H. Analysis of fuel pollutant migration in water flow through porous media. Int. J. Multiph. Flow 1989, 15, 135–154. [Google Scholar] [CrossRef]
  37. Pistiner, A.; Shapiro, M.; Rubin, H. Gravitational migration of fuel in porous media. Transp. Porous Med. 1992, 9, 187–205. [Google Scholar] [CrossRef]
  38. Chen, Z.X. Some invariant solutions to two-phase fluid displacement problem including capillary effect. SPE Reserv. Eng. 1988, 28, 691–700. [Google Scholar] [CrossRef]
  39. Weaver, J.W.; Charbeneau, R.J.; Tauxe, J.D.; Lien, B.K.; Provost, J.B. The Hydrocarbon Spill Screening Model (HSSM) Volume 1: User’s Guide; The United States Environmental Protection Agency: Washington, DC, USA, 1995.
  40. Charbeneau, R.J.; Weaver, J.W.; Lien, B.K. The Hydrocarbon Spill Screening Model (Hssm) Volume 2: Theoretical Background and Source Codes; The United States Environmental Protection Agency: Washington, DC, USA, 1995.
  41. Zheng, C.M. Recent developments and future directions for MT3DMS and related transport codes. Ground Water 2009, 47, 620–625. [Google Scholar] [CrossRef]
  42. Zheng, C.M. MT3DMS v53 Supplemental User’s Guide; Department of Geological, Sciences University of Alabama: Tuscaloosa, AL, USA, 2010. [Google Scholar]
  43. Ori, G.G.; Serafini, G.; Visentin, C.; Ricci Lucchi, F.; Casnedi, R.; Colalongo, M.; Mosna, S. The Pliocenee-Pleistocene foredeep (Marche and Abruzzo, Italy): An integrated approach to surface and subsurface geology. In Proceedings of the 3rd EAPG Conference, Adriatic Foredeep Fieldtrip Guidebook, Firenze, Italy, 26–30 May 1991; p. 85. [Google Scholar]
  44. Artoni, A.; Casero, P. Sequential balancing of growth structures, the late Tertiary example from the central Apennines. Bull. Société Géologique Fr. 1997, 168, 35–49. [Google Scholar]
  45. Carruba, S.; Casnedi, R.; Cesare, R.; Perotti, T.; Tornaghi, M.; Bolis, G. Tectonic and sedimentary evolution of the Lower Pliocene Periadriatic foredeep in Central Italy. Int. J. Earth Sci. 2006, 95, 665–683. [Google Scholar] [CrossRef]
  46. Centamore, E.; Nisio, S. Significative events in the Peryadriatic feredeeps evolution (Abruzzo e Italy). Studi Geol. Camerti 2003, 39–48. Available online: http://193.204.8.201:8080/jspui/handle/1336/862 (accessed on 1 January 2023).
  47. Di Celma, C.; Cantalamessa, G.; Didaskalou, P.; Lori, P. Sedimentology, architecture, and sequence stratigraphy of coarse-grained, submarine canyon fills from Pleistocene (Gelasiane-Calabrian) of the Peri-Adriatic basin, central Italy. Mar. Pet. Geol. 2010, 27, 1340–1365. [Google Scholar] [CrossRef]
  48. Ori, G.G.; Roveri, M.; Vannoni, F. Plio-Pleistocene Sedimentation in the Apenninic-Adriatic Foredeep (Central Adriatic Sea, Italy); Allen, P.A., Homewood, P., Eds.; IAS Special Publication: Gent, Belgium, 1986; Volume 8, pp. 183–198. [Google Scholar]
  49. Dattilo, P.; Pasi, R.; Bertozzi, G. Depositional and structural dynamics of the Pliocene peri-adriatic foredeep, NE Italy. J. Pet. Geol. 1999, 22, 19–36. [Google Scholar] [CrossRef]
  50. Tinterri, R.; Lipparini, L. Seismo-stratigraphic study of the Plio-Pleistocene foredeep deposits of the Central Adriatic Sea (Italy): Geometry and characteristics of deep-water channels and sediment waves. Mar. Pet. Geol. 2013, 42, 30–49. [Google Scholar] [CrossRef]
  51. Regione Abbruzzo. “Piano di Tutela delle Acque”. Relazione Generale e Allegati, Carta dei Complessi Idrogeologici. 2010. Available online: http://www.regione.abruzzo.it/pianoTutelaacque/docs/elaboratiPiano/CartografiaPiano/1_4.pdf (accessed on 1 January 2023).
  52. 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]
  53. Lax, P.; Wendroff, B. Systems of conservation laws. Commun. Pure Appl. Math. 1960, 3, 217–237. [Google Scholar] [CrossRef]
  54. Hou, T.Y.; LeFloch, P.G. Why nonconservative schemes converge to wrong solutions: Error analysis. Math. Comp. 1994, 62, 497–530. [Google Scholar] [CrossRef]
  55. 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).
  56. Cactus Developers. Cactus Computational Toolkit. Available online: http://www.cactuscode.org/ (accessed on 1 January 2023).
  57. Goodale, T.; Allen, G.; Lanfermann, G.; Massó, J.; Radke, T.; Seidel, E.; Shalf, J. The Cactus Framework and Toolkit: Design and Applications. In High Performance Computing for Computational Science—VECPAR 2002; Springer: Berlin, Germany, 2003. [Google Scholar]
  58. Schnetter, E.; Hawley, S.H.; Hawke, I. Evolutions in 3D numerical relativity using fixed mesh refinement. Class. Quantum Grav. 2004, 21, 1465–1488. [Google Scholar] [CrossRef]
  59. Schnetter, E.; Diener, P.; Dorband, E.N.; Tiglio, M. A multi-block infrastructure for three-dimensional time-dependent numerical relativity. Class. Quantum Grav. 2006, 23, S553. [Google Scholar] [CrossRef]
  60. 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]
  61. 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]
  62. Freeze, R.A.; Cherry, J.A. Groundwater Book; Prentice-Hall Inc.: Englewood Cliffs, NJ, USA, 1979. [Google Scholar]
  63. 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]
  64. Hamutoko, J.T.; Post, V.E.A.; Wanke, H.; Beyer, M.; Houben, G.; Mapani, B. The role of local perched aquifers in regional groundwater recharge in semi-arid environments: Evidence from the Cuvelai-Etosha Basin, Namibia. Hydrogeol. J. 2019, 27, 2399–2413. [Google Scholar] [CrossRef]
  65. Hernàndez-Diaz, R.; Petrella, E.; Bucci, A.; Naclerio, G.; Feo, A.; Sferra, G.; Chelli, A.; Zanini, A.; Gonzalez-Hernandez, P.; Celico, F. Integrating Hydrogeological and Microbiological Data and Modelling to Characterize the Hydraulic Features and Behaviour of Coastal Carbonate Aquifers: A Case in Western Cuba. Water 2019, 11, 1989. [Google Scholar] [CrossRef]
  66. Celico, F.; Petrella, E.; Celico, P. Hydrogeological behaviour of some fault zones in a carbonate aquifer of Southern Italy: An experimentally based model. Terra Nova 2006, 18, 308–313. [Google Scholar] [CrossRef]
  67. Andersson, J.E.; Ekman, L.; Nordqvist, R.; Winberg, A. Hydraulic testing and modeling of a low-angle fracture zone at Finnsjon, Sweden. J. Hydrol. 1991, 126, 45–77. [Google Scholar] [CrossRef]
  68. Antonellini, M.; Aydin, A. Effect of faulting on fluid flow in porous sandstones: Petrophysical properties. Am. Assoc. Petrol. Geol. Bull. 1994, 78, 355–377. [Google Scholar]
  69. Bense, V.F.; Van den Berg, E.H.; Van Balen, R.T. Deformation mechanisms and hydraulic properties of fault zones in unconsolidated sediments; the Roer Valley Rift System, The Netherlands. Hydrogeol. J. 2003, 11, 319–332. [Google Scholar] [CrossRef]
  70. Bense, V.F.; Gleeson, T.; Loveless, S.E.; Bour, O.; Scibek, J. Fault zone hydrogeology. Earth Sci. Rev. 2013, 127, 171–192. [Google Scholar] [CrossRef]
  71. Caine, J.S.; Evans, J.P.; Forster, C.B. Fault zone architecture and permeability structure. Geology 1996, 24, 1025–1028. [Google Scholar] [CrossRef]
  72. Chester, F.M.; Logan, J.M. Composite planar fabric of gouge from the Punchbowl fault, California. J. Struct. Geol. 1986, 9, 621–634. [Google Scholar] [CrossRef]
  73. Balsamo, F.; Storti, F. Grain size and permeability evolution of soft-sediment extensional sub-seismic and seismic fault zones in high-porosity sediments from the Crotone basin, southern Apennines, Italy. Mar. Pet. Geol. 2010, 27, 822–837. [Google Scholar] [CrossRef]
  74. Pizzati, M.; Balsamo, F.; Storti, F. Displacement-dependent microstructural and petrophysical properties of deformation bands and gouges in poorly lithified sandstone deformed at shallow burial depth (Crotone Basin, Italy). J. Struct. Geol. 2020, 137, 104069. [Google Scholar] [CrossRef]
  75. Acher, A.; Boiderie, P.; Yaron, B. Soil pollution by petroleum products. I. Multiphase migration of kerosene components in soil columns. J. Contam. Hydrol. 1989, 4, 333–345. [Google Scholar] [CrossRef]
Figure 1. Study area (from [51] Abruzzo, Regione, et al. “Piano di Tutela delle Acque.” Relazione Generale e Allegati, Carta dei Complessi Idrogeologici (2010), modified. Red dots on the map indicate debris sediment).
Figure 1. Study area (from [51] Abruzzo, Regione, et al. “Piano di Tutela delle Acque.” Relazione Generale e Allegati, Carta dei Complessi Idrogeologici (2010), modified. Red dots on the map indicate debris sediment).
Water 15 01900 g001
Figure 2. Hydrogeological section.
Figure 2. Hydrogeological section.
Water 15 01900 g002
Figure 3. Example of the three-dimensional grid geometry used in the numerical simulation of a three-phase fluid flow (water, LNAPL and air) with a spatial grid resolution of 0.50   m and a grid dimension of 160 × 80 × 18   m , at the initial time t   = 0   s . The immiscible contaminant is situated at the top of the parallelepiped on the z x plane (left-hand side) and the z y plane (right-hand side), respectively. The green rectangle zone in each panel is amplified in the upper area at each time indicated.
Figure 3. Example of the three-dimensional grid geometry used in the numerical simulation of a three-phase fluid flow (water, LNAPL and air) with a spatial grid resolution of 0.50   m and a grid dimension of 160 × 80 × 18   m , at the initial time t   = 0   s . The immiscible contaminant is situated at the top of the parallelepiped on the z x plane (left-hand side) and the z y plane (right-hand side), respectively. The green rectangle zone in each panel is amplified in the upper area at each time indicated.
Water 15 01900 g003
Figure 4. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) in a dry soil using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 2,080,768 s. A hydraulic gradient of 0.04. Left-hand side shows the saturation contours in the z x plane. Right-hand side shows the saturation contours on the y x one. The spill is released from an oil pipeline at x , y , z = 0 , 0 , 1   m .
Figure 4. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) in a dry soil using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 2,080,768 s. A hydraulic gradient of 0.04. Left-hand side shows the saturation contours in the z x plane. Right-hand side shows the saturation contours on the y x one. The spill is released from an oil pipeline at x , y , z = 0 , 0 , 1   m .
Water 15 01900 g004
Figure 5. Three-dimensional numerical simulation results of the depth as a function of the water saturation S w (blue points), gasoline saturation S n (red points), and air saturation S a (green points) at various times for a gasoline leak in Figure 4. Initially, at t = 0 s, there is a sharp front of contaminant saturation situated on top of the grid, rapidly going to zero when height decreases. At the same time, it is filled by the air saturation (green point) in the unsaturated zone and the water saturation in the saturated zone. At later times, the contaminant will have already reached the aquifer zone. However, a small part enters the saturated zone (center) and remains floating while moving toward the groundwater flow.
Figure 5. Three-dimensional numerical simulation results of the depth as a function of the water saturation S w (blue points), gasoline saturation S n (red points), and air saturation S a (green points) at various times for a gasoline leak in Figure 4. Initially, at t = 0 s, there is a sharp front of contaminant saturation situated on top of the grid, rapidly going to zero when height decreases. At the same time, it is filled by the air saturation (green point) in the unsaturated zone and the water saturation in the saturated zone. At later times, the contaminant will have already reached the aquifer zone. However, a small part enters the saturated zone (center) and remains floating while moving toward the groundwater flow.
Water 15 01900 g005
Figure 6. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) in a dry soil using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 3686 s; (b) t = 7168 s; (c) t = 102,400 s; (d) t = 1,040,998 s. A hydraulic gradient of 0.04. The spill is released from an oil pipeline at x , y , z = 0 , 0 , 10   m . Notice that the first panel corresponds to a time equal to 3686 s rather than the initial time equal to zero.
Figure 6. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) in a dry soil using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 3686 s; (b) t = 7168 s; (c) t = 102,400 s; (d) t = 1,040,998 s. A hydraulic gradient of 0.04. The spill is released from an oil pipeline at x , y , z = 0 , 0 , 10   m . Notice that the first panel corresponds to a time equal to 3686 s rather than the initial time equal to zero.
Water 15 01900 g006
Figure 7. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) in a dry soil, using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times on the y   x plane. A hydraulic gradient of 0.004. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m .
Figure 7. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) in a dry soil, using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times on the y   x plane. A hydraulic gradient of 0.004. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m .
Water 15 01900 g007
Figure 8. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 482,508 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.20 .
Figure 8. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 482,508 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.20 .
Water 15 01900 g008
Figure 9. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 409,600 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.50 .
Figure 9. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 409,600 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.50 .
Water 15 01900 g009
Figure 10. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 1228 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 285,900 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.20 .
Figure 10. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 1228 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 285,900 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.20 .
Water 15 01900 g010
Figure 11. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 301,875 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.50 .
Figure 11. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 301,875 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.50 .
Water 15 01900 g011
Figure 12. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) in a dry soil, using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 780,492 s. A hydraulic gradient of 0.04. The spill is released from an oil pipeline at x , y , z = 0 , 0 , 1   m .
Figure 12. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) in a dry soil, using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 0 s; (b) t = 204 s; (c) t = 102,400 s; (d) t = 780,492 s. A hydraulic gradient of 0.04. The spill is released from an oil pipeline at x , y , z = 0 , 0 , 1   m .
Water 15 01900 g012
Figure 13. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) in a dry soil, using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 3686 s; (b) t = 102,400 s; (c) t = 307,200 s; (d) t = 794,419 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m .
Figure 13. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) in a dry soil, using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 3686 s; (b) t = 102,400 s; (c) t = 307,200 s; (d) t = 794,419 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m .
Water 15 01900 g013
Figure 14. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 482,304 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.20 .
Figure 14. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 482,304 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.20 .
Water 15 01900 g014
Figure 15. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 409,600 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.50 .
Figure 15. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 18 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 409,600 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 1   m . The unsaturated zone has a S w = 0.50 .
Water 15 01900 g015
Figure 16. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times (a) t = 3686 s; (b) t = 20,480 s; (c) t = 102,400 s; (d) t = 285,900 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.20 .
Figure 16. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times (a) t = 3686 s; (b) t = 20,480 s; (c) t = 102,400 s; (d) t = 285,900 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.20 .
Water 15 01900 g016
Figure 17. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 300,851 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.50 .
Figure 17. Three-dimensional numerical results on the saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, diesel oil, and air) using a spatial grid resolution of 0.50 m and a grid dimension of 160 × 80 × 27 m, at different times: (a) t = 204 s; (b) t = 3686 s; (c) t = 102,400 s; (d) t = 300,851 s. A hydraulic gradient of 0.04. The spill was released from an oil pipeline at x , y , z = 0 , 0 , 10   m . The unsaturated zone has a S w = 0.50 .
Water 15 01900 g017
Figure 18. Depth versus arrival time at the groundwater table for dry soil, for two different LNAPLs. The data was taken from the numerical simulations of the saturation contours previously shown in Table 3.
Figure 18. Depth versus arrival time at the groundwater table for dry soil, for two different LNAPLs. The data was taken from the numerical simulations of the saturation contours previously shown in Table 3.
Water 15 01900 g018
Figure 19. Water saturation in the unsaturated zone versus arrival time at the groundwater table, for two different species of LNAPL. Red line corresponds to gasoline and green line corresponds to diesel oil.
Figure 19. Water saturation in the unsaturated zone versus arrival time at the groundwater table, for two different species of LNAPL. Red line corresponds to gasoline and green line corresponds to diesel oil.
Water 15 01900 g019
Figure 20. Arrival time at the groundwater table vs. water saturation of the unsaturated zone for two different pressure spills of the gasoline.
Figure 20. Arrival time at the groundwater table vs. water saturation of the unsaturated zone for two different pressure spills of the gasoline.
Water 15 01900 g020
Figure 21. Model refinement. Saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a grid resolution of 0.25 m ((a) t = 204 s; (c) t = 614 s), and a grid resolution 0.50 m ((b) t = 204 s; (d) t = 604 s, see also Figure 4) in the case of a gasoline spill from the oil pipeline.
Figure 21. Model refinement. Saturation contours σ n = S n ϕ of a three-phase immiscible fluid flow (water, gasoline, and air) using a grid resolution of 0.25 m ((a) t = 204 s; (c) t = 614 s), and a grid resolution 0.50 m ((b) t = 204 s; (d) t = 604 s, see also Figure 4) in the case of a gasoline spill from the oil pipeline.
Water 15 01900 g021
Table 1. Definitions of the parameters used in the numerical simulations of a gasoline spill from an oil pipeline.
Table 1. Definitions of the parameters used in the numerical simulations of a gasoline spill from an oil pipeline.
ParameterSymbolValue
Absolute permeability k 2.059 × 10 11   m 2
Rock compressibility c R 4.35 × 10 7   Pa 1
Porosity ϕ 0 0.43
Water viscosity μ w 10 3   kg / ms
Water density ρ w 10 3 kg / m 3
Oil (gasoline) viscosity μ n 4.5 × 10 4   kg / ms
Oil (gasoline) density ρ n 750   kg / m 3
Air viscosity μ a 1.8 × 10 5   kg / ms
Air density ρ a 1.225   kg / m 3
Van Genuchten n , m 2.68 , 1 1 2.68
Irreducible wetting phase saturation S w i r 0.045
Superficial tension air-water σ a w 6.5 × 10 2   N / m
Interfacial tension in nonaqueous water σ n w 2.6 × 10 2   N / m
Capillary pressure of air-water at zero saturation p c a w 0 676.55   Pa
Capillary pressure air-nonaqueous at zero saturation p c a n 0 405.93   Pa
Table 2. Definitions of the parameters used in the numerical simulations of a diesel oil spill from an oil pipeline.
Table 2. Definitions of the parameters used in the numerical simulations of a diesel oil spill from an oil pipeline.
ParameterSymbolValue
Absolute permeability k 2.059 × 10 11   m 2
Rock compressibility c R 4.35 × 10 7   Pa 1
Porosity ϕ 0 0.43
Water viscosity μ w 10 3   kg / ms
Water density ρ w 10 3 kg / m 3
Oil (diesel oil) viscosity μ n 3.61 × 10 3   kg / ms
Oil (diesel oil) density ρ n 830   kg / m 3
Air viscosity μ a 1.8 × 10 5   kg / ms
Air density ρ a 1.225   kg / m 3
Van Genuchten n , m 2.68 , 1 1 2.68
Irreducible wetting phase saturation S w i r 0.045
Superficial tension air-water σ a w 6.5 × 10 2   N / m
Interfacial tension in nonaqueous water σ n w 3.0 × 10 2   N / m
Capillary pressure of air-water at zero saturation p c a w 0 676.55   Pa
Capillary pressure air-nonaqueous at zero saturation p c a n 0 374.68   Pa
Table 3. Analysis of arrival time at the groundwater table and the position in the x coordinate under different parameter conditions: two different densities, two different hydraulic gradients, and five different depths from the oil pipeline spill and dry soil. The “-” sign indicates no available results, i.e., the contaminant is still moving in the unsaturated zone and has not reached the groundwater table.
Table 3. Analysis of arrival time at the groundwater table and the position in the x coordinate under different parameter conditions: two different densities, two different hydraulic gradients, and five different depths from the oil pipeline spill and dry soil. The “-” sign indicates no available results, i.e., the contaminant is still moving in the unsaturated zone and has not reached the groundwater table.
Type (Density)Thickness of the Unsaturated Zone (m)Hydraulic GradientArrival Time at the Groundwater Table (s)Position in x after One Day and 4.4 h (m)
Gasoline1.00.04 204.8 −16.0
0.004 204.8 −14.0
Diesel oil1.00.04 204.8 −6.0
0.004 204.8 −5.5
Gasoline 2.00.04 204.8 −16.5
0.004 204.8 −15.0
Diesel oil2.00.04614.4−6.0
0.004614.4−5.0
Gasoline5.00.041228.8−16.5
0.0041228.8−15.0
Diesel oil5.00.0410,854.4−3.5
0.00412,288.0−3.5
Gasoline10.00.046348.8−14.0
0.0046348.8−12.5
Diesel oil10.00.04286,720.0-
0.004276,480.0-
Gasoline20.00.0492,160.0−2.5
0.00492,160.0−2.5
Diesel oil20.00.04--
0.004--
Table 4. Analysis of arrival time at the groundwater table and the position in the x coordinate under different parameter conditions: two different types of densities, two different depths from the oil pipeline spill and, water saturation of the unsaturated zone equal to 0.20 and 0.50. The hydraulic gradient is 0.04.
Table 4. Analysis of arrival time at the groundwater table and the position in the x coordinate under different parameter conditions: two different types of densities, two different depths from the oil pipeline spill and, water saturation of the unsaturated zone equal to 0.20 and 0.50. The hydraulic gradient is 0.04.
Unsaturated Zone Depth (m)Type of
Contaminant
Water Saturation in the Unsaturated ZoneArrival Time of to the Groundwater Table (s)Position in x after 1 Day and 4.4 h (s)
1.0Gasoline0.0 204.8 −16.0
0.2 204.8 −16.0
0.5 204.8 −16.0
1.0Diesel oil0.0 204.8 −6.0
0.2 204.8 −6.0
0.5 204.8 −6.0
10.0Gasoline0.06348.8−14.0
0.24915.0−16.0
0.54710.4−16.0
10.0Diesel oil 0.0286,720.0-
0.2151,552.0-
0.573,720.0−1.5
Table 5. Water saturation of the unsaturated zone vs. arrival time at the groundwater table for a gasoline spill at the atmospheric pressure.
Table 5. Water saturation of the unsaturated zone vs. arrival time at the groundwater table for a gasoline spill at the atmospheric pressure.
Unsaturated Zone Depth (m)Type of ContaminantWater Saturation in the Unsaturated Zone Arrival Time at the Groundwater Table (s)
1.0Gasoline0.01024.0
1.0Gasoline0.20819.2
1.0Gasoline0.50614.4
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.; Pinardi, R.; Scanferla, E.; Celico, F. 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. Water 2023, 15, 1900. https://doi.org/10.3390/w15101900

AMA Style

Feo A, Pinardi R, Scanferla E, Celico F. 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. Water. 2023; 15(10):1900. https://doi.org/10.3390/w15101900

Chicago/Turabian Style

Feo, Alessandra, Riccardo Pinardi, Emanuele Scanferla, and Fulvio Celico. 2023. "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" Water 15, no. 10: 1900. https://doi.org/10.3390/w15101900

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