SIMULATION OF NAPL VERTICAL INFILTRATION IN A
HETEROGENEOUS SOIL
1Boris N. Chetverushkin*, Natalia G. Churbanova*, Nikolaj V. Isupov† and
Marina A. Trapeznikova*
* Institute for Mathematical Modeling RAS 4-A Miusskaya Square, Moscow 125047, Russia
e-mail: marina@imamod.ru
† Moscow Institute of Physics and Technology
9, Institutskii per., Dolgoprudny, Moscow Region, 141700,Russia
e-mail: n_isupov@7ka.mipt.ru
Key words: Multiphase Flow, Light and Dense Non-Aqueous Phase Liquids, Infiltration,
Heterogeneous Porous Media
Abstract. The problem of numerical simulation of contaminant infiltration and spreading in
heterogeneous soils consisting of layers with different porosities and permeabilities is investigated. As contaminants Light and Dense Non-Aqueous Phase Liquids are considered. The governing model takes into account capillary and gravity forces and includes a capillary pressure – saturation relation. Some interface conditions are provided on boundaries between different materials. A number of tests are predicted.
1 INTRODUCTION
The research deals with numerical simulation of contaminant infiltration and spreading in soils. The problem has substantial ecological sense to predict possible contamination of the soil air and the groundwater. As contaminants mineral fuels, solvents and detergents (in other words Non-Aqueous Phase Liquids - NAPLs) can be considered. There are Light and Dense NAPLs depending on the relation of their density to the water density. For example, petrol and fuel oil are LNAPLs but tetrachloroethylene is DNAPL. These liquids are immiscible with gas and water therefore one has to treat multiphase flows. Governing models take into account capillary and gravity forces and include some capillary pressure – saturation relations1, 2.
A number of test problems were solved. Tests concerning the DNAPL infiltration into a fully saturated reservoir and the LNAPL infiltration into an unsaturated reservoir require the
two-phase problem statement. In the present work this statement is transformed to the phase pressure – saturation formulations2, 3, numerically implemented by the finite differences with consequent solution by the IMplicit Pressure – Explicit Saturation (IMPES) method1. Algorithms developed by the authors earlier4, 5 are extended for heterogeneous media. In the future the authors intend to develop an efficient fully implicit algorithm to avoid the strong time-step restriction.
All the problems are treated in the 2D vertical section. It is assumed that porous media are heterogeneous, i.e. they consist of layers with different porosities and permeabilities or contain low-permeable lenses. As there are saturation discontinuities on boundaries between different materials special interface conditions3 should be realized to ensure physically correct solutions. The goal of analysis is investigation of the time-evolution of the contamination domain.
The applied problem that will be under consideration in the nearest future is connected with the petrol infiltration from the earth surface into the relatively dry ground until reaching the water table. Realistic data obtaining by measuring in the locality will be used. Investigation of the oil-products propagation is of great practical importance due to the potential danger of drinking water contamination. In this case the three-phase flow has to be considered. Software tools developed and tested up to the present time will be implemented for solving such kind of ecological problems.
2 MATHEMATICAL MODEL
The process of NAPL infiltration into the ground happens under the influence of the gravitation. Then the Darcy velocity is low therefore the capillarity can not be neglected. Let us consider in these conditions incompressible two-phase flow of wetting and non-wetting fluids: for example, the coupled flow of LNAPL and soil air or the flow of water and DNAPL. Such flow is governed by the given below system of equations taking into account capillary and gravity forces (as usual index 1 corresponds to the wetting phase and index 2
– to the non-wetting phase):
i= - i
( ) (
i i)
, 1,2,( )
0,1 , 1(
here and everywhere below)
i k s k P ρ g i s s µ − = = = W grad e e (1)
( )
div , 1,2 i i i s m q i t ∂ + = = ∂ W (2) 1 2 1 s + =s (3)( )
2 1 c P − =P P s (4)The next notations are used: si is the phase saturation, Pi is the phase pressure, Pc is the
is the gravitational acceleration vector, qi is the source term. Note that in general the absolute
permeability is a tensor but in the present paper it is considered as a function of coordinates, the anisotropy is not taken into account yet.
In addition to equations (1)-(4) the following relationship for the relative permeability must be satisfied:
( )
, =1,2 =k s iki i (5)
It should be noted that the relationships for the capillary pressure and the relative permeability are strongly non-linear functions of the saturation. An analytical determination of the capillary pressure – saturation relation is impossible because of the irregular pore geometry. The most famous correlations fitted to experimental data are models by Leverett1, Brooks & Corey and Van Genuchten2, 3. The relative permeability functions can be derived using a connection between the capillary pressure – saturation relation and the relative permeability.
The general form of the two-phase flow equations (1)-(5) transforms to different mathematical formulations depending on the individual problem setup. In the current investigation the phase pressure – saturation formulation2, 3 is employed. One of the pressures and one of the saturations are eliminated using algebraic constraints. As a result one has to solve the system of two partial differential equations - the elliptic equation for the phase pressure and the transport equation for the saturation.
3 TEST PREDICTIONS
3.1 LNAPL infiltration into an unsaturated reservoir
As one of test problems LNAPL (for example, oil) infiltration and propagation in the unsaturated zone consisting of layers with different absolute permeabilities are predicted (see Figure 1).
Here the capillary pressure is presented by the Leverett model with the next approximation of the Leverett function J(s):
( )
cos( )
c m P s J s k σ θ = (6)( )
0 0 0 1 , 0 1 , 0, 1 0.1, 0.1, 0.9 0.05, cos 1 s J s s s s J s J s s s s s s J s s σ θ × − ≤ ≤ = × − < ≤ < ≤ = = = = = (7)Relative permeabilities are defined by the following relations:
( )
( )
(
)
3.5 3.5 1 2 0, 0 1 3 , 0.1 , 1 0, 1 s s s s s s s k s s s k s s s s s s s ≤ ≤ − × + ≤ < = − = < ≤ ≤ ≤ (8)Everywheres is the residual saturation and s is the critical saturation. Initial conditions:
( )0 ( )0
2 2
, atm
s = s P =P +ρ gy (9)
Boundary conditions can be written like this (in notations of Figure 1):
1
Γ − the principal plane:
1 2 1 2 1 2 0 0 0 x x P P s s W W x x x x ∂ ∂ ∂ ∂ = = ⇒ = = ⇒ = = ∂ ∂ ∂ ∂ (10) 2
Γ − the impenetrable bed:
1 2 1 2 0 , 0 y y i P P s W W g y y ρ y ∂ ∂ ∂ = = ⇒ = = = ∂ ∂ ∂ (11) 3
2 atm 2 , 1
P =P +ρ gy P =const (12) 4
Γ − the earth surface:
2 1 1 1 , 2 atm, '( ) k P g P g P P s y y y P s ρ ρ ∂ − ∂ = = ∂ = ∂ ∂ ∂ (13) 5
Γ − the NAPL spot surface: a) till the spot exists:
2 atm 1 ( ),
P =P +ρgH t s s= (14)
b) after the complete infiltration ( ( ) 0)H t ≤ - the same condition as at Γ 4 The spot height is calculated using the next relations:
2 2 1 1 0 0 ( ( )) H R init R H H t m s dxdy ρ × × − =ρ
∫ ∫
(15) 2 2 1 0 0 1 , : ( ) 0 ( ) 0 : ( ) 0 H R init H m s dxdy t H t H t R t H t − > = ≤ ∫ ∫
(16)Thus, the NAPL spot has a fixed volume and can spread from the surface into the underground completely. It has an effect on boundary conditions. The condition on Γ5 is being changed during computations (see (14), (13)).
For computations the above general system of equations is reduced to the (P2, s1
)-formulation. Remember that s1 is the wetting (NAPL) phase saturation, P2 is the non-wetting
(soil air) phase pressure. Let us denote s = s1, P = P2.
( )
( )
( )
( )
( )
( )
1 2 1 1 2 1 1 2 1 2 1 2 div kk s kk s P div kk s P sc k s k s g k g k y y µ µ µ ρ ρ µ µ − + = − + ∂ ∂ + + ∂ ∂ grad grad (17)( )
( )
( )
( )
1 1 1 1 1 1 1 div k s c div k s k s s m k P s k P g k t µ µ ρ y µ ∂ + = + ∂ ∂ grad grad ∂ (18)Equations (17)-(18) are approximated by finite difference schemes over rectangular grids of the MAC type and implemented by the IMPES method1. Numerical algorithms developed by the authors earlier for solution of oil recovery problems4, 5 are mainly used. These are implicit red-black ad hoc SOR (the local relaxation method) for equation (17) and the explicit method with upwind differences for equation (18).
Results obtained in the case of a homogeneous medium when the absolute permeability is constant (k = 10-12 m2) are depicted in Figures 2. Results obtained in the case of a heterogeneous medium consisting of three layers with absolute permeabilities equal 10-13 m2, 10-14 m2 and 10-12 m2 correspondently are shown in Figures 3. In computations the next values of constants were used: R1 = 2.5 m,R2 = 10 m, Hinit = 10 cm, m = 0.1, ρ1 = 103 kg/m3 ,
ρ2 = 1 kg/m3, µ1 = 5×10-4 kg/(m·s), µ2 = 1.84×10-5 kg/(m·s).
Figure 2: LNAPL saturation in a homogeneous medium at the final moment.
Figure 3: LNAPL saturation in a heterogeneous medium at the final moment.
means that the entry pressure equals zero and the NAPL penetrates to the lower permeable layer not damming up above the interface.
3.2 DNAPL infiltration into a fully saturated reservoir
Another test problem treated by the authors is the DNAPL (for example, tetrachloroethylene) infiltration into a reservoir filled with two kinds of sand and saturated with water. A low permeability lens (fine sand) is placed into its interior (see Figure 4). The material properties and the model parameters were taken from works by R. Helmig et al 2, 6.
It should be noted that DNAPL is a non-wetting phase with respect to water in contrast to the previous test problem where LNAPL is the wetting phase with respect to air. Therefore the (P1, s2)-formulation of the model is used, i. e. the primary variables are the wetting (water)
phase pressure and the non-wetting (NAPL) phase saturation.
In this investigation the capillary pressure – saturation relation by Brooks and Corey2, 3 are chosen: λ 1 − = d e c PS P (19)
where S is the effective saturation seems like this: e
s s s Se − − = 1 (20) d
P is the entry pressure - the capillary pressure required for a non-wetting fluid phase to displace a wetting fluid phase inside the corresponding sand, λ is the pore size distribution index of the sand.
The corresponding relative permeability functions after Brooks and Corey are expressed as follows:
( )
(2 3λ)/λ 1 + = e e S S k (21)( ) (
)
2(
(2 λ)/λ)
2 1 1 + − − = e e e S S S k (22)where S is the effective saturation (20). e
of the first sand is equal to the entry pressure of the second sand. In the term of the water saturation this interface condition can be written as follows:
[ ]
< ≥ = − Γ Γ Γ * 1 1 * 1 2 s if if 1 2 G S P S G s G s G c (23) Here[ ]
1 2 − G cp is the inverse function of
(
Pc − relation. s)
For accurate numerical implementation of interface conditions the computational grid nodes should fall exactly on the interface. The same algorithm as for solution of the LNAPL infiltration problem is used.
Figure 4: Setup of the problem of DNAPL infiltration into the medium including a low permeability lens.
Figures 5, 6 demonstrate the obtained time evolution of the DNAPL saturation field. As evident from these Figures DNAPL can not penetrate inside the lens under the given conditions. After reaching the interface (see Figure 5) DNAPL spreads in the horizontal direction and flows around the lens (see Figure 6) what agrees with experimental data.
Figure 6: DNAPL saturation in the medium with a low permeability lens at the moment of flowing around the lens.
5 CONCLUSIONS
- The problem of NAPL vertical infiltration and propagation in saturated and unsaturated heterogeneous porous media has been numerically investigated.
- Algorithms for accurate implementation of the transition across the interface between zones with different properties have been developed.
- In the future the algorithms will be extended for solution of an applied ecological problem concerning three-phase flow in a heterogeneous soil.
- An efficient fully implicit algorithm will be developed to avoid the strong time-step restriction.
REFERENCES
[1] K. Aziz and A. Settari, Petroleum reservoir simulation, Applied Science Publ., Lmt.. London, (1979).
[2] R. Helmig, Multiphase flow and transport processes in the subsurface – A contribution to the modelling of hydrosystems, Berlin, Springer, (1997).
[3] P. Bastian, Numerical computation of multiphase flows in porous media, Habilitation thesis, Christian-Albrechts-Universitaet Kiel, (1999).
parallel computers with distributed memory”, In: CFD’98, Proc. of the 4th European CFD Conf., Wiley, Chichester, Vol. 1, Pt. 2, pp. 929-934, (1998).
[5] M.A. Trapeznikova, N.G. Churbanova and B.N. Chetverushkin, “Parallel elliptic solvers and their application to oil recovery simulation”, In: HPC‘2000, Grand Challenges in Computer Simulation, Proc. of Simulation Multi Conf., SCS, San Diego, CA, pp. 213-218, (2000).
