Developing a Sustainable Petroleum Transportation Framework to Reduce Environmental Impact During Emergency Scenarios

1. Introduction

Petroleum, also known as crude oil, remains one of the most essential sources of energy. However, because of their non-sustainable and non-renewable nature, petroleum resources eventually deplete, while their extraction or transportation can cause environmental pollution and human health concerns [1,2,3,4,5,6]. Moreover, environmental conditions are crucial to the fate and impact of a crude oil spill [7,8]. Therefore, it is necessary to study ways to safely transport it into the environment [9,10].

Crude oil is classified by geographic location, sulfur content, and API gravity [11], a measure of oil density that indicates whether the crude oil is heavy or light compared to water. If API gravity is greater than 10, then the petroleum floats; if it is less than 10, it sinks. Although most crude oils are light nonaqueous phase liquids (LNAPLs), with densities lower than water, they may contain asphaltenes, which are heavy organic macromolecules classified as dense nonaqueous phase liquids (DNAPLs) [12]. It can also be considered its viscosity, which typically ranges approximately from $0.01kg/\left(ms\right)$ to $100kg/\left(ms\right)$, and relies on the values of the American Petroleum Institute (API) gravity and temperature [13],[14]. Usually, to transport crude oil with elevated viscosity through a pipeline, it is necessary to increase the temperature using a hot-pumping method. When pumping is shut down, the system temperature is decreased, and therefore the crude oil viscosity increases [14,15].

Oil seepage and spills are environmental disasters with lasting, harmful effects on marine ecosystems. Recently, the introduction of nanoparticles into contaminated areas to enhance bioremediation of oil spills has emerged as a successful method, especially for marine microorganisms [16,17]. Other oil spills cleanup techniques include mechanical methods using containment booms, and an underflow dam [18], in situ chemical oxidation of petroleum hydrocarbons [19], thermal methods for a rapid recovery of oil and microplastics in marine oil spill co-contamination [20], bioremediation techniques [21] using green composites [22], or natural mineral sorbents as green materials [23]. All these techniques have advantages and disadvantages, as the efficacy of each depends on factors such as viscosity, volume, density, location, and weather conditions [16].

Numerical modeling of crude oil leakage and migration from pipelines and storage facilities is challenging. Risk analysis using different models and numerical methods has been studied for petroleum/crude oil migration due to stress-induced corrosion, corroded ions, and microorganisms [24]. The length of parallel pipelines for transporting petroleum is continuously increasing. Ref. [25] uses numerical simulations to investigate the impact of oil spill and fires accidents nearby natural gas pipelines and indicates that vapor concentration and the size of are the most critical factors. Ref. [26], analyzes the migration from an accidental hole in a tank containing petroleum and the influence of hole diameter, tank height, and tank width. Ref. [27] uses the lattice Boltzmann method to investigate hydrocarbon migration and accumulation and multi-scale quantitative numerical simulation method and a comparison with a sandbox experiment results. Other numerical results include FEM (finite element method) [28] with COMSOL, FDM (finite difference method) [29] investigates BTEX plume using MT3DMS, lattice Boltzmann method [30,31], smooth particle hydrodynamics (SPH) [32].

This paper investigates the influence of hydrogeological features, such as the porosity of the porous medium, and different climatic conditions, including rainfall-induced water saturation, on petroleum migration in the variable-saturated zone. The objective is to select a better location based on its hydrogeological features for the possible construction of onshore petroleum pipelines, as well as to suggest possible technical solutions in designing the pipeline features and give useful information to optimize remediation activities in emergency scenarios. Saturation is linked to hydraulic conductivity and provides information on the type of porous soil that facilitates contaminant migration. Also, porosity can influence the velocity of the vertical migration in the unsaturated zone. We analyze the risk analysis of various scenarios that mimic a pipeline leakage (where the pipeline is included in the numerical simulation) and the consequent spread of the contaminant material in the environment.

We use the numerical CactusHydro code [33,34], that has been validated for several scenarios, in particular, using a sandbox laboratory experiment [35], and various applications including PCE extraction in potential emergency scenarios [36]. The CactusHydro code accurately resolves the sharp discontinuities and temporal dynamics of three-phase fluid flow, thereby accurately tracking contaminant migration in space and time through a variably saturated zone. We present sixteen numerical simulation results that mimic petroleum outflow from a broken pipeline (the most extreme event) and migration toward the saturated zone under the forces of gravity and capillary pressures, to quantify the environmental damage that can be caused over time. Due to its elevated viscosity, crude oil has very limited mobility; for this reason, it migrates very little vertically, and most of it moves under the pressure from the pipeline outflow. The broken pipeline is included in the numerical simulations and modeled as a low-permeability material. Previous investigations using pipelines, but not with petroleum, and with different geometries have been reported in Refs. [37,38].

2. Materials and Methods

2.1. Hydrogeological and Geological Data

