Simulation of oil spill infiltration and redistribution in a shallow aquifer

Simulations of oil spill infiltration and redistribution in a shallow aquifer were conducted in this work to quantify their effects on groundwater. Results showed that flow of oil spills can be adequately described by infiltration and redistribution patterns within the domain. Without the influence of rainfall, the contaminant plume pervaded the water table within one year of the spill through the influence of dissolution, diffusive and convective mechanisms. The plume advanced in the domain first by domain wall trailing in a manner suggestive of oil-wetting characteristics followed by bulk plume movement. Simultaneous infiltration with rainfall increased the rate of contaminant penetration and aided the plume front further into the domain. Without influence of rain, uncontaminated water can be found at a depth beyond 30 m from the spill surface at 30 years after the spill, while the possibility of clean water with influence of rainfall, at the same period of time, became remote even up to the bottom of the aquifer as the rainfall drives the plume into the furthest depth though, this is accompanied by cleaner water body at the upper portion of the aquifer. It was found that the impacts of contaminants on aquifer water are influenced by contaminant components, wettability of the aquifer, annual rates of rainfall, boundary conditions of the source and possibly size of the domain. Findings in this report point to the needs for quicker clean up exercises and identification of multi-approach procedures for different spill scenarios.


INTRODUCTION
Worldwide, oil spills are recurrent problems with serious consequences both on the living and non-living in the ecosystem. Soil or geological contaminations, owing to spills, are global phenomenon which are eyesores at many sites worldwide. The atmosphere is similarly not spared as volatile and less dense components easily evaporate. Many organic liquids existing as a separate phase are in fact often slightly miscible with water and their solubility can exceed the drinking water standard by orders of magnitude (Van Genuchten, 1980). The dire consequences of these oil spills are, majorly, borne by the host communities where exploration takes place. Also, remote from the exploration sites, the pipe networks for transporting oil to different terminals have further drawn the danger closer to several communities as the cases of ruptured pipes and vandalisation have resulted in casualties from these remote places.
In the Niger Delta area of Nigeria, recurrent cases of spills have drawn international attention. This area is the flood plain area of the country which accounts for over 7.5% of total Nigeria land mass. Kadafa (2012) reported estimates of 9 to 13 million barrels of spill since the beginning of exploration in 1958. Many of these cases have been attributed largely to, pipeline corrosion and faulty tankers. The area consists of highly networked pipelines connecting several small oil fields.
These pipelines are highly prone to rupture as a consequence of poor material structure and maintenance. Vandals and saboteurs are highly involved in this menace as they engage in vandalisation for financial gains or to prove political points by holding government to ransom. Many households, in the exploration areas, whose means of livelihood are based on subsistent farming, have been reduced to paupers due to contamination of their land by persistent spills of oil which destroy the soil fauna. A further and major consequence of the spills is the health impairment due to the presence of the contaminants in the soil, air and drinking water. Court cases have been initiated against the oil-producing companies, at both local and international courts, with some still remaining inconclusive till date.
Nwankwo and Emujakporue (2012), made it clear that even years after the clean-up of oil-contaminated sites, the contaminants proved difficult to eliminate. Problems of water contamination are particularly more serious in Niger Delta areas owing to shallow aquifers with the topmost groundwater levels occurring anywhere between close to the surface and a depth of 10 m 3 . Nwankwo and Emujakporue (2012) carried out electrical resistivity surveys on a spill site in the Niger Delta Area of Nigeria to evaluate subsoil conditions and impacts of contamination on groundwater quality.
This was done as a follow up to the post-spill clean-up exercise. The results showed that, 3 years after clean-up was done, the plume was still strongly present up to a depth of 49 m although they concluded that clean aquifer could still be accessed at an adjacent depth of about 30 m.
This work presents a numerical simulation to investigate the infiltration and redistribution of oil, a non-aqueous phase liquid (NAPL), into the shallow aquifer as found in the Niger Delta area of Nigeria. The Ahoada region of the Niger Delta area of Nigeria has an annual rainfall of about 2400 mm and mean temperature of 27°C (Ehirim et al., 2009). The resistivity values of the sedimentary rock are controlled by parameters such as water contents, salinity, texture, matrix conductivity and the presence of clay materials (Nwankwo and Emujakporue, 2012).
Niger Delta basin is characterised by over 90% sandstone with shale intercalations and thickness of about 2100 m 4 . The domain is composed of coarse gravel and locally fine grained, poorly sorted, sub angular to wellrounded sand samples which give the formation good porosity and effective permeability (Nwankwo and Emujakporue, 2012). Edwards and Santogrossi (1990) described the primary Niger Delta reservoirs as Miocene paralic sandstones with 40% porosity and permeability of 2 Darcy's having 100 mm thickness. However, Molua et al. (2011) found a porosity of 0.22 for the Niger-Delta basin using Wyllie Time-Average equation. Our work made use of intermediate porosity between these reported values.
Niger Delta aquifers are highly vulnerable to oil spills in the geological environment and in area where there is heavy dependence on groundwater for its water supply, Abidoye and Wairagu 547 effects of spill can be catastrophic. Flores et al. (2007) assessed the suitability of image analysis method to measure LNAPL saturation distribution in the fluctuating groundwater. They evaluated the effects of fluctuating groundwater on the oil migration. They used a 1-dimensional column in a small domain (3.5 × 3.5 × 50 cm). In this work, a 1000 × 1 × 50 m aquifer dimension was utilised to simulate the site at Ahoada where the spill history was known.