In Italy, natural hydrocarbon emissions are mostly located in the Southern Apennine thrust belt [39]. Based on the structural units present in the area and their relative permeabilities, distinct hydrogeological complexes are identified. For example, carbonate complex and high permeability due to fracturing and karst phenomena, and alluvional complex with variable permeability, among others.

We investigate the possible consequences of an emergency scenario in the area where a possible network of onshore pipelines could be constructed to transport crude oil, with hydrogeological features including an alluvial complex and a lithology of alternating clay, gravel, and sand. Currently, comparison with experimental results is not possible because such data are not publicly available. However, this study employs real hydrogeological values consistent with previously described lithological features. Therefore, the results are applicable and reproducible in any region exhibiting similar hydrogeological characteristics.

In this area, a risk analysis is developed to estimate the quantity of petroleum released into the environment in the event of an accidental loss of containment from an onshore pipeline. We are also interested in the time required for crude oil to migrate from the leak site toward the groundwater table, and the quantity of petroleum mass released, under different climate conditions, to plan the time needed to begin possible remediation. Since crude oil has very limited mobility due to its high dynamic viscosity, we are interested in determining the minimum time beyond which crude oil does not reach the groundwater table under various climate conditions, including the dry season and rainfall-induced saturation, for different values of the porous medium’s porosity. We are interested in providing guidelines to minimize the damage of crude oil spills in emergency situations.

2.2. Mathematical Model

For a three-phase immiscible fluid flow composed of nonaqueous (n), i.e., petroleum, water (w), and air (a) in a variably saturated zone, CactusHydro resolves the system of three partial differential equations as a function of the saturations of the three-phase fluid flow, and pressure variables: $S_n$, $S_w$, $S_a$, and $p_a$. The conservation relation, $S_n+S_w+S_a=1$ holds. See Refs. [33,34,37]. The rock compressibility $c_R$ as a function of the porosity and pressure is given by: $c_R={\displaystyle \frac{1}{\varphi }}{\displaystyle \frac{\partial \varphi }{\partial p}}$. We use a linear approximation of this expression, obtained by a Taylor expansion up to order one in $c_R$. The resulting porosity is given by: $\varphi ={\varphi }_0[1+c_R(p-p_0)]$, where ${\varphi }_0$ is the porosity value at the reference pressure $p_0$ (the atmospheric). We use the product of the saturation times the porosity: ${\sigma }_n=\varphi S_n$, ${\sigma }_w=\varphi S_w$, and ${\sigma }_a=\varphi S_a$. The system of equations can be written as (assuming constant density and viscosity for each phase fluid flow),

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

(1)

and

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

(2)

where $\alpha =(n,w,a)$, and

$$F_{\left(\alpha \right)}^i(S_n,S_w,S_a,p)=-{\displaystyle \frac{k_{r\left(\alpha \right)}(S_n,S_w,S_a)}{{\mu }_{\left(\alpha \right)}}}{k}^{ij}\left({\displaystyle \frac{\partial p}{\partial x^i}}+{\rho }_{n}g{\displaystyle \frac{\partial z}{\partial x^i}}\right)$$

(3)

do not depend on the spatial derivative of the saturation, while

$$Q_{\left(\alpha \right)}^i(S_n,S_w,S_a,p)=-{\displaystyle \frac{k_{r\left(\alpha \right)}(S_n,S_w,S_a)}{{\mu }_{\left(\alpha \right)}}}{k}^{ij}{\displaystyle \frac{\partial p_{ca\left(\alpha \right)}(S_n,S_w,S_a)}{\partial x^i}}$$

(4)

depend on the spatial derivative of the saturation. In the regime of interest, the dominant part is the hyperbolic one (Eq.(3)), and is the one that produces and propagating the shock fronts. Eq.(4) instead, is the parabolic part that depends on the capillary pressures, and it is also responsible for the lateral migration of nonacqueous contaminant. CactusHydro uses a High-Resolution Shock-Capturing (HRSC) flux conservative method [40,41,42] that treats the hyperbolic (advective) part of the flow equations better than other numerical methods. See Refs. [33,34] for details. This technique requires a time step size sufficiently small and for this reason the use of High-Performance Computing (HPC). The numerical code is based on an open source software framework to develop parallel HPC simulation codes called Cactus computational toolkit [43,44,45]. Data are evolved on a cartesian mesh wit six refinement level using Carpet [46,47].

The relative permeabilities, $k_{rn}$, $k_{rw}$, and $k_{ra}$, are given by Parker et al. [48]:

$$\begin{array}{cc}\hfill k_{rw}& =S_{ew}^{1/2}[1-{(1-S_{ew}^{1/m})}^m]2\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill k_{ra}& ={(1-S_{et})}^{1/2}{(1-S_{et}^{1/m})}^{2m}\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill k_{rn}& ={(S_{et}-S_{ew})}^{1/2}{[{(1-S_{ew}^{1/m})}^m-{(1-S_{et}^{1/m})}^m]}^2,\hfill \end{array}$$

(7)