GOVERNING EQUATIONS
In multiphase flow models, a key element that signifies them all is the functional relationship between relative permeability, k , fluid saturation, S and capillary pressure, c P . Leverett (1941) was among the first to discuss the theory of capillary pressure proposing a Leverett J function attempting to convert all capillary pressure data to a universal curve. In practice, capillary pressure and its relations with fluid-fluid and fluid-solid properties is shown in Equation 1. This is commonly known as the Young Laplace equation.
where, n P is the pressure in non-wetting phase, w P is the pressure in wetting phase,  is the interfacial tension,  is the contact angle and r is the effective pore radius.
Generally, wettability is related to the interface between the solids and fluid or fluid and fluid as shown in Figure 1. The location, flow and distribution of a fluid in porous media is governed by wettability (Anderson, 1987). The wettability of a NAPL-water system is rather complex in nature due to particle size, the disposition of the NAPL as well as the presence of surfactants (Mercer, 1990).
Models describing capillary pressure/saturation relationships have been in existence but the most established ones were those proposed by Brooks and Corey (Brooks, 1964) and Van Genuchten (Van Genuchten, 1980) as seen in Equations 2 and 3, respectively.
where, e S is the effective saturation, w S is the water saturation, rw S is the residual water saturation, d h is the displacement pressure head, c h is the pressure head Figure 1. Contact angle, wetting and non-wetting interface (Mercer, 1990). and  the fitting parameter for the pore size index.
h is the bubbling pressure and n m, ,  are soil parameters.

Relative permeability-saturation relationship
It is inherently difficult to experimentally determine this relationship. This has led to numerous propositions with the popular ones being those proposed by Brooks and Corey (1964) and Van Genuchten (1980). This work utilised the Van Genuchten model as shown in Equations 4 and 5, respectively.
Where, rw k is the aqueous relative permeability, rn k the NAPL relative permeability, w S the effective water saturation and n S the NAPL saturation.

Saturation
Saturation is the fraction of pore space that is engrossed by the fluid phase as represented in Equation 6. (6)

Mass and momentum balances
The main equation for consideration was the continuity equation that is shown in Equation 7.
Where,  is the porosity,  the density ( ML -3 ), S saturation, q the fluid flow velocity (LT -1 ), r k the relative permeability, k the intrinsic permeability ( L 2 ),  the viscosity (ML -1 T -1 ), P the pressure (ML -1 T -2 ), g the magnitude of gravitational acceleration ( LT -2 ) and g e is a unit vector in vertical direction. The mass balance equation and the Van Genuchten model were solved using a multivariable residual based Newton Raphson iteration technique with maximum number of Newton iterations set at 16 and convergence at 10 -6 .

MODEL DESCRIPTION
The subsurface transport over multiple phases (White and Oostrom, 1996) was used for simulation in this paper. It was chosen due to its versatility in multiphase flows and has been successfully validated and applied in the simulation of a wide range of multi-phase flow problems (White and Oostrom, 1996;Schroth et al., 1998;Ataie-Ashtiani et al., 2001;Das et al., 2004). It is composed of FORTRAN subroutines that solve the governing equations. A series of output files were postprocessed with Perl scripts and Tecplot was finally used to generate the plots.

Simulation domain
The non-hysteretic two phase version of water-oil mode of STOMP was used for simulation purposes where the governing equations are discretised into algebraic form using the integrated finite difference method (Patankar et al., 1980). A 2-D Cartesian coordinate system was used to identify the model domain of the Niger Delta region being studied. The 1000 × 1 × 50 m domain was divided into a grid containing 25 and 20 grid blocks in the x and z direction, respectively. Figure 2 shows the 2-D domain of the modelled aquifer. An automatic time step that was dependent on relative component property changes was used to discretise the total time step.   (White and Oostrom, 1996) Residual water, (Srw) 0.1 (White and Oostrom, 1996) Properties for homogenous soil system Generally, density, permeability, porosity and saturation are input soil parameters in the STOMP application. Table 1 lists the properties adapted for the current numerical simulation. It is worth noting that a constant porosity of 0.37 was used for the entire simulation which was within the average values reported for the Niger Delta region (Nwankwo and Emujakporue, 2012; Molua et al., 2011).