where $m=1-1/n$, and $S_{et}={\displaystyle \frac{S_w+S_n-S_{wir}}{1-S_{wir}}}$ is the total effective liquid saturation, and $S_{wir}$ is the irreducible wetting phase saturation. For the capillary pressure we use the van Genuchten model [49], which for the effective saturation has the following form: $S_e={[1+{\left(\alpha h_c\right)}^n]}^{(1-1/n)}$, where $h_c={\displaystyle \frac{p_c}{{\rho }_{w}g}}$ is the capillary pressure head and $\alpha$, m, and n are model parameters.

2.3. Hydrogeological Features and Crude Oil Parameters

The 3D numerical model investigated in this paper simulates an emergency scenario in which an onshore pipeline transporting high-pressure petroleum ruptures, releasing petroleum into the environment. The system is then turned off after one hour, and the spatial and temporal distribution of the petroleum is examined. Since crude oil has a high viscosity, its mobility is very limited. Therefore, we are interested in defining a temporal window during which the petroleum remains in the unsaturated zone and can be mechanically remediated before reaching the groundwater table. We applied the numerical code CactusHydro, originally introduced in Refs. [33,34], which has been checked and expanded over the years using several different initial and boundary conditions and simulates different scenarios of the petroleum plume migration in a variably saturated zone. We briefly summarized the mathematical model we use, which will help fix notation and clarify parameter notation.

The hydraulic conductivity is considered homogeneous and isotropic in the porous medium, and consistent with the hydrogeological features previously mentioned, and it is given by $K=5.0\times {10}^{-4}m/s$. Consequently, the permeability k, where $k={\displaystyle \frac{K{\mu }_w}{{\rho }_{w}g}}$, is given by $k=5.102\times {10}^{-11}{m}^2$, where $g=9.8m/s^2$ is the gravity constant on Earth, ${\rho }_w={10}^{3}kg/m^3$ is the density of the water, and ${\mu }_w={10}^{-3}kg/\left(ms\right)$ its dynamic viscosity. The porosity is fixed to ${\varphi }_0=0.43$, at the reference pressure $p_0$ (the atmospheric, $p_0=1.013\times {10}^{5}Pa$), and it is compatible with ’sand’ [50]. The rock compressibility, $c_R$ for ’sand’ is given in Ref. [51]. The hydrogeological parameter features are listed in Table 1.

The petroleum crude oil investigated in this work has a density of ${\rho }_n=831.8kg/m^3$ and an elevated dynamic viscosity, ${\mu }_n=3.107kg/\left(ms\right)\equiv 3.1\times {10}^{3}cp$ [52] (compared with previous investigations, as for example, ${\mu }_n=3.61\times {10}^{-3}kg/\left(ms\right)$ for the diesel oil [37], and ${\mu }_n=4.5\times {10}^{-4}kg/\left(ms\right)$ for the gasoline [38]). The dynamic viscosity limits its displacement, which is why a pressure of $2.0\times {10}^{6}Pa$ is required along the pipeline. The numerical simulation mimics a scenario in which the onshore pipeline is completely broken, and crude oil is released at a pressure of $2.0\times {10}^6$ Pascal for 3600 seconds. This is how long it takes to intervene and shut down the pipeline.

The numerical model also includes the broken pipeline (not just the volume cell where the contaminant is ejected), represented as an impermeable parallelepiped with a very low permeability, $k_p=5.102\times {10}^{-15}{m}^2$. We also assumed a hydraulic gradient of 0.01, so that crude oil, as an LNAPL, would remain near the groundwater table and eventually be transported with the flow. In this work, we aim to investigate the time at which remediation can be initiated before the contaminant reaches the groundwater table. We are not interested in letting the crude oil reach the groundwater table, also because the time is very long. In any case, the bottom part of the parallelepiped also has a low permeability, $k_b=5.102\times {10}^{-15}{m}^2$, similar to $k_p$.

The capillary pressures in Table 1 are calculated from the superficial tensions. The superficial tension air-nonaqueous is ${\sigma }_{an}=0.034N/m$ [53]. The interfacial tension air-water is ${\sigma }_{aw}=0.065N/m$ [54]. From here we get, ${\beta }_{an}={\displaystyle \frac{{\sigma }_{aw}}{{\sigma }_{an}}}={\displaystyle \frac{0.065}{0.034}}=1.912$, and $p_{can}={\displaystyle \frac{p_{caw}}{{\beta }_{an}}}=353.84Pa$, where $p_{caw}={\rho }_{w}g/\alpha$, and $\alpha =14.5m^{-1}$ (from here, $n=2.68$ [50]). Using the relation that for three-phase fluid flow only two capillary pressures are independent, we have $p_{cnw}=p_{caw}-p_{can}=322.71Pa$.

2.4. The Numerical Model, Initial and Boundary Conditions of the Petroleum Leakage

We consider a petroleum leakage with the features given in Table 1 from an onshore pipeline located in a variably saturated zone at one meter depth. We analyze two different initial positions for the petroleum spill: a) $(x,y,z)=(0,0,2)m$, and b) $(x,y,z)=(0,0,5)m$, from the groundwater table surface, which crosses the point $(x,y,z)=(0,0,0)m$. For case a), the grid geometry shown is a parallelepiped with an extension of $60m$ long from $x=[-40,+20]m$, $30m$ wide from $y=[-15,15]m$, and $19m$ depth from $z=[+3,-16]m$. For case b), the grid geometry shown is a parallelepiped with similar measures in x and y, and a depth of $22m$ from $z=[+6,-16]m$. The petroleum spill is released with a pressure of $2.0265\times {10}^{6}Pa$ for one hour, and has an initial saturation of one, $S_n(0,0,2,t=0)=1;S_n(0,0,5,t=0)=1$. The remaining vadose zone is at atmospheric pressure. A low permeability zone ($k_b$) is situated between $z=[-15,-16]m$ in both cases. We also constructed a real onshore pipeline (depicted in the Figures with light grey) using a low permeability value ($k_p$) that goes from $DDx=[-40,-0.25]m$, $DDz=[2.0,2.5]m$, and $LLx=[0.25,20]m$, $LLz=[2.0,2.5]m$, for case a). For case b), the onshore pipeline goes from $DDx=[-40,-0.25]m$, $DDz=[5.0,5.5]m$, and $LLx=[0.25,20]m$, $LLz=[5.0,5.5]m$. Boundary conditions are no-flow, except at the top of the parallelepiped in the infiltration zone.

In this work we performed the following transient numerical simulations of crude oil,

• petroleum spill at $(x,y,z)=(0,0,2)$ from the groundwater table, and $S_w=(0.0,0.10,0.20,0.30,0.40)$ in the unsaturated zone,
• petroleum spill at $(x,y,z)=(0,0,5)$ from the groundwater table, and $S_w=(0.0,0.10,0.20,0.30,0.40)$ in the unsaturated zone,
• petroleum spill at $(x,y,z)=(0,0,5)$ from the groundwater table, $S_w=0.0$, and initial porosity values ${\varphi }_0=(0.10,0.25,0.60,0.75,0.90)$ in the variably saturated zone.

where these scenarios represent conditions such as dry season (corresponding to zero water saturation), rainfall-induced saturation and different porosities texture.

3. Results and Discussions

In this section, we present the results of sixteen transient numerical simulations of petroleum crude oil migration in a variably saturated zone following a spill from an onshore pipeline. The spill migration is investigated at two different distances from the groundwater table, $(0,0,2)m$ and $(0,0,5)m$. Moreover, the migration of the leaked petroleum is investigated using different values of water saturation in the unsaturated zone to mimic different climatic conditions, $S_w=[0.0,0.10,0.20,0.30,0.40]$, corresponding to different situations starting with dry soil and increasing constant rainfall percolation. The migration of the petroleum spill is also investigated using different values of the porosity in the porous medium, ${\varphi }_0=[0.10,0.25,0.43,0.60,0.75,0.90]$, with constant $S_w=0.0$, and a petroleum spill at $(0,0,5)m$.

3.1. 3D Transient Numerical Simulations Results for Water Saturation Values $S_w=[0.0,0.1,0.2,0.3,0.4]$ with a Spill Situated at $(0,0,2.0)m$, and ${\varphi }_0=0.43$

Figure 1 shows the migration of a spilled crude oil at $(0,0,2.0)m$, in a variably saturated zone from an onshore pipeline. The unsaturated zone is dry soil, $S_w=0.0$, and the transient numerical simulations are at different times. The left-hand side shows the 3D numerical saturation contours, $({\sigma }_n=S_n\varphi ,{\sigma }_w=S_w\varphi ,{\sigma }_a=S_a\varphi ,)$, of the three-phase immiscible fluid flow (crude oil, water, air) versus the position, at different times, and corresponds to the $z-x$ plane. The right-hand side instead shows similar results, but in the $z-y$ plane. We show a parallelepiped with dimensions $60m\times 30m\times 19m$, with a grid spatial resolution of $\Delta x=\Delta y=\Delta z=0.50$. Zoomed inset on the $(z-x)$ plane. Color bars on the right state ${\sigma }_n=S_n\varphi$, and ${\sigma }_w=S_w\varphi$. On the left-hand side, the onshore pipeline is depicted in light grey and simulated with a low permeability value (see Table 1).

After the release of petroleum from the onshore pipeline, it moves downward toward the saturated zone due to the gravity force, see Figure 1 a), which shows the migration at $t=204s$. Since petroleum has very limited mobility after 31 days and 8.1 hours, it remains in the unsaturated zone (and has not reached the saturated zone). See Figure 1 b).

A similar behavior happens if the water saturation of the unsaturated zone increases. See for example, Figure 2 with $S_w=0.10$, Figure 3 with $S_w=0.20$, Figure 4 with $S_w=0.30$, Figure 5 with $S_w=0.40$, respectively. They have similar behavior after 204 seconds, see Figure 2 a), Figure 3 a),Figure 4 a), Figure 5 a), and basically the only difference is the water saturation of the unsaturated zone.