Fluid properties
Data on the properties of the various fluids involved in the contaminant transport is required for the modelling of LNAPL flow. The fluid properties are listed in Table 2.

Initial and boundary conditions
In the Niger Delta, majority of the water table is found at around 10 m from the surface (UNEP, 2011), thus the simulation was conducted by setting the water table at the   Abidoye and Wairagu  549 value. Subsequently, aquifer thickness was set at 40 m and was used to calculate the initial boundary NAPL and aqueous pressures with the addition of atmospheric pressure to a value of 493725 Pa. The boundary conditions were set up as follows. At the top of the domain, a constant atmospheric pressure was applied. The spill zone was along the entire x direction at the top boundary. To conserve the mass of the oil and make its movement finite within the domain, the west and east boundaries were designated as zero flux for NAPL. A table of boundary conditions used for the simulation is shown in Table 3.

Infiltration without influence of rain
At the beginning of the simulation with no oil saturation present, oil infiltration began at all the specified surface nodes as specified in input simulation card. After about one year of infiltration it was observed that the oil front was on the surface of the water table. A similar trend was observed by Nwankwo and Emujakporue (2012). The saturation at this instance is below 0.4 and the highest point of concentration lies behind the front plume as seen in Figure 3. It can also be seen that the oil plume travelled in phases ranging from a very light zone at the front phase and the thickest zone close to the spill source. These layers were thought to be of varying oil properties such as density and viscosity as well as relative permeability due to the different components of oil. As a result, each component travels at a different rate.
It is widely acknowledged that infiltration and redistribution at spill sites occurs over a long period of time. In the current simulation, it was seen that at one year of infiltration, the oil was present in the water table. This can be understood from the nature of aquifers in the Niger Delta area of Nigeria. They are typically shallow whereby boreholes and wells are able to obtain fresh water within a few meters of digging (UNEP, 2011).
In the present simulation, it was shown that within a short period of infiltration, the oil was able to invade the water table. This would impact negatively on the suitability of the water for human consumption. It can be seen that the concentration of oil close to the water table was about 0.1. The more concentrated front was observed close to the spill surface as infiltration proceeded.
The diffusive nature of the plume is discernible as the concentration gradient progresses outwardly from the infiltration zone. While lighter zone components can be found below the water table, the water table retards the movement of the heavier insoluble constituents. However, it was later shown that the diffusive behaviour of the plume broke resistance to penetration as lighter or less viscous components separated out. Five years after

Property
Value Molecular weight of oil (g mol   -1   ) 165.83 (White and Oostrom, 1996) Oil density, ρ (kg m -3 ) 830 (Isehunwa and Olubukola, 2012;1996) Oil viscosity, µ (Pa s) 0.002 (Isehunwa and Olubukola, 2012;1996) Critical compressibility of oil, bar 0.28 (White and Oostrom, 1996) Critical pressure of oil, bar 47.6 (White and Oostrom, 1996) Boiling temperature of oil, K 394 (White and Oostrom, 1996) Freezing temperature of oil, K 251 (White and Oostrom, 1996)   the spill, the water table level was observed to be fully contaminated and further infiltration into the water table level was visible as seen in Figure 4. This would suggest typical behaviour in shallow aquifer regions such as the Niger Delta of Nigeria. However, this behaviour might be different in deep lying aquifers. Clear evidence of the influence of time on the percolation mechanism is evident as the 1 year of infiltration showed the oil plume below the 10 m depth from surface as seen in Figure 3. In contrast, Figure 4 showed the plume to be a few me-  ters below the water level table mark. As redistribution proceeds, it is observed that the penetration force of the advancing plume starts reducing as is seen in Figure 5 10 years into the spill. This is due to oil infiltration at the surface of the spill stopping after one year. However, the fringe redistribution in the domain showed the tendency of the plume to trail the domain wall at the east and west ends which precedes the bulk plume movement into the water body. This could be as a result of the changing wettability of the medium which tends to be more oil-wetting or weak-wetting as oil saturation Abidoye and Wairagu 551 increases. It was seen that there was a marked drastic movement of the plume in the domain after 20 years of simulation as shown in Figure 6. The movement of the plume at the wall boundaries was seen to be moving at a faster rate into the domain. This resulted in more bulk oil penetrating into the water body. It is noteworthy that the concentration of the oil in these areas had become weaker owing to the discontinuity of infiltration of oil at the surface. This resulted in limited clean water 20 m below the surface of the oil after 30 years of the spill. It was seen that oil infiltration advanced slowly and the concentration reduced as seen in Figure 7. The wall boundaries of the plume still continued advancing thus showing the possibility of oil wetting conditions or a weak wetted medium.
Nwankwo and Emujakporue (2012) reported three major resistivity zones in their domain; low, medium and high as seen in Figure 8. They considered the low resistive zones as uncontaminated. Contamination was found as far as thirty meters below the surface though clean water was still available ten to fifteen metres towards the west of the aquifer. This was three years after the actual spill as reported in their paper. In the same line, the simulations carried out in this paper highlight the danger of delayed remediation which would result in contaminant presence within the entire water body over a period of time.