In Figure 2 b), the petroleum did not reach the groundwater table after 20 days and 2.2 hours. In Figure 3 b), the petroleum just reaches the groundwater table after 20 days and 6.4 hours. In Figure 4 b), the contaminant arrived at the groundwater table after 18 days and 11.1 hours. In Figure 5 b), the contaminant arrives at the groundwater table after 15 days and 22.7 hours. It seems that the water content (water saturation) in the vadose zone increases the contaminant’s velocity of arrival, even though a smaller quantity of contaminant is released as $S_w$ increases, because water occupies part of the porous medium. The same happens for the right-hand side of each of these figures. Also, the contaminant concentration ${\sigma }_n$, decreases as $S_w$ increases. Notice that the right-hand side of each figure is symmetric, since it represents the $z-y$ plane. We will conduct a quantitative study of this issue to determine the mass released in each case.

It is worth noting that the groundwater surface level rises as water saturation increases, since the boundary conditions are held constant in the head. As a result, we observe that petroleum is moving more rapidly toward the groundwater table as the groundwater table rises. From the previous figures, one notices that as far as the water saturation of the unsaturated zone increases, it takes less time for crude oil to arrive near the groundwater table. To resolve this issue, the following subsection presents a migration for a scenario similar to the previous one, but using a spill located 5 meters from the groundwater table.

Table 2,Table 3,Table 4,Table 5,Table 6 show the elevation of the spilled crude oil at the initial position $(x,y,z)=(0,0,2)m$, water saturation $S_w=(0.0,0.10,0.20,0.30,0.40)$, respectively, and different times. The data are the one used for Figure 1,Figure 2,Figure 3,Figure 4,Figure 5, respectively. We show the distribution of the petroleum mass as a function of elevation for different values of $S_w$ and times: 1 hour, 5 hours, 1 day, 5 days, 15 days, etc, and the percentage of trapped contaminant spread all over the space.

Figure 6 shows the elevation as a function of the crude oil mass and the spill at $(x,y,z)=(0,0,2)m$, at fixed water saturation in the unsaturated zone and at different times. To reconstruct the spilled crude oil mass, we use Roe’s Riemann Solver [55] and the data in the previous Table 2,Table 3,Table 4,Table 5,Table 6. Basically, the amount of contaminant initially located on top of the grid slowly moves downward for all values of water saturation in the unsaturated zone.

After 1 hour, the contaminant is almost at the top of the pipeline, and the situation is the same for all values of $S_w$. As time increases, the contaminant slowly goes downward due to the gravitational force. After 31 days, for $S_w=0.0$, the contaminant has not yet reached the groundwater table, which crosses the point $(x,y,z)=(0,0,0)m$. See Table 2.

For the case $S_w=0.40$ and after 15 days (bottom panel), the petroleum has reached the same elevation in meters as in the case $S_w=0.0$ after 31 days. This implies that the contaminant is moving downward faster as water saturation in the unsaturated zone increases. See also Table 6. Indeed, for $S_w=0.0$, after 31 days the petroleum has reached the region $0.75<z<1.25$ with a percentage of trapped mass contaminant of 8.93, while for $S_w=0.40$ the petroleum reaches the same region with a percentage of trapped mass contaminant of 10.34 just after 15 days.

Figure 7 shows the elevation as a function of the crude oil mass (spill at $(x,y,z)=(0,0,2)m$) at fixed time and different water-saturation values in the unsaturated zone. Basically, there is no substantial difference in the values of $S_w$ after 1 hour. Similarly, after 5 hours.

For the 5-day and 10-day cases, it is observed that petroleum is moving downward similarly at all saturation values, although its mobility is very limited due to its high viscosity, compared with previous studies using diesel oil [37]. After 15 days, the situation is very similar across all $S_w$ values, but for $S_w=0.40$, the contaminant arrives faster than in the other cases. See the line read in the last panel, which shows a greater quantity of contaminant mass at the elevation of $0.75m$.

3.2. 3D Transient Numerical Simulations Results for Water Saturation Values $S_w=[0.0,0.1,0.2,0.3,0.4]$ with a Spill Situated at $(0,0,5.0)m$, and ${\varphi }_0=0.43$

Let us investigate a similar scenario as before, but using a spilled crude oil situated at $(x,y,z)=(0,0,5)m$ instead of $(x,y,z)=(0,0,2)m$. This is to avoid potential interference between the pipeline location and the rising groundwater surface due to rainfall recharge (see Figure 5).

Figure 8 shows the migration of a spilled crude oil at $(0,0,5.0)m$, in a variably saturated zone from an onshore pipeline. The unsaturated zone is dry soil ($S_w=0.0$), and the transient numerical simulations are conducted at different times. The left-hand side shows the 3D numerical saturation contours (${\sigma }_n=S_n\varphi ,{\sigma }_w=S_w\varphi ,{\sigma }_a=S_a\varphi$) of the th three-phase immiscible fluid flow (crude oil, water, air) versus the position, at different times, and corresponds to the $z-x$ plane. The right-hand side, instead, shows similar results but in the $z-y$ plane.

We show a parallelepiped grid with dimensions $60m\times 30m\times 11m$, with a spatial resolution of $\Delta x=\Delta y=\Delta z=0.50m$. Zoomed inset on the $(z-x)$ plane. Color bars on the right state ${\sigma }_n=S_n\varphi$, and ${\sigma }_w=S_w\varphi$. On the left-hand side, the onshore pipeline is depicted in light grey along the grid and is simulated with a very low permeability, $k_p$ (see Table 1).

Figure 8 and Figure 1, which differ only in the initial position of the spilled crude oil, are very similar in the contaminant’s spatial and temporal migration. See Table 2, Table 7 for $S_w=0.0$. This result is expected, since the contaminant released from the onshore pipeline does not know the distance to the groundwater surface. After 18 days, the petroleum has moved very little and is far from the groundwater surface.

Figure 9,Figure 10,Figure 11,Figure 12 show similar scenarios with $S_w=0.10,S_w=0.20,S_w=0.30$, and $S_w=0.40$, respectively. Since the distance from the spilled contaminant to the groundwater table has increased compared to previous cases, we can now analyze in more detail the impact of vadose-zone water saturation on petroleum mobility. For these figures, we show results for contour saturations of the three-phase fluid flow up to six days. Indeed, after 6 days, the scenario is very similar among the three of them, due to the petroleum’s limited mobility. Nevertheless, it shows that the contaminant moves more easily as $S_w$ increases. See both the left- and right-hand sides of these figures. Also, the contaminant concentration, ${\sigma }_n$, decreases as $S_w$ increases. This is because the contaminant can spread more easily when there is water in the unsaturated zone, and fewer contaminant particles occupy the void spaces (which are now occupied by water). Similar situation happens in the case of spill at $(x,y,z)=(0,0,0)$ (see Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5).

Table 7,Table 8,Table 9,Table 10,Table 11, show the elevation of the spilled crude oil at the initial position $(x,y,z)=(0,0,5)m$, water saturation $S_w=(0.0,0.10,0.20,0.30,0.40)$, respectively, and different times. The data are the one used for Figure 8,Figure 9,Figure 10,Figure 11,Figure 12, respectively. We show the distribution of the petroleum mass as a function of elevation for different values of $S_w$ and times: 1 hour, 5 hours, 1 day, 5 days, etc, and the percentage of trapped contaminant spread all over the space.

From the previous tables, we observe that after 5 days and for the range $3.75<z<4.25$, the spilled contaminant mass is 14.8 Kg and 2.67 % of trapped contaminant; (15.5 Kg, 2.82 %) for $S_w=0.10$; (17.0 Kg, 3.10 %) for $S_w=0.20$; (20.4 Kg, 3.74 %) for $S_w=0.30$; and (22.3 Kg, 4.09 %) for $S_w=0.40$. Moreover, the spilled contaminant mass and the percentage of trapped contaminant mass left after 5 days, in the range $5.25<z<5.75$ is 212.0 Kg and 38.2 % for $S_w=0.0$; (207.4 Kg, 37.7 %) for $S_w=0.10$, (199.1 Kg, 36.34 %) for $S_w=0.20$; (196.7 Kg, 35.9 %) for $S_w=0.30$; and (191.7 Kg, 35.1 %) for $S_w=0.40$. These results suggest that as $S_w$ increases, less contaminant is trapped on top of the grid where the petroleum was released. At the same time, more contaminants are spread in the saturated zone. This suggests that the presence of water (due to rainfall) in the unsaturated zone facilitates the migration of the contaminant toward the saturated zone.

Figure 13 shows the elevation as a function of the crude oil mass at the spill, $(x,y,z)=(0,0,5)m$, at fixed values of $S_w$ in the unsaturated zone and at different times (15 days only for $S_w=0.0$). To reconstruct the spilled crude oil mass, we use Roe’s Riemann Solver [55] and the data in the previous Table 2, Table 3. Table 4, Table 5, Table 6. For dry soil ($S_w=0.0$), the contaminant is reduced over time. Initially, after 1 hour, the contaminant is almost at the top of the pipeline, and over time, it slowly moves downward due to gravity. After 5 days, the contaminant has not yet reached the groundwater table, even though $S_w$ is increasing. In all scenarios in Figure 13, the contaminant never reaches the groundwater table before 15 days. For $S_w=0.40$ (bottom panel), there is more contaminant in the lower part near the saturated zone (see the green line).

Figure 14 shows the elevation as a function of the crude oil mass at $(x,y,z)=(0,0,5)m$, at fixed time values, and for different water-saturation levels in the unsaturated zone. Basically, there is no substantial difference in the value of $S_w$ across the four panels, which correspond to different time intervals (1 hour, 5 hours, 1 day, 5 days). The only difference is noticed in the last panel (after 5 days). Here, we can clearly see that the contaminant is released faster for $S_w=0.40$. See the red line. Also, it arrives downward faster than the other water saturation values. This corroborated the previous figures and data.