Influence of rainfall on redistribution
In investigating the behaviour and effects of pollutants,  one should also investigate the effects of certain environmental conditions particular on the area of study which have direct and indirect influence on the behaviour and effects. The Niger Delta is reported to have a mean rainfall of 2400 mm per annum (Ehirim et al., 2009). Influence of rainfall was investigated on the infiltration and redistribution of the oil plume. This condition was implemented in our simulation.
It was observed that five years after the spill and with inclusion of rainfall, the movement of the plume within the domain was more pronounced as compared to dry conditions. For the same length of time, that is, 5 years, under dry conditions ( Figure 4) the plume was at about 10 m from the spill surface while under the influence of rainfall (Figure 9), the plume was already at about 30 m from the same surface. The concentration of the oil had thinned out from 0.1 in Figure 4, to 0.001 in Figure 9. This was as a result of advanced dispersion of the oil phase by the rainfall. Concentration distributions were observed within the domain as shown in Figure 9. The region closest to the spill surface remained light and uniform in concentration  at about 0.09 while the region immediately following this, shown with yellow band has concentration in the region of 0.1. The thickest region of the plume had a concentration of about 0.13. Thus, it was deduced that the concentration increased away from the spill surface. The region closest to the spill surface remained lighter owing to the discontinuity of oil infiltration after a year. Meanwhile, the least oil-saturated region of the plume occurs at the front where the oil body interfaces with the fresh water body. This is thought to occur due to the rapid Abidoye and Wairagu 553 dissolution of the oil components at this interface or differential speed of oil components as it enters different medium. This behaviour can also be explained to arise from the effect of convective transport in which the incoming rainfall transports the oil components further into the domain. The viscous forces could be higher for the light components of the oil thus they are transported further into the domain leaving behind the more viscous samples which tend to be evenly distributed throughout the top portion of the aquifer.
The significant influence of rainfall was more evident 10 years after the spill. Redistribution of the oil was under the influence of rainfall as seen in Figure 10. The layer closest to the spill surface had further thinned out owing to the washing off by rainfall into the domain. This led to aggregation which increased the concentration of the oil plume further into the domain. It was also observed that aggregates of oil were found at the bottom of the domain close to the bottom boundary surfaces with further advancement of the main plume in comparison with Figure 5. Rainfall was shown to speed up the contamination within the domain. At later simulations times, it was observed that rainfall led to further infiltration within the domain as seen in Figure 11. The oil plume had advanced to the bottom of the domain as compared to Figure 7 at the same thirty years period. Thus, it is overwhelmingly evident that rainfall speeds up aquifer contamination making the rain-prone areas more susceptible.

CONCLUSION
Simulations of oil spill infiltration and redistribution in a shallow aquifer were conducted in this work to quantify their effects on groundwater. A case study of the Niger Delta in Nigeria was used to obtain the input conditions for the simulation. Results showed that flow of oil spills can be adequately described by infiltration and redistribution patterns within the domain.
Without influence of rainfall, the infiltration profile showed that the penetrating plume of oil had moved up to the water table within one year of active penetration with average NAPL saturation of about 0.3 in the unsaturated portion of the domain. The redistribution of oil contaminants, at the end of infiltration, was mainly characterised by domain wall-plume movement which later drew the bulk fluid of the plume further down into the domain. Bulk oil was found at around 20 m from the spill surface 30 years after the spill.
The plume was observed to be thinning out after the stoppage of infiltration. This showed that redistribution continues long after the infiltration stops. This served to carry the contaminants further into the water body. It was also observed that the medium showed oil-wetting or weak-wetting condition as evidenced with the plume trailing along the domain walls which continued to look like norm before the bulk plume advanced. Experimental work is needed to investigate this further and could point direction on the control of wettability to check contaminant progress in an aquifer.
Rainfall was shown to drive the contaminant plume further into the water zone. The interpretation of this is that contaminant transport was faster with rainfall when compared with the infiltration without rain at the same length of time.
Thus, the impacts of contaminants on aquifer water are influenced by contaminant components, wettability of the aquifer, annual rates of rainfall, boundary conditions of the source and possibly size of the domain.