3.3. 3D Transient Numerical Simulations Results for Different Initial Porosity Values ${\varphi }_0=[0.10,0.25,0.60,0.75,0.90]$ with a Spill Situated at $(0,0,5.0)m$

Up to here, we were using a fixed initial value of the porosity (${\varphi }_0=0.43$ at the atmospheric pressure) given in Table 1, for all the simulations investigated, since we were interested in describing a porous medium composed of ’sand’ with the hydrogeological features described in the section methods and materials.

Let us now investigate the influence on the porosity value for the scenario $S_w=0.0$ and a petroleum spill at $(0,0,5)m$. The porosity values investigated are ${\varphi }_0=(0.10,0.25,0.60,0.75,0.90)$, where we included the value ${\varphi }_0=0.43$ (previously shown) to facilitate comparisons, and ${\varphi }_0$ is the porosity at $p_0$ (the atmospheric pressure).

Figure 15 shows the saturation contours of three immiscible phases (water, crude oil, and air) as a function of position at different times, for an initial porosity of ${\varphi }_0=0.10$. The first panel shows the initial time $t=204s$, and the second panel shows $t=5$ days and 1.7 hours. As can be seen, the quantity of contaminant mass expelled from the onshore pipeline is smaller than the previous cases with ${\varphi }_0=0.43$, since there are fewer void spaces in the porous medium. However, it migrates more rapidly toward the groundwater table as it spills out with a very high pressure from the pipeline, and the available space to move is reduced compared to the previous cases with ${\varphi }_0=0.43$. See Figure 8 for $S_w=0.0$.

Figure 16, Figure 17, Figure 18, and Figure 19 show a similar scenario as before with ${\varphi }_0=0.10$, but with ${\varphi }_0=(0.25,0.60,0.75,0.90)$, respectively. These results suggest that as ${\varphi }_0$ increases, the amount of expelled contaminant increases, but it moves less. Indeed, the contaminant concentration ${\sigma }_n$ increases accordingly, as shown in these figures. This is due to its highly dynamic viscosity.

Table 12 refers to data used in Figure 15 with ${\varphi }_0=0.10$, and Table 13, Table 14, Table 15, Table 16, shows the data used Figure 16 with ${\varphi }_0=0.25$, Figure 17 with ${\varphi }_0=0.60$, Figure 18 with ${\varphi }_0=0.75$, Figure 19 with ${\varphi }_0=0.90$, respectively. We also include the data in Table 7 for porosity ${\varphi }_0=0.43$. As can be seen from these data, when the porous medium porosity is very low, such as ${\varphi }_0=0.10$ in Table 12, the spilled contaminant mass, after one hour, and in the range $5.25<z<5.75$ meters is 207.5 Kg and 44.8 % of trapped contaminant. As porosity increases, the amount of petroleum that can be accommodated outside the pipeline increases. See for example, Table 13 for ${\varphi }_0=0.25$ with 235.9 Kg (45.7 %), Table 7 for ${\varphi }_0=0.43$ with 261.0 Kg (46.2 %), Table 14 for ${\varphi }_0=0.60$ with 278.7 Kg (46.4 %), Table 15 for ${\varphi }_0=0.75$ with 294.3 Kg (46.6 %), Table 16 for ${\varphi }_0=0.90$ with 309.9 Kg (46.7 %).

Regarding the migration time, after 5 days, a very small but nonzero quantity of contaminant (0.34 Kg, 0.08 %) arrives in the range $2.75<z<3.25$ m, and (17.4 Kg, 4.2 %) in the range $3.25<z<3.75$ for the porosity value ${\varphi }_0=0.10$. See Table 12. For the other porosity values the most extended contaminant migration are situated in the range $3.25<z<3.75$ with 0.093 kg for ${\varphi }_0=0.25$ (Table 13); $3.25<z<3.75$ with 0.0006 kg for ${\varphi }_0=0.43$ (Table 7); $3.75<z<4.25$ with 10.64 kg for ${\varphi }_0=0.60$ (Table 14); $3.75<z<4.25$ with 4.47 kg for ${\varphi }_0=0.75$ (Table 15); $3.25<z<3.75$ with 1.51 kg for ${\varphi }_0=0.90$ (Table 16).

Figure 20 shows the z elevation as a function of the spilled petroleum mass at $(x,y,z)=(0,0,5)$ for $S_w=0.0$, different values of the porosities, ${\varphi }_0=(0.10,0.25,0.60,0.75,0.90)$, and different times. As the porosity increases, the amount of contaminant expelled also increases. This can be clearly seen, for example, in the black line (representing one hour) that runs across all the panels. At the same time, as porosity increases, the contaminant migration velocity decreases, since the upper part empties more quickly for low porosity values, such as ${\varphi }_0=0.10$ and ${\varphi }_0=0.25$, whereas for higher porosity values the mobility of the contaminant decreases drastically. In a certain sense, the low porosity values force the contaminant, which flows out of the pipeline at very high pressure, to find space as it moves downward. Indeed, for ${\varphi }_0=0.10$, after 5 days (see the green line), the contaminant quantity is higher than in the other porosity cases, indicating that the contaminant moves at a higher velocity. For the scenarios with ${\varphi }_0=0.60$, ${\varphi }_0=0.75$, and ${\varphi }_0=0.90$, the contaminant discharge is larger than in the previous cases, but its mobility is compromised by the high viscosity. Indeed, these three values are basically similar.

Figure 21 shows the z elevation as a function of the spilled petroleum mass at $(x,y,z)=(0,0,5)$ for $S_w=0.0$, different values of the time, and different porosities, ${\varphi }_0=(0.10,0.25,0.60,0.75,0.90)$. After one hour (top left panel), it is clearly seen that the biggest quantity of contaminant mass expelled corresponds to the biggest value of the porosity, ${\varphi }_0=\left(0.90\right)$. This mass decreases with decreasing porosity: ${\varphi }_0=(0.10,0.25,0.60,0.75)$. After 5 hours (top-right panel), the contaminant is moving downward due to gravity. The contaminant mass moving more rapidly is the one corresponding to ${\varphi }_0=0.10$ (see the blue line). Similar behavior is confirmed at later times, such as 1 day and 5 days (see bottom panels). That would imply that a low porosity porous medium tends to be contaminated more quickly, although the amount of contaminant is lower than for higher porosity values. This can be used as a recommendation when planning an onshore pipeline: it is better to build an oil pipeline in an area with higher porosity (and/or drier soil) to avoid potential spill migration into the aquifer and reduce water pollution.

3.4. Convergence Study for the Numerical Simulations

The results obtained using CactusHydro and the HRSC method were previously validated in [33,34], using different analytical models, and results from a sandbox experiment [35]. Here, we do not yet have experimental results to compare, since they are not yet available. We have relied on classical convergence tests performed with the same numerical code at different grid resolutions. Figure 22 shows a comparison between results on the saturation contours for petroleum leakage at three different times using two different grid resolutions: $0.25m$ and a time step of $0.00625s$ a),c),e), and $0.50m$ and a time step of $0.0025s$ b),d),f). As can be seen from Figure 22, both results using different times (0.3 hours, 1.0 hour, and 2 days and 8.9 hours are very similar within a $5\%$ error.

4. Conclusions

One of the safest methods of transporting petroleum is through onshore pipelines. However, they are subject to risks such as pipeline corrosion and environmental hazards, which can lead to groundwater pollution and soil poisoning. In this paper, we quantitatively investigate the migration of petroleum in a scenario resembling a broken pipeline (where the pipeline is included in the numerical simulation), and the consequent spread of the contaminant material in the environment.

We use the numerical code CactusHydro, which accurately resolves the sharp discontinuities and temporal dynamics of three-phase fluid flow and thereby accurately tracks contaminant migration in time through a variably saturated zone with variable permeability and porosity. We analyze various scenarios encompassing different climatic conditions, including rainfall-induced water saturation, in the variable-saturation porous medium zone, and different porosity values, for a total of sixteen transient numerical simulations. Hydraulic conductivity is linked to saturation and provides information on the type of porous soil that facilitates contaminant migration. Also, the porosity can convey vertical movement more or less rapidly through the unsaturated zone. The numerical results show that the vertical migration of the petroleum is primarily driven by the high pressure within the pipeline, since the contaminant has a very limited mobility due to its high viscosity. Moreover, hydraulic conductivity and porosity play an important role in petroleum vertical migration. Our results suggest that, for crude oil, the distribution of leaked oil may not extend as deep as two meters in elevation before 10 days.

Our conclusions show that the average seasonal water saturation of the porous medium, its related hydraulic conductivity, and effective porosity must be carefully evaluated when assessing the location of an onshore pipeline and the potential environmental impact of a leak. Our results also suggest that the distribution of leaked petroleum mass for the scenarios investigated may not extend deeper than 2 meters within 10 days. These results may serve as a guideline for constructing an onshore pipeline path that avoids or reduces extreme pollution in a groundwater system.

Our findings on building a sustainable transportation path for petroleum to minimize emergency scenarios include recommendations for constructing a pipeline in a region where a) the aquifer depth is sufficient to protect the groundwater. We find out that constructing a pipeline more than two meters from the groundwater table helps to mechanically remediate the harmful material within a week; b) the characteristics of the soil matrix are important and help to contain the spreading of the crude oil. Indeed, if possible, the choice of constructing a pipeline should prefer a porous medium with lower permeability over one with higher permeability. We have demonstrated that, with a saturation different from zero, crude oil migrates more rapidly downward, as previously reported in other papers using different hydrocarbons [37]. The difference here is that petroleum is highly viscous compared to previous studies, which is good news since its spread is contained. With these recommendations fulfilled, a safer onshore pipeline can be constructed, and in the case of an emergency, less damage is expected.