- Email: [email protected]
Multi-phase flow modeling of air sparging
Advances in Wukr Resources, Vol. 18, NO. 6, pp. 319-333, 1995 Copyright 0 1995 Elsevier Science Limited Printed in Great Britain. All rights reserved 0309-1708/95 $09.50 + 0.00
0309-1708(95)00028-3
ELSEVIER
Multi-phase flow modeling of air sparging M. I. J. van Dijke, S. E. A. T. M. van der Zee Department of Soil Science and Plant Nutrition, Wageningen Agricultural University, PO Box 8005, 6700 EC Wageningen, The Netherlantis &
C. J. van Duijn Department of Mathematics, De@ University of Technology, PO Box 5031, 2600 GA De@, The Netherlands (Received 14 February
1995; accepted 25 July 1995)
Air injection into groundwater (air sparging) in a homogeneous axially symmetric porous medium is modeled using a two-phase flow approach. A numerical method based on the mixed form of the Richards equation for both phases is presented. Furthermore two analytical approximations are discussed to explain the numerical results. One is a one-dimensional description explaining the occurrence of small air saturations. The other is a closed form approximation for the distribution of the air saturation in the resulting steady state. From the latter we can estimate the maximum radius of influence of air sparging, as a function of the physical parameters. The analytical approximation at steady state and the numerical results are in good agreement. Key words: groundwater, air sparging, multi-phase influence, analytical solution.
NOTATION
Hu (4,)
Dimensionless variables between brackets ( ). . Subscript denoting air. Dimensionless filter surface. D Dimensionless diffusion function for two-phase flow dependent on S. Dimensionless diffusion function for b single-phase flow dependent on k. Unit vector lin z-direction. Dimensionless position of free boundary. z-Component of flux function. Flux function. (Dimensionless) position of free boundary at water tab1.e [ml. Air fractional flow function. h Gravity [m s2]. g H Filter depth [ml. (Dimensionless) depth of bottom boundary Hb (hd
k K abs
Row, steady state, radius of
(Dimensionless) depth of upper side of filter
[ml.
1
HI (hJ
k]* rmensionless) depth of lower side of filter
Ht (ht)
E]* rmensionless) height of top boundary
krj A4
:l,
N2
N, Njg Ng pc Pjt
(PA
(Pjt>
Pt
Q R (4 & kb) S
si Sr St
[ml. 319
Redefined air relative permeability. Absolute permeability [m’]. Fluid j relative permeability. Mobility ratio. van Genuchten parameter. Dimensionless numbers determining steady-state situation. Capillary number. Fluid j gravity number Gravity number. (Dimensionless) capillary pressure [Pa]. (Dimensionless) fluid j pressure at top boundary pa]. Dimensionless capillary pressure at top boundary. Air injection rate [m3 s-t]. (Dimensionless) horizontal coordinate [ml. (Dimensionless) width of right boundary [ml. Redefined air saturation. Fluid i saturation. Residual water saturation. Air threshold saturation.
320 T, (0 uj <$> uin 4 V l (4
i(f) x & Pj 4
M. I. J. van Dijke et al.
(Dimensionless) time [s]. (Dimensionless) fluid j Darcy velocity [m s-l]. Air injection velocity [m s-l]. Dimensionless total fluid velocity. Dimensionless air volume in flow domain. Subscript denoting water. (Dimensionless) vertical coordinate [ml. van Genuchten parameter. (Dimensionless) filter radius [m]. Mobility function. Fluid j viscosity [Pa s]. Fluid j density [kg me3]. Porosity.
1 INTRODUCTION A method for remediating an aquifer, which is contaminated by organic liquids (solvents, gasoline) trapped in the saturated zone, is to inject air or oxygen into the aquifer. Injection of air may enhance microbial degradation and volatilization. For remediation of the unsaturated zone, gas venting has been studied both experimentally and numerically.637’11Besides the use of soil vapor extraction wells, injection of air at a small distance below the water table to remove contaminants at groundwater level has been considered.7 For remediation of the saturated zone the use of a vacuum vaporizer well has been described.g Emphasis was laid on assessing the sphere of influence of the water circulation, where dissolved oxygen enhances microbial degradation. Direct injection of air in the saturated zone, known as air sparging, together with vapor extraction in the unsaturated zone, has been put forward as an effective in-situ remediation technique. Air sparging has been studied experimentally on both field3>i5and laboratory’0’22 scales in the last few years, with emphasis on the region in the saturated zone where air is present (radius of influence). A numerical study of air sparging at steady state has been provided for the case of injection just below the original water table.16 Although laboratory experiments with air injection in a glass bead medium” showed that air flow can occur as isolated bubbles for bead diameters exceeding 2 mm, it is assumed that under natural subsurface conditions air flow is most likely to occur in small continuous channels.‘0112 A recent field study for a uniform soil with mean grain size of 0.25mm,i4 indicated that the density of channels in the main region of air flow must be very high. Therefore, it is reasonable to model air flow macroscopically as a continuum. Air sparging is only possible in relatively coarsegrained soils. The reported injection rates varied from about 3 m3 h-’ up to 100m3 h-’ for thick sandygravelly deposits.3112 Injection in fine textured soils
requires high air entry pressures, which may cause fracturing of the soil. This may result in the formation of a few channels through which air flows upwards.” Hence, knowledge of soil layering and heterogeneity is important because these affect the radius of influence.12 Besides flow continuity, air-phase compressibility may be a complicating factor in modeling air sparging. In the early stage of sparging when air pathways from the injection filter to the unsaturated zone have not yet been established, air density is not constant. However, emphasizing the steady-state situation in which continuous channels to the vadose zone exist, compressibility is expected to play a minor role. This expectation is supported by the field study,14 which describes sparging at a rate of 34.2 m3 h-l, resulting in a pressure increase of 10.5 kPa at a distance of 0.6m from the filter. This causes an air density increase of about 10% only. Important for dimensioning the technique of air sparging is quantitative knowledge of the effect of soil, fluid and filter parameters on the radius of influence of an injection filter. This motivated us to study an axially symmetric model for air injection through a vertically positioned filter in an initially saturated region below the vadose zone. In the model we consider air and water as two immiscible incompressible continuous phases. The interaction between the fluid phases and the soil matrix is described by the saturation dependent relative permeability and capillary pressure functions. Contaminants are assumed to be part of the soil matrix and do not affect the flow process. The basic equations are given and discussed in Section 2, where they are reformulated in terms of dimensionless numbers (that are combinations of the physical parameters). In Section 3 we discuss a numerical method, which is based on the mixed form of the Richards equation for both water and air. Results of numerical computations are presented in terms of the distribution of air saturations and the volume of air that is stored in the domain. In Section 4 we discuss two analytical approximations which are valid in the relevant part of the flow domain. One is an upper bound that explains the occurrence of small saturations. The other is an explicit solution for the steady-state situation, which is derived under the assumptions that water is immobile and that gravity is the dominant effect in the vertical direction. The explicit solution provides an expression for the radius of influence of sparging. In Section 5 we quantify the agreement between the analytical and numerical solutions at steady state.
2 MODEL We use Darcy’s Law for both air (a) and water (w) q = _5!E+(Pj J
+pjg
j= w,a
(2.1)
Multi-phase flow modeling of air sparging
and the mass balance e:quations aSj + $,t+V’Uj=O
the water table is situated at Z = 0. Hence,
j=w,a
(2.2)
In eqn (2.2) Sj denotes the effective fluid saturation and C#J the effective porosity, i.e. redefined according to s
SW-sr
=
w:
&
Sll := 1 _ $7
1-S’
C$:= (1 - S,)C$
where S, is the residual water saturation. Furthermore, Kabs is soil absolute permeability, CJ~soil porosity, Cj fluid Darcy velocity, Sj fluid saturation, krj fluid relative permeability, /Lj fluid viscosity, Pj fluid pressure, pj fluid density and g gravity. We assume that the soil is homogeneous and isotropic and that both fluids are incompressible. The set of equations., (2.1) and (2.2), is completed by the constitutive relati’ons S,,, + S, = 1, the capillary pressure PC = P, - P,,,, Sj = Sj(Pc) and krj = _k,i(Sj) and is solved for the unknowns S,,,, S,, P,, P,, U, and CO. The dependence of capillary pressure and relative permeability on the saturations is given by the wellknown expressions17 p&y,)
= +,lirn
- ,)I-,
Pj = Pj,
for j = w, a, where Pwr = -p,gHt and Pat = 0. Air is injected with velocity Vi, into this domain through a filter with radius E > 0 and length HI - H,, giving the t_otal injection rate Q = 2rE(Hl - HJUi,. Writing uj = (“j,r, V,,z), we have the flux boundary conditions UUJ,,= 0 uo,r =
(24
k&J
= Sf,(l - (1 - S0)1’m)2m
(2.5)
Equations (2.1) and (2.2) are solved in the twodimensional axially symmetric domain of Fig. 1. The level Z = 0 corresponds to the initial position of the water table, where P, := 0. Along the top boundary of the domain, Z = H,, the air pressure equals P, = 0. The water pressure along this boundary is selected such that
=
Pj = Pjt \
pi = pjr + 9
c
(h, .-
Z)
-c Z < Ht
Pj = Pjt + pjg(Hf - Z) for r > E, Z < Ht
(2.8)
We introduce the dimensionless variables
TUin
_ 3j uj=Y,y
R
r=-
t=-@,
H’
z=-
Z H
(2.9)
0Pj Pj = Pg w
j=w,a
h,,+‘,
/q+
h$,
pit=-- ffpjf P&T
The flow problem is now determined by the mobility ratio, the gravity numbers for water and air, the capillary number and the dimensionless filter surface, i.e.
&&CL,,
N,
Pa N
C
=
18
KabsPwg
=
~a Ui,Ha ’
KabsPjg Pa uin
(2.10)
Q A=Ui,H*
For both water and air, eqns (2.1) and (2.2) are combined into the mixed form of the Richards equation.41’3 Together with the boundary and initial conditions the resulting problem is given in (2.11)
t>O
+ N,,=))
for r = E, z < -hl,
0
-H,
The initial conditions at T = 0, say, are
for r = E, -h, < z < -h,, Uj,l
forR = E, Z < -H,,
(2.7)
in r > E, z < h,,
<
for R = E, -HI < Z < -H,, 1
Uj,, = 0
13,17 C’
= V - (k,V(Ncpa
uin
E E=---, H
= St,(l - (1 - Sy”)“)*
at
(24
(2.3)
k&Y,,,)
-
for R > E, Z = Ht
where H = (H, + Hl)/2, the depth of the filter center, is chosen as a characteristic length. Furthermore, we define the dimensionless constants
where 0 < m < 1 and (Y> 0 are given constants. Note that S,,,(P,) is dehned for PC 2 0. For P, < 0 we set S, = S,(O) = 1. Hence, 0 5 S, 5 1 for all values of p
321
-h,
t>O
< z < h,,
t> 0
for r > E, z = h,,
t>O
for r > E, z < h,,
t=O
(2.11)
M. I. J. van Dijke et al.
322
where
Mkr=
f, =
(2.14)
kw + Mkra denotes the air fractional flow function and
(2.15) denotes the mobility function. Ng = Nwg - Nag and e’, is the unit vector in the z-direction. Further, i& = u’, + il, denotes the total velocity, which satisfies the incompressibility condition v.u’,=o
Observe that in case gravitational and capillary effects are neglected, the air velocity ii, is the fraction f,of the total velocity. Setting S = S, we arrive at the initial-boundary value problem (2.17)
Fig. 1. Schematic of the domain (s2)for air sparging.
dt
+
V
(faz-it + N&”
.
f, - NJ%
= 1
- NJVp,)
(2.16)
= 0
for r > E, z < h,,
t>O
for r = E, -h, < z < -h,,
t>O
(2.17)
24t,r -1 -
ds 0 ar= U
for r = c, z < -h,,
-h,
< z < h,,
t>O
1,r -0 -
PC =Pt , Pc=P+-
N (h, - ‘) c
for r > C, z = h,,
t>O
for r > E, z < h,,
t=O
where pwr = -(N,,/N,)h, and par = 0. The dimensionless capillary pressure, which depends now on S,, is, according to the redefinition of the phase pressures pc(Sa) = ((1 - s,)-“m
- 1)1-m
(2.12)
with S, = 0 for pc < 0. Problem (2.11) is solved for one of the phase saturations and one of the phase pressures. In general, however, it cannot be solved explicitly due to the nonlinear nature of the equations. To develop an understanding of the behaviour of the solution, we consider a numerical solution technique based on the formulation of problem (2.11). In addition, we characterize some aspects of the solutions analytically. To achieve this we construct a single equation for the air saturation only. We find in terms of Sa:2T5
ii, = f,i& + N,Xe’, - NJV,,
(2.13)
where pr =pot -pwt = (N,,/N,)h, = cd!Zh,. Note that this problem is equivalent to the original description (2.11) in terms of the separate phase saturations and pressures. However, to understand the qualitative behaviour of the solution, one often considers iit in (2.17) as given, leading to a nonlinear advectiondiffusion equation in terms of S only, which is more accessible to analysis. In problem (2.17), f,(S), X(S) and p,(S) are nonlinear functions of S. These nonlinearities are grouped in the functions &S,
i&) = f,(S)&
+ N,X(S)e’,
(2.18) D(S)
= N,X(S)
g
(S)
In Fig. 2 we show some typical examples of (2.18). We note here that the r-component of the flux, F,, is increasing in S for all z.+ > 0. The z-component, F,,
Multi-phase flow modeling of air sparging 0.14 0.12 O.1o
323
0.030 ,_----__ F, --__ --__ I’ ___;; . . . . . . . . . . . . . . . . . . .
0.08
/: , : : :
0.06
:
j
!A
0.020
n
:
:
: :
: :
0.02
:’
~
0.00
--------_____
Fr
0.04
0.00 EIl”’
0.025
---_____
0.015
0.010
0.005
jst -I 0.06
0.12
0.18
0.24
0.30
0.36
0.42
S
(a)
0300
0.48
L 0.00
0.10
04
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
S
Fig. 2. Graphs of (a) the radial component F, of F’for z+ = 0.1 and the vertical component F, for z+ = 0.1, and of (b) D (M = 734, Ng = 664, N, = 1.48 and m = O-500). however, need not be monotone in S.19 This only occurs if NJu~,~ > 1. Under this condition we find a unique value S,, the threshold saturation, for which Fz satisfies (2.19)
J’,(S,, at,z) = I;,(l, ut,,) = ar,z
size is constant in the Z-direction (3 nodes per meter). Convergence is obtained for the Picard iterations by adjusting the time steps. For the hydraulic head for each phase j at nodal point i,
and for S, < S < 1
F,(S, al+) > I;,(17 Q)
(2.20) we require that
One easily verifies that (2.19) is equivalent to k,,,(S,) = %
(2.21) g
The formulation of problem (2.17) shows that the twophase flow process is determined by the dimensionless numbers M, Ng, N,, and A, and by the exponent m in the nonlinearities.
3 NUMERICAL
SOLUTION
Using a numerical two-phase flow problem (2.11) on the finite domain -Hb < Z < Ht (Hb > HJ. Hence, conditions at the lower (no-flow) and pressure) boundaries
model we solve E < R < Rb and we impose also right (hydrostatic
Uj,r = 0
for E < R < Rb, Z = -Hb
Pj=Pjl+%(Hl-Z) c
forR=Rb,
-Hb < Z < Ht
(3.1) Numerical computations are made with nontransformed physical variables. The flow domain is discretized by linear triangular finite elements and the time discretization is fully implicit. The resulting algebraic equations are solved by the modified Picard method,4 giving good mass balance. In the R-direction, the grid is finest close to the Z-axis (19 nodes in total, the width of an element is 1.08 times the width of the previous element), whereas the grid
is reached within 20 Picard iterations. Here k refers to time and n to the Picard iteration level. Due to large differences in density and viscosity between air and water, the flow problem is dominated by advection. This requires relatively small time steps to reach convergence in 20 iterations, and hence large computation times. The initial time step is 0.03 s and the maximum allowable time step is 15 s. For all computations we use the relative convergence error e, = 0.001. The following soil and fluid parameters are fixed during all computations: Kabs = 5.30 x lo-” m2. r$ = 0.390; b0 = 1.77 x 10m5Pas; pa = l*24kgm-3r pL, = 1.30 x low3 Pas 9 p,,, = 1.00 x lo3 kgmp3; g= 9.8 m sp2. The effective porosity 4 = 0.390 reflects a soil porosity of O-400 and a residual water content of O*OlO. Parameters involving the boundary conditions are: E = 5.00 x 10e2 m; Rb = 3.65 m; HI = 1.50m; Hl - H, = 1*OO m; PWtE - 1a47 x lo4 Pa; P,,t = 0.0 Pa, whereas Hb is adjusted for each computation, such that the lower boundary does not affect the air flow. The van Genuchten parameters Q! and m, the air injection velocity Uin and the filter depth H are varied, which implies the variation of the exponent m and the dimensionless numbers Ng, N, and A. The mobility ratio M = 73.4 is kept constant. In Table 1 the data for the reference case (case 1) and for the other cases are
324
M. I. J. van Dijke et al. L
0.17
0.17
0.00
0.00
0.00
0.00
-0.17
-0.17
-0.17
-0.33
-0.33
-0.33
-0.50
-0.50
-0.67
-0.67
-0.67
-0.63
-0.63
-0.83
-1.00
-1.00
- 1.00
-1.17
-1.17
-1.17
-0.50 N
N
-1.33 0.01
(6)
0.21
0.41
- 1.33 1 0.01
0.61
r
0.21
r
(b) 0
0.61
0.41
0.81 0.17
‘s
0.00
\
-0.67
0.21
0.41
0.61
-
-0.17
-
-0.33
-
-0.50
-
-0.67
-
-0.83
-
-1.00
-
-1.17
- -1.33 I.81
r
Fig. 3. Air saturation
contourplots
at (a) t = 0.302, (b) t = 1.21 and (c) t = 151 (case 1).
0.61
0.61
-1.33
Multi-phase flow modeling of air sparging
325
0.10
_..-..-..
case 10 _____ ____._____.. 0.08
0.009
0.06
case 14 _______ ______-_________________ ________________________---_____________-
k
m
case 3
0.006
0.04
i
0.003 j
0.02
i 1
o.ooo’
4.0
0.0
8.0
12.0
16.0
20.0
24.0
0.00 0.01
0.17
0.09
0.25
0.33
0.41
0.49
t Fig. 4. Air volume stored in the flow domain (cases 1,3, 10 and
14). summarized. Observe that in the reference case, the total injection rate Q = 5.00 m3 h-t. All results are presented in terms of the dimensionless variables (2.9). Typical contourplots of the air saturations for the reference (case are shown in Fig. 3. The time dependence is illustrated by showing the increase of the volumle of air that is stored in the domain
V(t)Jh,S =2~
-hb
- 27r p
c
blS(r,z,t)drdz c” rS(r, z, 0) dr dz
(3.2)
J-h6 JE
where rb = Rb/H. In Fig. 4, v(t) iS shown for cases 1, 3, 10 and 14. For some causesoscillations occur, which are caused by physical instabilities. When injected air reaches the water table an amount of stored air leaves the flow domain through the unsaturated zone, which takes place much faster than the injection of air. Hence, the stored air volume decreases. This process can occur several times. Note, however, that the resolution
Fig. 5. Cross-section of the air sat;ation (case 1).
profile at z = -0.2
of the calculations can also affect the extent of oscillation. The plots in Fig. 4 show that a steady situation is reached before 1= 20 for all cases considered. This is equivalent to about 2 h in real time. To characterize the numerical results, we show in Fig. 5 for the reference case a cross-section of the saturation profile at z = -0.2 at steady state. For a time tl > 20, we determine the saturation S, = S(e, -0.2, tl) and the radial position rl for which S(r,, -0.2, tt) = f&. Both S, and rl are shown in Fig. 6 as a function of the dimensionless numbers, which are normalized with respect to the reference case, indicated by the subscript 1
NC,,
=-$Cl
A,
(3.3)
=$
From the numerical results we conclude that the air saturation is much less than one in the entire flow domain and that the water pressure is
Table 1. Parameters and dimensionless numbers wed in computations. Relative to the reference case (case l), m is vrwied for casea 2-5, NE for cases 6-9, N, for cases lo-13 and A for cases 14-17 Case 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Q (m-l)
Vi, ( 10m3m s-l)
H (m)
m
2.00 2.00 2.00 2.00 2.00 0.800 1.oo 1.50 3.00 2.50 1.67 1*oo 0.500 1.20 1.46 1.64 2.57
4.42 4.42 4.42 4.42 4.42 11.1 8.84 5.90 2.95 4.42 4.42 4.42 4.42 442 4.42 4.42 442
4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 7.50 6.17 5.50 3.50
0.667 0.400 0.500 0.750 0+800 0.667 0.667 0.667 0.667 0.667 0,667 0.667 0.667 0.667 0.667 0.667 0.667
% 6.63 6.63 6.63 6.63 6.63 2.68 3.32 4.97 9.95 6.63 6.63 6.63 6.63 6.63 6.63 6.63 6.63
NC (10-3)
A
0.738 0.738 0.738 0.738 0.738 0.738 0.738 0.738 0.738 0.590 0.885 1.48 2.95 0.738 0.738 0.738 0.738
15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 5-59 8.26 10.4 25.6
M. I. J. van Dijke et al.
326 0.165
0.36
0.32
0.24
0.065
0.00
0.50
1.00
(4
1.50 2.00
2.50
3.00
3.50
4.00
0.16 0.00
4.50
0.50
(b)
normalized numbers
1.00
1.50 2.00
2.50
3.00
3.50 4.00
4.50
normalized numbers
Fig. 6. Dependence of (a) the maximum saturation S, and (b) the radial position, rl, on the dimensionless numbers. For instance, the m-curve reflects values of the parameters for cases 2-5, see Table 1.
almost hydrostatic at large times. This allows us to treat the steady-state air flow as a single-phase process. 4 BOUNDS AND APPROTUMATIONS 4.1 Onedimensional flow
To explain the small values of the air saturation in the flow domain, we consider a relatively simple auxiliary problem involving an air saturation, which is above the saturation of problem (2.17). For this auxiliary problem we assume, as in Fig. 7, that air is injected through a horizontal disk at z = -h, with a uniform rate u,,, = 1. This leads to a much higher flow rate at z = -h, than in the case of injection through the vertical filter. At z = 0 we impose the natural outflow condition &S/& = 0. Consequently, the air saturation in this simplified case serves as an upper bound for the air saturation in the case of injection through a vertical filter. The injection through the horizontal disk leads to a one-dimensional analogue of problem (2.17), for -h, < z < 0 in which u~,~3 1. We therefore consider the problem (4.1)
F(S,)
S(r, 0 < St,
for z = -h,,
t>O
as-0 az-
for z = 0,
t>O
s=o
for -h,
t=O
F(S)-D(S)g=
1
where the boundary condition at z = -h,, with l), ensures the injection of air P(S) = I;,@, %,r = only. Note that F is non-monotone if Ng > 1, see
D(S,)$
< 1
for - h, 5 z _<0,
t _>0
(4.2)
To check the validity of this upper bound we also make a comparison with the numerical solution. To this end we compare the saturation just above the injection filter, S,,,, = S(E, -h,, t), with the threshold value S, for ut,= = 1 and for all values of the parameters in a relevant range. For example, Fig. 8 shows S, as a function of m and Nr, as S, appears to depend on these numbers, see (2.21).
t>O
$+;(F(S)-D(S)$)
= If
implying S,,, < S,. Hence, we have shown
in -h,
I
=0
(2.20). We show, below, that the saturation of the auxiliary problem (4.1) is below the threshold saturation S,, see (2.19). By an elementary argument based on the maximum principle,20 one finds that a solution of eqn (4.1) attains its maximum value at z = -h,, with &S/az < 0. This statement can also be obtained on physical grounds, when remembering that eqn (4.1) describes injection of air from the bottom into a medium which is completely filled with water at t = 0. To assess the maximum value, denoted by S,, we find from the boundary condition at z=-h u
(4.1)
The saturation just above the filter is chosen because it appears to be stationary and maximum on the flow domain.
Multi-phase Jlow modeling of air sparging
327
Consequently, we obtain in terms of the air saturation the approximate problem ( V - kr,(NgZZ - N,Vp,) = 0
II
_Nk
c
M
s-1 &
in r > c, z < h, for r = E, -hl < z < -h,
-
for r = c, z < -h,, -h,
PC =Pt
I
Pc-$Pt-
< z < ht
for r > E, z = ht $(ht c
-z)
for
r-+cq r>c,
z
(4.5)
I
Fig. 7. Schematic of the domain (01) for air injection through a horizontal disk.
In view of the large difference between the densities of air and water we assume that upwards from the injection titer and below the water table, flow in the z-direction is dominated by advection, i.e.
4.2 Steady state flow
(4.6)
At steady state the tilme derivative in the saturation equation in problem (2.17) vanishes. The reduced equation must be solveid subject to the same boundary conditions. Inherited from the initial pressure distribution we impose in addition pc + pt - 2
c
(h, - z) for ’ ---f O”’
I > E,
z < h,
z + -_oo
(4.3)
The numerical results ;show that the water pressure is approximately hydrostatic at large times. This implies that the water velocity can be disregarded with respect to the air velocity, i.e. Zt = ii=. Using this in eqn (2.13) yields ii,, = N,k,,e’, - N,k,,Vp,
(4.4)
As a consequence we are left with a nonlinear diffusion equation in the variables r and z. When solving this equation, no boundary condition can be imposed at z = ht, representing a ‘future time’. Hence, we consider the flow domain above the well of infinite extent, represented by z > -h,. Interpreting solutions only up to z = 0 will produce a reasonable approximation in this region. Assuming furthermore that the well-radius 6 < 1, we arrive at the problem, setting k(r, z) = k,,(S(r,z))
in r > 0, z > -h, $O,z)
=o
for z > -h,
t k(cqz)
= 0
for z > -h,
0.45I--
0.80
0.40 .
0.70
0.35 -
0.60
0.30 : m
d ,I' /' /In /' ,' /' /' hux
0.15 -
(4
0.50 v)
St
0.20 -
0.050.20
/IP
/
0.25
0.10 -
(4.7)
0.40 0.30
P-
/
/'
0.30
0.40
0.20 0.10
0.50
0.60
m
0.70
0.80
0.90
0.0
1.00 (b)
%lU 2.0
4.0
6.0
8.0
10.0
._
12.0
N,
Fig. 8. Comparison of S, for u,,, = 1 and S,, as a function of (a) the exponent m and (b) the gravity number Ng. The m-curve reflects values of the parameters for cases 2-5 and the NE-curve for cases 6-9 of Table 1.
328
M. I. J. van Dijke et al.
Appendix B(k) - nDkPD
(4.11)
where 2(1 - d,(m-1)/(4m+l), 4m+l
nD=
2(1 -m) “= 0.00 0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
Fig. 9.
Note PD <
k
4m+l
that
in
n,
>
()
(4.12) o
’
most
cases
of
practical
interest
l.
Hence, following Ref. 1 and 18 we find that for sufficiently small k the solution of eqn (4.7) with the conditions (4.9) and (4.10) can be approximated by eqn (4.13)
Graph of B(k) with m = 0.750.
(4.13)
where N1 = N,/N, and D(k) = kdp,/dk. In Fig. 9 we show an example of b(k). At z = --Jr, we need to impose an ‘initial condition’, which should reflect the influence of the injection filter. First, we observe, by a straightforward continuity argument that ccl ru,,,dr = A for all z 2 -h, 2lr (4.8) s0 Using this approximation rewritten in terms of k
00 J
rk(r, z) dr = F
implied by (4.6) this can be
for all 2 2 -h,
0
(4.9)
where N2 = A/TN,. For computational purposes we assume that the injection of air only takes place at the centre point of the filter, i.e. r = 0, z = - 1. This leads us to consider (4.7) and (4.9) for all z > -1, subject to the ‘initial condition’ k(r, -1) = 0
for all r > 0
(4.10)
In other words, the injection of air is described by a delta-distribution at the point r = 0, z = - 1. Only for special diffusion coefficients b(k) eqn (4.7) admits an explicit solution. For example when B is a power law function a solution is known in terms of similarity variables. This solution is known as the Barenblatt-Pattle point source solution.11’8 As a result of the upper bound derived in the previous section, we may conclude that k is small throughout the flow domain. We use this observation to approximate B(k) by a power law function, see the
where f(z) = fo . (z + l)‘@(PD + ‘1)
forz>-1
(4.14)
and fo = r$Ny)
1”2(pDi”‘(?$)j
(4.15)
Observe that the function f(z) represents the interface that separates the region where k, S > 0 from the region where k, S = 0. In the mathematical jargon the set {(f(z>,z); z 2 -1). IS k nown as the free boundary of eqn (4.7). Strictly speaking, the Barenblatt-Pattle solution (4.13) is valid for all z > - 1. However, since eqn (4.7) is formulated for z above the injection filter, i.e. z > -h,, we restrict the solution to these values of z. Furthermore, since (4.11) is a good approximation for small k, we expect that the quality of the approximation (4.13) increases with z. Using eqn (A4) the expression for k is easily converted into an expression for the saturation S.
5 APPLICABILITY OF THE ANALYTICAL STEADY-STATE APPROXIMATION The approximation for the steady state depends only on the three numbers Ni, N2 and m and is valid for z > -h,. We check whether the numerically obtained steady-state solution also depends on these three numbers for z > -h,. This is done by simulations 18 and 19, for which the input parameters are shown in Table 2. The other parameters are the same as in Section 3. Observe that although the numbers Nz, NC and A are different Ni and NZ are the same.
Multi-phase flow modeling of air sparging
329
Table 2. Parameters and dimensionless numbers for checking dependence of tbe steady-state solution on N1, Nz and m Case
ct (m--l)
Ui* (10e3 ms-‘)
H (m)
m
%
NC (10-3)
A
4.42 2.21
7.50 5.30
0.667 0.667
6.63 13.3
0.885 1.77
5.59 11.2
1.00 1.41.
18 19
Computations result in almost identical contourplots of air saturation at steady-state, i.e. at large time. However, the transient behaviour is different. This is illustrated in Fig. 10, where the stored air volume V(t) is plotted. We conclude that both steady states are identical. Therefore, also the numerical solutions at steady state depend only on m and the (combined) numbers Nt and N2 for -h, < z < 0. We investigated the agreement between the analytical and numerical solutions at steady state in terms of the relative permeability, k(r, z). For the numerical solution, k is calculated from the air saturation by the van Genuchten relation (2.5). For -h, < z < 0, the relative permeability contourplots of both the analytical and numerical solution for case 1 are shown in Fig. 11. Taking the cross-section of the solutions for z = -0.2, we use k,,(r) and knum(r) to represent the analytical and numerical solutions, respectively, as shown in Fig. 12. We find that the two profiles give excellent qualitative agreement. To obtain a quantitative criterion, we compute the masses involved, i.e.
z,, =
J
see (4.9), and
Note that Z,, depends. only on N1, N2 and m. To plot the relative error AZ = (I,, - Z,,,)/Z,, for various values of m, Ni, and Nz, we make use of Table 3, which reflects the same computations as shown in Table 1. AZ is shown in Fig. 13 as a function of the dimensionless numbers, which are normalized with respect to the reference case, indicated by the
0.010
case 19 __-------
T-
1 , -* I’
0.006 F=
,I’
0.004
0.000 0.0
2.0
4.0
6.0
8.0
10.0
t Fig.
10.
m
m,=--,
Nln
NI Nit
N2
N2n = -
=-9
ml
N21
for those cases in which only one of the numbers is varied. We observe that the relative error is nearly constant and small (less than 2%) for most cases, except for values of the parameters taken from case 17. In this case the injection depth, h, is small and as a consequence both the geometry of the filter and the outflow boundary condition have a relatively large influence on the numerical approximation. We may conclude that at steady state (T > 2 h) the complex two-phase flow process can be accurately approximated by the relatively simple expressions (4.13) and (4.14). An important qualitative feature of the solution is its radius of influence, r = f (0) = fo, see (4.15), which denotes the outflow radius of the injected air.
6 CONCLUSIONS
O” rkan(4 dr 0
0.008
subscript 1:
Comparison of stored air volumes for cases 18 and 19.
We modeled continuous air injection below the groundwater levels as multi-phase flow of two immiscible fluids, with emphasis on the development of the region containing air and the air saturation profile within this region. We found that the solutions depend on four dimensionless numbers, composed of soil, fluid and injection parameters. Using a numerical model based on the mixed form of the Richards equation for two phases, we made Table 3. Dimensionless numbers characterizing tbe steady state Case
m
N1
N2 (10-4)
1 2 3 4 5 6 7
0.667 0.400 0.500 0.750 0.800 0.667 0.667 0.667 O-667 0.667 0.667 0.667 0.667 0.667 0.667 0.667 0.667
8.99 8.99 8.99 8.99 8.99 3.59 4.49 6.74 13.5 11.2 7.49 4.49 2.25 8.99 8.99 8.99 8.99
7.45 7.45 7.45 7.45 7.45 18.6 14.9 10.0 4.96 7.45 7.45 7.45 7.45 2.68 3.97 4.98 12.3
: 10 11 12 13 14 15 16 17
M.
330 -O.OO
0.01 I
-0.11
N
0.11
0.21
0.31
0.41
I. J. van Dijke et al.
0.51 -O.OO
-0.00
0.01 r
t-
-0.44
-0.44
-0.56
-0.56
-0.56
-0.56
-0.67
-0.67
-0.87
-0.67
-0.76
-0.78
-0.78
-0.76
-0.89
-0.89
-0.89 0.11
0.01
L
-I 0.31
0.21
0.41
0.51
N
E
0.01
r
(4
-0.44
-0.89
computations for different values of the dimensionless numbers. The computations provided accurate quantitative information about the development of an air cone in the flow domain and about the corresponding saturations. For some cases, the numerical computations showed instabilities during the development of the air cone. Within a short time (about 2 h) a steady-state situation was reached, which was approximately a single-phase (air) flow process. For all cases air
0.41
0.51
r
(b)
Fig. 11. Comparison of the relative permeability contourplots of
0.31
0.21
0.11
the analytical and (b) the numerical solution (case 1).
saturations were small in the entire flow domain. Unfortunately, the numerical computations were very time consuming. An alternative for the use of the two Richards equations is the fractional flow formulation. This formulation involves one equation for the air saturation and one time-independent equation for the total velocity. It presents an understanding of the
O.1° I 0.024
,
0.016
0.008 -0.02
’ 0.20
J 0.40
0.60
0.07
0.14
0.21
0.28
0.35
0.42
r
Fig. 12. Cross-section of the relative permeability profiles k,,(r) and knm(r) at z = -0.2 (case 1).
1.00
1.20
1.40
1.60
1.80
normalized number8
0.000 0.00
0.80
AI of the volume integrals versus the dimensionless numbers. The m-cuNe reflects values of the parameters for cases 2-5, the Nl-curve for cases lo-13 and the N2-curve for cases 14-17 of Table 3. Fig. 13. Relative errors
Multi-phase flow modeling of air sparging nonlinearities of the equations and is used for derivation of an upper bound on the air saturation in the flow domain and an approximate analytical solution for the steady state. The observation of small air saturations enabled us to approximate the complicated expressions for the relative permeability and capillary pressure functions by simple power law functions. Assuming that flow in the vertical direction is dominated by advection, we derived the analytical solution for the steady-state situation which is valid for the flow domain above the injection filter. This analytical solution was proven to be in good agreement with the numerical results and it depends on three dimensionless numbers, which are a combination of the four numbers mentioned before. For practical reference we give the expression of the approximate analytical solution for the air saturation and the expression for the dimensionless numbers involved in this approximation. The numbers are N1
=
(Pw- PaWa
(6.1)
Pw N
=
2
&aUinE(HI- Hu) &bshv
-
(6.4
dgH2
and the exponent m in the nonlinearities. The analytical solution is given by
I
ns(“@$$))-(1
S(r,z) =
IO
-h,LzlO
for r > f(z),
-h,
where f(z) =fo
* (Z +
Table 4. Oufflow radii
Case 1 2 : :
FO Cm)
Case
2.31
7 8 9 10 11 12
1.61 1.87 2.53 2.65 3.65
FO @I
Case
FO Cm)
13 14 15 16 17 18
4.15 3.56 3.02 2.74 1.87 3.85
3.27 2.67 1.89 2.10 250 3.10
Observe that the coefficient ps/po = l/( 1 - m) in eqn (6.3) satisfies ps/pD > 1 for all 0 < m < 1, which means that the saturation profile in the r-direction always has zero derivative at the free boundary. Hence, the profile shows a relatively large zone where S is almost zero. Since for remediation purposes a minimum air saturation is necessary, we need to know of the air saturation profile to determine at what distance from the injection axis this minimum saturation is reached. In our calculations we imposed at the filter an influx, which is equally distributed over the length of the filter. However, in practical situations most air leaves the filter near the top.15 Because for the approximate analytical solution we assumed point source injection, the use of this approximation can be improved by setting the injection level higher than the centre of the filter.
ACKNOWLEDGEMENT
- (&)2r’pD
forOIrIf(z),
331
1)1’(2(pD+1)) for - h, 5 z 5 0
This project was sponsored by the Netherlands Organization for Scientific Research (NW0 project NLS 61251) and carried out in the framework of the Netherlands Integrated Soil Research Programme (PCTB project 5029), with the use of super computing facilities of the Dutch National Computing Facilities Foundation (NCF).
(6.4) with REFERENCES
represents the outer botmdary of the air cone. Note, that in eqn (6.3) dimensionless spatial coordinates r and z need to be inserted. To obtain FO, the radius of influence in a real physical dimension, f0 must be multiplied by the characteristic length H. Furthermore, 2 ns = m4ml(4m+l), Ps=4m+ 1
31 -m)m(m-k)j(4m+l),pD =W -4 4m+l Q= 4m+l For reasonable values of m, such as 0.30 < m < 0.90, we have O-52 < ns < O-92, 0.43 < ps < @90, O-043 < nD < 0.93 and 0.043 < po < 064. To give an impression of the width of air cones, we list the outflow radii F0 in Table 4 for the computations of Section 3.
1. Barenblatt, G. I., On some unsteady motions of a liquid or a gas in a porous medium. Prikl. Mat. Mekh., 16 (1952) 67-78 (in Russian). 2. Bear, J., Dynamics of Fluids in Porous Media. Elsevier, New York, 1972. 3. Bohler, U., Brauns, J., H&l, H. & Nahold, M., Air injection and soil-air extraction as a combined method for cleaning contaminated sites - observations from test sites in sediments and solid rocks. In Contaminated Soil ‘90, ed. F. Ahrendt. Kluwer Academic Publishers, Dordrecht, 1990, pp. 1039-44. 4. Celia, M. A., Bouloutas, E. T. & Zarba, R. L., A general mass-conservative numerical solution for the unsaturated flow equation. Water Resow. Res., 26 (1990) 1483-96. 5. Chavent, G. & Jaffre, J., Mathematical Models and Finite Elements for Reservoir Simulation. ed. J. L. Lions, Studies
in Mathematics and its Applications, Vol. 17. NorthHolland, Amsterdam, 1986. 6. Crow, W. L., Anderson, E. P. & Mint@, E. M., Subsurface
venting
of vapors
emanating
from hydro-
M. I. J. van Dijke et al.
332 carbon
product
on ground
water. Ground Water Mon.
Rev., 7 (1987) 51-7. 7. Forsyth, P. A. & Shao, B. Y., Numerical simulation of gas venting for NAPL site remediation. Adv. Water Resour., 14 (1991) 354-67. 8. Gilding, B. H., Qualitative mathematical analysis of the Richards equation. Transport in Porous Media, 5 (1991) 651-66. 9. Herrling, B. L Stamm, J., Numerical results of calculated
3D vertical circulation flows around wells with two screen sections for in situ or on-site aquifer remediation. In
nSkPs and nokPD, respectively,
where nj andpj (j = S, D) are constants depending on m. The relations (2.5) and (2.12) for the capillary pressure and the relative permeability are:
PC(S) = ((1 - S)+
k(S) = Si(l - (1 - S)1’“‘)2m Expanding
Computational Methods in Water Resources IX, Vol. 1: Numerical Methods in Water Resources, ed. T. F. Russel.
Elsevier Applied Science, London, 1992, pp. 483-92. 10. Ji. W., Dahmani. A.. Ahlfeld. D. P.. Lin. J. D. & Hill III. E: HI, Laboratory. study of air’ sparging: Air flow visualization. Ground Water Mon. Rev., 13 (1993) 115-26. 11, Johnson, P. C., Kemblowski, M. W. & Colthart, J. D., Quantitative analysis for the cleanup of hydrocarboncontaminated soils by in-situ soil venting. Ground Water, Johnson, P. C., McWorther, D. B., Goodman, I., An overview of in-situ air Water Mon. Rev., 13 (1993) 127-35.
J. & Parker, J. C., An efficient finite for modeling multiphase flow. Water
(1 - S)“”
(l@=l-f-S+;
in a binomial k-1
m
(
k(S) = m -2y+m+l+ = @k(S))
For &(S) Q(S)
(1995) l-27.
15. Marley, M. C., Hazebrouch, D. J. & Walsh, M. T., The application of in-situ air sparging as an innovative soils and ground water remediation technology. Ground Water Mon. Rev., 12 (1992) 137-45. 16. Mohtar, R. H., Wallace, R. B. & Segerlind, L. J., Finite element simulation of oil spill cleanup using air sparging. In Computational Methodr in Water Resources X, Vol. 2, ed. A. Peters. Kluwer Academic Publishers, Dordrecht, 1994, pp. 967-74. 17. Parker, J. C., Lenhard, R. J. & Kuppusamy, T., A parametric model for constitutive properties governing multiphase flow in porous media. Water Resour. Res., 23 (1987) 618-24.
18. Pattle, R. E., Diffusion from an instananeous point source with a concentration-dependent coefficient. Quart. J.
;S2+0(S3) >
O(S4”‘1)
W)
= k(dp,/dS)/(clk/dS)
we find
2(1- m)
=
m
Resour. Res., 25 (1989) 43-54.
14. Lundegard, P. D. & LaBrecque, D., Air sparging in a sandy aquifer (Florence, Oregon, USA): Actual and apparent radius of influence. J. Contam. Hydrol., 19
series
for k
we obtain
28 (1990) 413-29.
12. Johnson, R. L., Hinchee, R. E. & sparging. Ground 13. Kaluarachchi, J. element method
- l)‘-m
S(l - (1 - s>+r-” x (1 -S)““(l Again
- (1 -S)‘/m+4S(l
using the series (Al),
&(S)
=
-,)++I)
we obtain
2~~~~)mm-1S1pm + O(S2-2m)
w
which is only accurate for S close to zero, Taking first order terms we get
S(k) N m 4”l/(4m+
1) k2/(4t?Z+
1) (A41
Observe that these approximations are only good for small m, because when m r 1 the first and the higher
Mech. Appl. Math., 12 (1959) 407-9.
19. Proskurowski, W., A note on solving the Buckley-Leveret equation in the presence of gravity. J. Comput. Phys., 41
0.24 ,
(1981) 136-41.
20. Protter, M. H. & Weinberger, H. F., Maximum Principles in Dt&ential Equations. Prentice-Hall, London, 1967. 21. Smoller, J., Shock Waves and Reaction D@usion Equations, Springer Verlag, New York, 1983. 22. Wehrle, K., In-situ cleaning of CHC contaminated sites: Model-scale experiments using the air injection (in-situ stripping) method in granular soils. In Contaminated Soil ‘90, ed. F. Ahrendt. Kluwer Academic Publishers, Dordrecht, 1990, pp. 1061-2.
0.20 -
0.16 ________-----
________-----n
0.12 ____-_-__-
___---
’
0.00 0.00
APPENDIX According to the observation that for z > -h,, S(r, z) is significantly smaller than one, we approximate S(k) and b(k) = kdp,/dk by power law functions of the form
0.03
0.06
0.09
0.12
regression
0.15
0.18
0.21
k
Fig. Al. Approximations of b(k) with m = 0.750 for k < 0.151. Approximation by binomial series yields (nn,pn) = (0*127,0.125), whereas regression analysis yields (nn,pn) = (0.182,0:146).
Multi-phaseflow modeling of air sparging order terms of eqns (AZ!) and (A3) become of the same order of magnitude. We get ns = m h/(4m+ 2 “=4m+l
1)
7
2(1 - “),(m-:l)/(4m+l), nD= 4m+l “=
2(1 -m) 4m+l ’
O 0
0
WI
333
Note, that using the upper bound for S, i.e. S < S,, see (4.2), or equivalently k < l/N,, it is possible to find other power law functions by regression analysis, which approximate D(k) for 0 5 k 5 l/N,. Figure Al shows the function D(k) with approximations for m = 0.750, which must be valid for k 5 1/IVr = 0.15 1. Because the approximation by binomial series shows the correct limiting behaviour of b(k) for k JO, we use this approximation for the derivative of the analytical solution.
0309-1708(95)00028-3
ELSEVIER
Multi-phase flow modeling of air sparging M. I. J. van Dijke, S. E. A. T. M. van der Zee Department of Soil Science and Plant Nutrition, Wageningen Agricultural University, PO Box 8005, 6700 EC Wageningen, The Netherlantis &
C. J. van Duijn Department of Mathematics, De@ University of Technology, PO Box 5031, 2600 GA De@, The Netherlands (Received 14 February
1995; accepted 25 July 1995)
Air injection into groundwater (air sparging) in a homogeneous axially symmetric porous medium is modeled using a two-phase flow approach. A numerical method based on the mixed form of the Richards equation for both phases is presented. Furthermore two analytical approximations are discussed to explain the numerical results. One is a one-dimensional description explaining the occurrence of small air saturations. The other is a closed form approximation for the distribution of the air saturation in the resulting steady state. From the latter we can estimate the maximum radius of influence of air sparging, as a function of the physical parameters. The analytical approximation at steady state and the numerical results are in good agreement. Key words: groundwater, air sparging, multi-phase influence, analytical solution.
NOTATION
Hu (4,)
Dimensionless variables between brackets ( ). . Subscript denoting air. Dimensionless filter surface. D Dimensionless diffusion function for two-phase flow dependent on S. Dimensionless diffusion function for b single-phase flow dependent on k. Unit vector lin z-direction. Dimensionless position of free boundary. z-Component of flux function. Flux function. (Dimensionless) position of free boundary at water tab1.e [ml. Air fractional flow function. h Gravity [m s2]. g H Filter depth [ml. (Dimensionless) depth of bottom boundary Hb (hd
k K abs
Row, steady state, radius of
(Dimensionless) depth of upper side of filter
[ml.
1
HI (hJ
k]* rmensionless) depth of lower side of filter
Ht (ht)
E]* rmensionless) height of top boundary
krj A4
:l,
N2
N, Njg Ng pc Pjt
(PA
(Pjt>
Pt
Q R (4 & kb) S
si Sr St
[ml. 319
Redefined air relative permeability. Absolute permeability [m’]. Fluid j relative permeability. Mobility ratio. van Genuchten parameter. Dimensionless numbers determining steady-state situation. Capillary number. Fluid j gravity number Gravity number. (Dimensionless) capillary pressure [Pa]. (Dimensionless) fluid j pressure at top boundary pa]. Dimensionless capillary pressure at top boundary. Air injection rate [m3 s-t]. (Dimensionless) horizontal coordinate [ml. (Dimensionless) width of right boundary [ml. Redefined air saturation. Fluid i saturation. Residual water saturation. Air threshold saturation.
320 T, (0 uj <$> uin 4 V l (4
i(f) x & Pj 4
M. I. J. van Dijke et al.
(Dimensionless) time [s]. (Dimensionless) fluid j Darcy velocity [m s-l]. Air injection velocity [m s-l]. Dimensionless total fluid velocity. Dimensionless air volume in flow domain. Subscript denoting water. (Dimensionless) vertical coordinate [ml. van Genuchten parameter. (Dimensionless) filter radius [m]. Mobility function. Fluid j viscosity [Pa s]. Fluid j density [kg me3]. Porosity.
1 INTRODUCTION A method for remediating an aquifer, which is contaminated by organic liquids (solvents, gasoline) trapped in the saturated zone, is to inject air or oxygen into the aquifer. Injection of air may enhance microbial degradation and volatilization. For remediation of the unsaturated zone, gas venting has been studied both experimentally and numerically.637’11Besides the use of soil vapor extraction wells, injection of air at a small distance below the water table to remove contaminants at groundwater level has been considered.7 For remediation of the saturated zone the use of a vacuum vaporizer well has been described.g Emphasis was laid on assessing the sphere of influence of the water circulation, where dissolved oxygen enhances microbial degradation. Direct injection of air in the saturated zone, known as air sparging, together with vapor extraction in the unsaturated zone, has been put forward as an effective in-situ remediation technique. Air sparging has been studied experimentally on both field3>i5and laboratory’0’22 scales in the last few years, with emphasis on the region in the saturated zone where air is present (radius of influence). A numerical study of air sparging at steady state has been provided for the case of injection just below the original water table.16 Although laboratory experiments with air injection in a glass bead medium” showed that air flow can occur as isolated bubbles for bead diameters exceeding 2 mm, it is assumed that under natural subsurface conditions air flow is most likely to occur in small continuous channels.‘0112 A recent field study for a uniform soil with mean grain size of 0.25mm,i4 indicated that the density of channels in the main region of air flow must be very high. Therefore, it is reasonable to model air flow macroscopically as a continuum. Air sparging is only possible in relatively coarsegrained soils. The reported injection rates varied from about 3 m3 h-’ up to 100m3 h-’ for thick sandygravelly deposits.3112 Injection in fine textured soils
requires high air entry pressures, which may cause fracturing of the soil. This may result in the formation of a few channels through which air flows upwards.” Hence, knowledge of soil layering and heterogeneity is important because these affect the radius of influence.12 Besides flow continuity, air-phase compressibility may be a complicating factor in modeling air sparging. In the early stage of sparging when air pathways from the injection filter to the unsaturated zone have not yet been established, air density is not constant. However, emphasizing the steady-state situation in which continuous channels to the vadose zone exist, compressibility is expected to play a minor role. This expectation is supported by the field study,14 which describes sparging at a rate of 34.2 m3 h-l, resulting in a pressure increase of 10.5 kPa at a distance of 0.6m from the filter. This causes an air density increase of about 10% only. Important for dimensioning the technique of air sparging is quantitative knowledge of the effect of soil, fluid and filter parameters on the radius of influence of an injection filter. This motivated us to study an axially symmetric model for air injection through a vertically positioned filter in an initially saturated region below the vadose zone. In the model we consider air and water as two immiscible incompressible continuous phases. The interaction between the fluid phases and the soil matrix is described by the saturation dependent relative permeability and capillary pressure functions. Contaminants are assumed to be part of the soil matrix and do not affect the flow process. The basic equations are given and discussed in Section 2, where they are reformulated in terms of dimensionless numbers (that are combinations of the physical parameters). In Section 3 we discuss a numerical method, which is based on the mixed form of the Richards equation for both water and air. Results of numerical computations are presented in terms of the distribution of air saturations and the volume of air that is stored in the domain. In Section 4 we discuss two analytical approximations which are valid in the relevant part of the flow domain. One is an upper bound that explains the occurrence of small saturations. The other is an explicit solution for the steady-state situation, which is derived under the assumptions that water is immobile and that gravity is the dominant effect in the vertical direction. The explicit solution provides an expression for the radius of influence of sparging. In Section 5 we quantify the agreement between the analytical and numerical solutions at steady state.
2 MODEL We use Darcy’s Law for both air (a) and water (w) q = _5!E+(Pj J
+pjg
j= w,a
(2.1)
Multi-phase flow modeling of air sparging
and the mass balance e:quations aSj + $,t+V’Uj=O
the water table is situated at Z = 0. Hence,
j=w,a
(2.2)
In eqn (2.2) Sj denotes the effective fluid saturation and C#J the effective porosity, i.e. redefined according to s
SW-sr
=
w:
&
Sll := 1 _ $7
1-S’
C$:= (1 - S,)C$
where S, is the residual water saturation. Furthermore, Kabs is soil absolute permeability, CJ~soil porosity, Cj fluid Darcy velocity, Sj fluid saturation, krj fluid relative permeability, /Lj fluid viscosity, Pj fluid pressure, pj fluid density and g gravity. We assume that the soil is homogeneous and isotropic and that both fluids are incompressible. The set of equations., (2.1) and (2.2), is completed by the constitutive relati’ons S,,, + S, = 1, the capillary pressure PC = P, - P,,,, Sj = Sj(Pc) and krj = _k,i(Sj) and is solved for the unknowns S,,,, S,, P,, P,, U, and CO. The dependence of capillary pressure and relative permeability on the saturations is given by the wellknown expressions17 p&y,)
= +,lirn
- ,)I-,
Pj = Pj,
for j = w, a, where Pwr = -p,gHt and Pat = 0. Air is injected with velocity Vi, into this domain through a filter with radius E > 0 and length HI - H,, giving the t_otal injection rate Q = 2rE(Hl - HJUi,. Writing uj = (“j,r, V,,z), we have the flux boundary conditions UUJ,,= 0 uo,r =
(24
k&J
= Sf,(l - (1 - S0)1’m)2m
(2.5)
Equations (2.1) and (2.2) are solved in the twodimensional axially symmetric domain of Fig. 1. The level Z = 0 corresponds to the initial position of the water table, where P, := 0. Along the top boundary of the domain, Z = H,, the air pressure equals P, = 0. The water pressure along this boundary is selected such that
=
Pj = Pjt \
pi = pjr + 9
c
(h, .-
Z)
-c Z < Ht
Pj = Pjt + pjg(Hf - Z) for r > E, Z < Ht
(2.8)
We introduce the dimensionless variables
TUin
_ 3j uj=Y,y
R
r=-
t=-@,
H’
z=-
Z H
(2.9)
0Pj Pj = Pg w
j=w,a
h,,+‘,
/q+
h$,
pit=-- ffpjf P&T
The flow problem is now determined by the mobility ratio, the gravity numbers for water and air, the capillary number and the dimensionless filter surface, i.e.
&&CL,,
N,
Pa N
C
=
18
KabsPwg
=
~a Ui,Ha ’
KabsPjg Pa uin
(2.10)
Q A=Ui,H*
For both water and air, eqns (2.1) and (2.2) are combined into the mixed form of the Richards equation.41’3 Together with the boundary and initial conditions the resulting problem is given in (2.11)
t>O
+ N,,=))
for r = E, z < -hl,
0
-H,
The initial conditions at T = 0, say, are
for r = E, -h, < z < -h,, Uj,l
forR = E, Z < -H,,
(2.7)
in r > E, z < h,,
<
for R = E, -HI < Z < -H,, 1
Uj,, = 0
13,17 C’
= V - (k,V(Ncpa
uin
E E=---, H
= St,(l - (1 - Sy”)“)*
at
(24
(2.3)
k&Y,,,)
-
for R > E, Z = Ht
where H = (H, + Hl)/2, the depth of the filter center, is chosen as a characteristic length. Furthermore, we define the dimensionless constants
where 0 < m < 1 and (Y> 0 are given constants. Note that S,,,(P,) is dehned for PC 2 0. For P, < 0 we set S, = S,(O) = 1. Hence, 0 5 S, 5 1 for all values of p
321
-h,
t>O
< z < h,,
t> 0
for r > E, z = h,,
t>O
for r > E, z < h,,
t=O
(2.11)
M. I. J. van Dijke et al.
322
where
Mkr=
f, =
(2.14)
kw + Mkra denotes the air fractional flow function and
(2.15) denotes the mobility function. Ng = Nwg - Nag and e’, is the unit vector in the z-direction. Further, i& = u’, + il, denotes the total velocity, which satisfies the incompressibility condition v.u’,=o
Observe that in case gravitational and capillary effects are neglected, the air velocity ii, is the fraction f,of the total velocity. Setting S = S, we arrive at the initial-boundary value problem (2.17)
Fig. 1. Schematic of the domain (s2)for air sparging.
dt
+
V
(faz-it + N&”
.
f, - NJ%
= 1
- NJVp,)
(2.16)
= 0
for r > E, z < h,,
t>O
for r = E, -h, < z < -h,,
t>O
(2.17)
24t,r -1 -
ds 0 ar= U
for r = c, z < -h,,
-h,
< z < h,,
t>O
1,r -0 -
PC =Pt , Pc=P+-
N (h, - ‘) c
for r > C, z = h,,
t>O
for r > E, z < h,,
t=O
where pwr = -(N,,/N,)h, and par = 0. The dimensionless capillary pressure, which depends now on S,, is, according to the redefinition of the phase pressures pc(Sa) = ((1 - s,)-“m
- 1)1-m
(2.12)
with S, = 0 for pc < 0. Problem (2.11) is solved for one of the phase saturations and one of the phase pressures. In general, however, it cannot be solved explicitly due to the nonlinear nature of the equations. To develop an understanding of the behaviour of the solution, we consider a numerical solution technique based on the formulation of problem (2.11). In addition, we characterize some aspects of the solutions analytically. To achieve this we construct a single equation for the air saturation only. We find in terms of Sa:2T5
ii, = f,i& + N,Xe’, - NJV,,
(2.13)
where pr =pot -pwt = (N,,/N,)h, = cd!Zh,. Note that this problem is equivalent to the original description (2.11) in terms of the separate phase saturations and pressures. However, to understand the qualitative behaviour of the solution, one often considers iit in (2.17) as given, leading to a nonlinear advectiondiffusion equation in terms of S only, which is more accessible to analysis. In problem (2.17), f,(S), X(S) and p,(S) are nonlinear functions of S. These nonlinearities are grouped in the functions &S,
i&) = f,(S)&
+ N,X(S)e’,
(2.18) D(S)
= N,X(S)
g
(S)
In Fig. 2 we show some typical examples of (2.18). We note here that the r-component of the flux, F,, is increasing in S for all z.+ > 0. The z-component, F,,
Multi-phase flow modeling of air sparging 0.14 0.12 O.1o
323
0.030 ,_----__ F, --__ --__ I’ ___;; . . . . . . . . . . . . . . . . . . .
0.08
/: , : : :
0.06
:
j
!A
0.020
n
:
:
: :
: :
0.02
:’
~
0.00
--------_____
Fr
0.04
0.00 EIl”’
0.025
---_____
0.015
0.010
0.005
jst -I 0.06
0.12
0.18
0.24
0.30
0.36
0.42
S
(a)
0300
0.48
L 0.00
0.10
04
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
1.00
S
Fig. 2. Graphs of (a) the radial component F, of F’for z+ = 0.1 and the vertical component F, for z+ = 0.1, and of (b) D (M = 734, Ng = 664, N, = 1.48 and m = O-500). however, need not be monotone in S.19 This only occurs if NJu~,~ > 1. Under this condition we find a unique value S,, the threshold saturation, for which Fz satisfies (2.19)
J’,(S,, at,z) = I;,(l, ut,,) = ar,z
size is constant in the Z-direction (3 nodes per meter). Convergence is obtained for the Picard iterations by adjusting the time steps. For the hydraulic head for each phase j at nodal point i,
and for S, < S < 1
F,(S, al+) > I;,(17 Q)
(2.20) we require that
One easily verifies that (2.19) is equivalent to k,,,(S,) = %
(2.21) g
The formulation of problem (2.17) shows that the twophase flow process is determined by the dimensionless numbers M, Ng, N,, and A, and by the exponent m in the nonlinearities.
3 NUMERICAL
SOLUTION
Using a numerical two-phase flow problem (2.11) on the finite domain -Hb < Z < Ht (Hb > HJ. Hence, conditions at the lower (no-flow) and pressure) boundaries
model we solve E < R < Rb and we impose also right (hydrostatic
Uj,r = 0
for E < R < Rb, Z = -Hb
Pj=Pjl+%(Hl-Z) c
forR=Rb,
-Hb < Z < Ht
(3.1) Numerical computations are made with nontransformed physical variables. The flow domain is discretized by linear triangular finite elements and the time discretization is fully implicit. The resulting algebraic equations are solved by the modified Picard method,4 giving good mass balance. In the R-direction, the grid is finest close to the Z-axis (19 nodes in total, the width of an element is 1.08 times the width of the previous element), whereas the grid
is reached within 20 Picard iterations. Here k refers to time and n to the Picard iteration level. Due to large differences in density and viscosity between air and water, the flow problem is dominated by advection. This requires relatively small time steps to reach convergence in 20 iterations, and hence large computation times. The initial time step is 0.03 s and the maximum allowable time step is 15 s. For all computations we use the relative convergence error e, = 0.001. The following soil and fluid parameters are fixed during all computations: Kabs = 5.30 x lo-” m2. r$ = 0.390; b0 = 1.77 x 10m5Pas; pa = l*24kgm-3r pL, = 1.30 x low3 Pas 9 p,,, = 1.00 x lo3 kgmp3; g= 9.8 m sp2. The effective porosity 4 = 0.390 reflects a soil porosity of O-400 and a residual water content of O*OlO. Parameters involving the boundary conditions are: E = 5.00 x 10e2 m; Rb = 3.65 m; HI = 1.50m; Hl - H, = 1*OO m; PWtE - 1a47 x lo4 Pa; P,,t = 0.0 Pa, whereas Hb is adjusted for each computation, such that the lower boundary does not affect the air flow. The van Genuchten parameters Q! and m, the air injection velocity Uin and the filter depth H are varied, which implies the variation of the exponent m and the dimensionless numbers Ng, N, and A. The mobility ratio M = 73.4 is kept constant. In Table 1 the data for the reference case (case 1) and for the other cases are
324
M. I. J. van Dijke et al. L
0.17
0.17
0.00
0.00
0.00
0.00
-0.17
-0.17
-0.17
-0.33
-0.33
-0.33
-0.50
-0.50
-0.67
-0.67
-0.67
-0.63
-0.63
-0.83
-1.00
-1.00
- 1.00
-1.17
-1.17
-1.17
-0.50 N
N
-1.33 0.01
(6)
0.21
0.41
- 1.33 1 0.01
0.61
r
0.21
r
(b) 0
0.61
0.41
0.81 0.17
‘s
0.00
\
-0.67
0.21
0.41
0.61
-
-0.17
-
-0.33
-
-0.50
-
-0.67
-
-0.83
-
-1.00
-
-1.17
- -1.33 I.81
r
Fig. 3. Air saturation
contourplots
at (a) t = 0.302, (b) t = 1.21 and (c) t = 151 (case 1).
0.61
0.61
-1.33
Multi-phase flow modeling of air sparging
325
0.10
_..-..-..
case 10 _____ ____._____.. 0.08
0.009
0.06
case 14 _______ ______-_________________ ________________________---_____________-
k
m
case 3
0.006
0.04
i
0.003 j
0.02
i 1
o.ooo’
4.0
0.0
8.0
12.0
16.0
20.0
24.0
0.00 0.01
0.17
0.09
0.25
0.33
0.41
0.49
t Fig. 4. Air volume stored in the flow domain (cases 1,3, 10 and
14). summarized. Observe that in the reference case, the total injection rate Q = 5.00 m3 h-t. All results are presented in terms of the dimensionless variables (2.9). Typical contourplots of the air saturations for the reference (case are shown in Fig. 3. The time dependence is illustrated by showing the increase of the volumle of air that is stored in the domain
V(t)Jh,S =2~
-hb
- 27r p
c
blS(r,z,t)drdz c” rS(r, z, 0) dr dz
(3.2)
J-h6 JE
where rb = Rb/H. In Fig. 4, v(t) iS shown for cases 1, 3, 10 and 14. For some causesoscillations occur, which are caused by physical instabilities. When injected air reaches the water table an amount of stored air leaves the flow domain through the unsaturated zone, which takes place much faster than the injection of air. Hence, the stored air volume decreases. This process can occur several times. Note, however, that the resolution
Fig. 5. Cross-section of the air sat;ation (case 1).
profile at z = -0.2
of the calculations can also affect the extent of oscillation. The plots in Fig. 4 show that a steady situation is reached before 1= 20 for all cases considered. This is equivalent to about 2 h in real time. To characterize the numerical results, we show in Fig. 5 for the reference case a cross-section of the saturation profile at z = -0.2 at steady state. For a time tl > 20, we determine the saturation S, = S(e, -0.2, tl) and the radial position rl for which S(r,, -0.2, tt) = f&. Both S, and rl are shown in Fig. 6 as a function of the dimensionless numbers, which are normalized with respect to the reference case, indicated by the subscript 1
NC,,
=-$Cl
A,
(3.3)
=$
From the numerical results we conclude that the air saturation is much less than one in the entire flow domain and that the water pressure is
Table 1. Parameters and dimensionless numbers wed in computations. Relative to the reference case (case l), m is vrwied for casea 2-5, NE for cases 6-9, N, for cases lo-13 and A for cases 14-17 Case 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Q (m-l)
Vi, ( 10m3m s-l)
H (m)
m
2.00 2.00 2.00 2.00 2.00 0.800 1.oo 1.50 3.00 2.50 1.67 1*oo 0.500 1.20 1.46 1.64 2.57
4.42 4.42 4.42 4.42 4.42 11.1 8.84 5.90 2.95 4.42 4.42 4.42 4.42 442 4.42 4.42 442
4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 4.50 7.50 6.17 5.50 3.50
0.667 0.400 0.500 0.750 0+800 0.667 0.667 0.667 0.667 0.667 0,667 0.667 0.667 0.667 0.667 0.667 0.667
% 6.63 6.63 6.63 6.63 6.63 2.68 3.32 4.97 9.95 6.63 6.63 6.63 6.63 6.63 6.63 6.63 6.63
NC (10-3)
A
0.738 0.738 0.738 0.738 0.738 0.738 0.738 0.738 0.738 0.590 0.885 1.48 2.95 0.738 0.738 0.738 0.738
15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 15.5 5-59 8.26 10.4 25.6
M. I. J. van Dijke et al.
326 0.165
0.36
0.32
0.24
0.065
0.00
0.50
1.00
(4
1.50 2.00
2.50
3.00
3.50
4.00
0.16 0.00
4.50
0.50
(b)
normalized numbers
1.00
1.50 2.00
2.50
3.00
3.50 4.00
4.50
normalized numbers
Fig. 6. Dependence of (a) the maximum saturation S, and (b) the radial position, rl, on the dimensionless numbers. For instance, the m-curve reflects values of the parameters for cases 2-5, see Table 1.
almost hydrostatic at large times. This allows us to treat the steady-state air flow as a single-phase process. 4 BOUNDS AND APPROTUMATIONS 4.1 Onedimensional flow
To explain the small values of the air saturation in the flow domain, we consider a relatively simple auxiliary problem involving an air saturation, which is above the saturation of problem (2.17). For this auxiliary problem we assume, as in Fig. 7, that air is injected through a horizontal disk at z = -h, with a uniform rate u,,, = 1. This leads to a much higher flow rate at z = -h, than in the case of injection through the vertical filter. At z = 0 we impose the natural outflow condition &S/& = 0. Consequently, the air saturation in this simplified case serves as an upper bound for the air saturation in the case of injection through a vertical filter. The injection through the horizontal disk leads to a one-dimensional analogue of problem (2.17), for -h, < z < 0 in which u~,~3 1. We therefore consider the problem (4.1)
F(S,)
S(r, 0 < St,
for z = -h,,
t>O
as-0 az-
for z = 0,
t>O
s=o
for -h,
t=O
F(S)-D(S)g=
1
where the boundary condition at z = -h,, with l), ensures the injection of air P(S) = I;,@, %,r = only. Note that F is non-monotone if Ng > 1, see
D(S,)$
< 1
for - h, 5 z _<0,
t _>0
(4.2)
To check the validity of this upper bound we also make a comparison with the numerical solution. To this end we compare the saturation just above the injection filter, S,,,, = S(E, -h,, t), with the threshold value S, for ut,= = 1 and for all values of the parameters in a relevant range. For example, Fig. 8 shows S, as a function of m and Nr, as S, appears to depend on these numbers, see (2.21).
t>O
$+;(F(S)-D(S)$)
= If
implying S,,, < S,. Hence, we have shown
in -h,
I
=0
(2.20). We show, below, that the saturation of the auxiliary problem (4.1) is below the threshold saturation S,, see (2.19). By an elementary argument based on the maximum principle,20 one finds that a solution of eqn (4.1) attains its maximum value at z = -h,, with &S/az < 0. This statement can also be obtained on physical grounds, when remembering that eqn (4.1) describes injection of air from the bottom into a medium which is completely filled with water at t = 0. To assess the maximum value, denoted by S,, we find from the boundary condition at z=-h u
(4.1)
The saturation just above the filter is chosen because it appears to be stationary and maximum on the flow domain.
Multi-phase Jlow modeling of air sparging
327
Consequently, we obtain in terms of the air saturation the approximate problem ( V - kr,(NgZZ - N,Vp,) = 0
II
_Nk
c
M
s-1 &
in r > c, z < h, for r = E, -hl < z < -h,
-
for r = c, z < -h,, -h,
PC =Pt
I
Pc-$Pt-
< z < ht
for r > E, z = ht $(ht c
-z)
for
r-+cq r>c,
z
(4.5)
I
Fig. 7. Schematic of the domain (01) for air injection through a horizontal disk.
In view of the large difference between the densities of air and water we assume that upwards from the injection titer and below the water table, flow in the z-direction is dominated by advection, i.e.
4.2 Steady state flow
(4.6)
At steady state the tilme derivative in the saturation equation in problem (2.17) vanishes. The reduced equation must be solveid subject to the same boundary conditions. Inherited from the initial pressure distribution we impose in addition pc + pt - 2
c
(h, - z) for ’ ---f O”’
I > E,
z < h,
z + -_oo
(4.3)
The numerical results ;show that the water pressure is approximately hydrostatic at large times. This implies that the water velocity can be disregarded with respect to the air velocity, i.e. Zt = ii=. Using this in eqn (2.13) yields ii,, = N,k,,e’, - N,k,,Vp,
(4.4)
As a consequence we are left with a nonlinear diffusion equation in the variables r and z. When solving this equation, no boundary condition can be imposed at z = ht, representing a ‘future time’. Hence, we consider the flow domain above the well of infinite extent, represented by z > -h,. Interpreting solutions only up to z = 0 will produce a reasonable approximation in this region. Assuming furthermore that the well-radius 6 < 1, we arrive at the problem, setting k(r, z) = k,,(S(r,z))
in r > 0, z > -h, $O,z)
=o
for z > -h,
t k(cqz)
= 0
for z > -h,
0.45I--
0.80
0.40 .
0.70
0.35 -
0.60
0.30 : m
d ,I' /' /In /' ,' /' /' hux
0.15 -
(4
0.50 v)
St
0.20 -
0.050.20
/IP
/
0.25
0.10 -
(4.7)
0.40 0.30
P-
/
/'
0.30
0.40
0.20 0.10
0.50
0.60
m
0.70
0.80
0.90
0.0
1.00 (b)
%lU 2.0
4.0
6.0
8.0
10.0
._
12.0
N,
Fig. 8. Comparison of S, for u,,, = 1 and S,, as a function of (a) the exponent m and (b) the gravity number Ng. The m-curve reflects values of the parameters for cases 2-5 and the NE-curve for cases 6-9 of Table 1.
328
M. I. J. van Dijke et al.
Appendix B(k) - nDkPD
(4.11)
where 2(1 - d,(m-1)/(4m+l), 4m+l
nD=
2(1 -m) “= 0.00 0.00
0.10
0.20
0.30
0.40
0.50
0.60
0.70
0.80
0.90
Fig. 9.
Note PD <
k
4m+l
that
in
n,
>
()
(4.12) o
’
most
cases
of
practical
interest
l.
Hence, following Ref. 1 and 18 we find that for sufficiently small k the solution of eqn (4.7) with the conditions (4.9) and (4.10) can be approximated by eqn (4.13)
Graph of B(k) with m = 0.750.
(4.13)
where N1 = N,/N, and D(k) = kdp,/dk. In Fig. 9 we show an example of b(k). At z = --Jr, we need to impose an ‘initial condition’, which should reflect the influence of the injection filter. First, we observe, by a straightforward continuity argument that ccl ru,,,dr = A for all z 2 -h, 2lr (4.8) s0 Using this approximation rewritten in terms of k
00 J
rk(r, z) dr = F
implied by (4.6) this can be
for all 2 2 -h,
0
(4.9)
where N2 = A/TN,. For computational purposes we assume that the injection of air only takes place at the centre point of the filter, i.e. r = 0, z = - 1. This leads us to consider (4.7) and (4.9) for all z > -1, subject to the ‘initial condition’ k(r, -1) = 0
for all r > 0
(4.10)
In other words, the injection of air is described by a delta-distribution at the point r = 0, z = - 1. Only for special diffusion coefficients b(k) eqn (4.7) admits an explicit solution. For example when B is a power law function a solution is known in terms of similarity variables. This solution is known as the Barenblatt-Pattle point source solution.11’8 As a result of the upper bound derived in the previous section, we may conclude that k is small throughout the flow domain. We use this observation to approximate B(k) by a power law function, see the
where f(z) = fo . (z + l)‘@(PD + ‘1)
forz>-1
(4.14)
and fo = r$Ny)
1”2(pDi”‘(?$)j
(4.15)
Observe that the function f(z) represents the interface that separates the region where k, S > 0 from the region where k, S = 0. In the mathematical jargon the set {(f(z>,z); z 2 -1). IS k nown as the free boundary of eqn (4.7). Strictly speaking, the Barenblatt-Pattle solution (4.13) is valid for all z > - 1. However, since eqn (4.7) is formulated for z above the injection filter, i.e. z > -h,, we restrict the solution to these values of z. Furthermore, since (4.11) is a good approximation for small k, we expect that the quality of the approximation (4.13) increases with z. Using eqn (A4) the expression for k is easily converted into an expression for the saturation S.
5 APPLICABILITY OF THE ANALYTICAL STEADY-STATE APPROXIMATION The approximation for the steady state depends only on the three numbers Ni, N2 and m and is valid for z > -h,. We check whether the numerically obtained steady-state solution also depends on these three numbers for z > -h,. This is done by simulations 18 and 19, for which the input parameters are shown in Table 2. The other parameters are the same as in Section 3. Observe that although the numbers Nz, NC and A are different Ni and NZ are the same.
Multi-phase flow modeling of air sparging
329
Table 2. Parameters and dimensionless numbers for checking dependence of tbe steady-state solution on N1, Nz and m Case
ct (m--l)
Ui* (10e3 ms-‘)
H (m)
m
%
NC (10-3)
A
4.42 2.21
7.50 5.30
0.667 0.667
6.63 13.3
0.885 1.77
5.59 11.2
1.00 1.41.
18 19
Computations result in almost identical contourplots of air saturation at steady-state, i.e. at large time. However, the transient behaviour is different. This is illustrated in Fig. 10, where the stored air volume V(t) is plotted. We conclude that both steady states are identical. Therefore, also the numerical solutions at steady state depend only on m and the (combined) numbers Nt and N2 for -h, < z < 0. We investigated the agreement between the analytical and numerical solutions at steady state in terms of the relative permeability, k(r, z). For the numerical solution, k is calculated from the air saturation by the van Genuchten relation (2.5). For -h, < z < 0, the relative permeability contourplots of both the analytical and numerical solution for case 1 are shown in Fig. 11. Taking the cross-section of the solutions for z = -0.2, we use k,,(r) and knum(r) to represent the analytical and numerical solutions, respectively, as shown in Fig. 12. We find that the two profiles give excellent qualitative agreement. To obtain a quantitative criterion, we compute the masses involved, i.e.
z,, =
J
see (4.9), and
Note that Z,, depends. only on N1, N2 and m. To plot the relative error AZ = (I,, - Z,,,)/Z,, for various values of m, Ni, and Nz, we make use of Table 3, which reflects the same computations as shown in Table 1. AZ is shown in Fig. 13 as a function of the dimensionless numbers, which are normalized with respect to the reference case, indicated by the
0.010
case 19 __-------
T-
1 , -* I’
0.006 F=
,I’
0.004
0.000 0.0
2.0
4.0
6.0
8.0
10.0
t Fig.
10.
m
m,=--,
Nln
NI Nit
N2
N2n = -
=-9
ml
N21
for those cases in which only one of the numbers is varied. We observe that the relative error is nearly constant and small (less than 2%) for most cases, except for values of the parameters taken from case 17. In this case the injection depth, h, is small and as a consequence both the geometry of the filter and the outflow boundary condition have a relatively large influence on the numerical approximation. We may conclude that at steady state (T > 2 h) the complex two-phase flow process can be accurately approximated by the relatively simple expressions (4.13) and (4.14). An important qualitative feature of the solution is its radius of influence, r = f (0) = fo, see (4.15), which denotes the outflow radius of the injected air.
6 CONCLUSIONS
O” rkan(4 dr 0
0.008
subscript 1:
Comparison of stored air volumes for cases 18 and 19.
We modeled continuous air injection below the groundwater levels as multi-phase flow of two immiscible fluids, with emphasis on the development of the region containing air and the air saturation profile within this region. We found that the solutions depend on four dimensionless numbers, composed of soil, fluid and injection parameters. Using a numerical model based on the mixed form of the Richards equation for two phases, we made Table 3. Dimensionless numbers characterizing tbe steady state Case
m
N1
N2 (10-4)
1 2 3 4 5 6 7
0.667 0.400 0.500 0.750 0.800 0.667 0.667 0.667 O-667 0.667 0.667 0.667 0.667 0.667 0.667 0.667 0.667
8.99 8.99 8.99 8.99 8.99 3.59 4.49 6.74 13.5 11.2 7.49 4.49 2.25 8.99 8.99 8.99 8.99
7.45 7.45 7.45 7.45 7.45 18.6 14.9 10.0 4.96 7.45 7.45 7.45 7.45 2.68 3.97 4.98 12.3
: 10 11 12 13 14 15 16 17
M.
330 -O.OO
0.01 I
-0.11
N
0.11
0.21
0.31
0.41
I. J. van Dijke et al.
0.51 -O.OO
-0.00
0.01 r
t-
-0.44
-0.44
-0.56
-0.56
-0.56
-0.56
-0.67
-0.67
-0.87
-0.67
-0.76
-0.78
-0.78
-0.76
-0.89
-0.89
-0.89 0.11
0.01
L
-I 0.31
0.21
0.41
0.51
N
E
0.01
r
(4
-0.44
-0.89
computations for different values of the dimensionless numbers. The computations provided accurate quantitative information about the development of an air cone in the flow domain and about the corresponding saturations. For some cases, the numerical computations showed instabilities during the development of the air cone. Within a short time (about 2 h) a steady-state situation was reached, which was approximately a single-phase (air) flow process. For all cases air
0.41
0.51
r
(b)
Fig. 11. Comparison of the relative permeability contourplots of
0.31
0.21
0.11
the analytical and (b) the numerical solution (case 1).
saturations were small in the entire flow domain. Unfortunately, the numerical computations were very time consuming. An alternative for the use of the two Richards equations is the fractional flow formulation. This formulation involves one equation for the air saturation and one time-independent equation for the total velocity. It presents an understanding of the
O.1° I 0.024
,
0.016
0.008 -0.02
’ 0.20
J 0.40
0.60
0.07
0.14
0.21
0.28
0.35
0.42
r
Fig. 12. Cross-section of the relative permeability profiles k,,(r) and knm(r) at z = -0.2 (case 1).
1.00
1.20
1.40
1.60
1.80
normalized number8
0.000 0.00
0.80
AI of the volume integrals versus the dimensionless numbers. The m-cuNe reflects values of the parameters for cases 2-5, the Nl-curve for cases lo-13 and the N2-curve for cases 14-17 of Table 3. Fig. 13. Relative errors
Multi-phase flow modeling of air sparging nonlinearities of the equations and is used for derivation of an upper bound on the air saturation in the flow domain and an approximate analytical solution for the steady state. The observation of small air saturations enabled us to approximate the complicated expressions for the relative permeability and capillary pressure functions by simple power law functions. Assuming that flow in the vertical direction is dominated by advection, we derived the analytical solution for the steady-state situation which is valid for the flow domain above the injection filter. This analytical solution was proven to be in good agreement with the numerical results and it depends on three dimensionless numbers, which are a combination of the four numbers mentioned before. For practical reference we give the expression of the approximate analytical solution for the air saturation and the expression for the dimensionless numbers involved in this approximation. The numbers are N1
=
(Pw- PaWa
(6.1)
Pw N
=
2
&aUinE(HI- Hu) &bshv
-
(6.4
dgH2
and the exponent m in the nonlinearities. The analytical solution is given by
I
ns(“@$$))-(1
S(r,z) =
IO
-h,LzlO
for r > f(z),
-h,
where f(z) =fo
* (Z +
Table 4. Oufflow radii
Case 1 2 : :
FO Cm)
Case
2.31
7 8 9 10 11 12
1.61 1.87 2.53 2.65 3.65
FO @I
Case
FO Cm)
13 14 15 16 17 18
4.15 3.56 3.02 2.74 1.87 3.85
3.27 2.67 1.89 2.10 250 3.10
Observe that the coefficient ps/po = l/( 1 - m) in eqn (6.3) satisfies ps/pD > 1 for all 0 < m < 1, which means that the saturation profile in the r-direction always has zero derivative at the free boundary. Hence, the profile shows a relatively large zone where S is almost zero. Since for remediation purposes a minimum air saturation is necessary, we need to know of the air saturation profile to determine at what distance from the injection axis this minimum saturation is reached. In our calculations we imposed at the filter an influx, which is equally distributed over the length of the filter. However, in practical situations most air leaves the filter near the top.15 Because for the approximate analytical solution we assumed point source injection, the use of this approximation can be improved by setting the injection level higher than the centre of the filter.
ACKNOWLEDGEMENT
- (&)2r’pD
forOIrIf(z),
331
1)1’(2(pD+1)) for - h, 5 z 5 0
This project was sponsored by the Netherlands Organization for Scientific Research (NW0 project NLS 61251) and carried out in the framework of the Netherlands Integrated Soil Research Programme (PCTB project 5029), with the use of super computing facilities of the Dutch National Computing Facilities Foundation (NCF).
(6.4) with REFERENCES
represents the outer botmdary of the air cone. Note, that in eqn (6.3) dimensionless spatial coordinates r and z need to be inserted. To obtain FO, the radius of influence in a real physical dimension, f0 must be multiplied by the characteristic length H. Furthermore, 2 ns = m4ml(4m+l), Ps=4m+ 1
31 -m)m(m-k)j(4m+l),pD =W -4 4m+l Q= 4m+l For reasonable values of m, such as 0.30 < m < 0.90, we have O-52 < ns < O-92, 0.43 < ps < @90, O-043 < nD < 0.93 and 0.043 < po < 064. To give an impression of the width of air cones, we list the outflow radii F0 in Table 4 for the computations of Section 3.
1. Barenblatt, G. I., On some unsteady motions of a liquid or a gas in a porous medium. Prikl. Mat. Mekh., 16 (1952) 67-78 (in Russian). 2. Bear, J., Dynamics of Fluids in Porous Media. Elsevier, New York, 1972. 3. Bohler, U., Brauns, J., H&l, H. & Nahold, M., Air injection and soil-air extraction as a combined method for cleaning contaminated sites - observations from test sites in sediments and solid rocks. In Contaminated Soil ‘90, ed. F. Ahrendt. Kluwer Academic Publishers, Dordrecht, 1990, pp. 1039-44. 4. Celia, M. A., Bouloutas, E. T. & Zarba, R. L., A general mass-conservative numerical solution for the unsaturated flow equation. Water Resow. Res., 26 (1990) 1483-96. 5. Chavent, G. & Jaffre, J., Mathematical Models and Finite Elements for Reservoir Simulation. ed. J. L. Lions, Studies
in Mathematics and its Applications, Vol. 17. NorthHolland, Amsterdam, 1986. 6. Crow, W. L., Anderson, E. P. & Mint@, E. M., Subsurface
venting
of vapors
emanating
from hydro-
M. I. J. van Dijke et al.
332 carbon
product
on ground
water. Ground Water Mon.
Rev., 7 (1987) 51-7. 7. Forsyth, P. A. & Shao, B. Y., Numerical simulation of gas venting for NAPL site remediation. Adv. Water Resour., 14 (1991) 354-67. 8. Gilding, B. H., Qualitative mathematical analysis of the Richards equation. Transport in Porous Media, 5 (1991) 651-66. 9. Herrling, B. L Stamm, J., Numerical results of calculated
3D vertical circulation flows around wells with two screen sections for in situ or on-site aquifer remediation. In
nSkPs and nokPD, respectively,
where nj andpj (j = S, D) are constants depending on m. The relations (2.5) and (2.12) for the capillary pressure and the relative permeability are:
PC(S) = ((1 - S)+
k(S) = Si(l - (1 - S)1’“‘)2m Expanding
Computational Methods in Water Resources IX, Vol. 1: Numerical Methods in Water Resources, ed. T. F. Russel.
Elsevier Applied Science, London, 1992, pp. 483-92. 10. Ji. W., Dahmani. A.. Ahlfeld. D. P.. Lin. J. D. & Hill III. E: HI, Laboratory. study of air’ sparging: Air flow visualization. Ground Water Mon. Rev., 13 (1993) 115-26. 11, Johnson, P. C., Kemblowski, M. W. & Colthart, J. D., Quantitative analysis for the cleanup of hydrocarboncontaminated soils by in-situ soil venting. Ground Water, Johnson, P. C., McWorther, D. B., Goodman, I., An overview of in-situ air Water Mon. Rev., 13 (1993) 127-35.
J. & Parker, J. C., An efficient finite for modeling multiphase flow. Water
(1 - S)“”
(l@=l-f-S+;
in a binomial k-1
m
(
k(S) = m -2y+m+l+ = @k(S))
For &(S) Q(S)
(1995) l-27.
15. Marley, M. C., Hazebrouch, D. J. & Walsh, M. T., The application of in-situ air sparging as an innovative soils and ground water remediation technology. Ground Water Mon. Rev., 12 (1992) 137-45. 16. Mohtar, R. H., Wallace, R. B. & Segerlind, L. J., Finite element simulation of oil spill cleanup using air sparging. In Computational Methodr in Water Resources X, Vol. 2, ed. A. Peters. Kluwer Academic Publishers, Dordrecht, 1994, pp. 967-74. 17. Parker, J. C., Lenhard, R. J. & Kuppusamy, T., A parametric model for constitutive properties governing multiphase flow in porous media. Water Resour. Res., 23 (1987) 618-24.
18. Pattle, R. E., Diffusion from an instananeous point source with a concentration-dependent coefficient. Quart. J.
;S2+0(S3) >
O(S4”‘1)
W)
= k(dp,/dS)/(clk/dS)
we find
2(1- m)
=
m
Resour. Res., 25 (1989) 43-54.
14. Lundegard, P. D. & LaBrecque, D., Air sparging in a sandy aquifer (Florence, Oregon, USA): Actual and apparent radius of influence. J. Contam. Hydrol., 19
series
for k
we obtain
28 (1990) 413-29.
12. Johnson, R. L., Hinchee, R. E. & sparging. Ground 13. Kaluarachchi, J. element method
- l)‘-m
S(l - (1 - s>+r-” x (1 -S)““(l Again
- (1 -S)‘/m+4S(l
using the series (Al),
&(S)
=
-,)++I)
we obtain
2~~~~)mm-1S1pm + O(S2-2m)
w
which is only accurate for S close to zero, Taking first order terms we get
S(k) N m 4”l/(4m+
1) k2/(4t?Z+
1) (A41
Observe that these approximations are only good for small m, because when m r 1 the first and the higher
Mech. Appl. Math., 12 (1959) 407-9.
19. Proskurowski, W., A note on solving the Buckley-Leveret equation in the presence of gravity. J. Comput. Phys., 41
0.24 ,
(1981) 136-41.
20. Protter, M. H. & Weinberger, H. F., Maximum Principles in Dt&ential Equations. Prentice-Hall, London, 1967. 21. Smoller, J., Shock Waves and Reaction D@usion Equations, Springer Verlag, New York, 1983. 22. Wehrle, K., In-situ cleaning of CHC contaminated sites: Model-scale experiments using the air injection (in-situ stripping) method in granular soils. In Contaminated Soil ‘90, ed. F. Ahrendt. Kluwer Academic Publishers, Dordrecht, 1990, pp. 1061-2.
0.20 -
0.16 ________-----
________-----n
0.12 ____-_-__-
___---
’
0.00 0.00
APPENDIX According to the observation that for z > -h,, S(r, z) is significantly smaller than one, we approximate S(k) and b(k) = kdp,/dk by power law functions of the form
0.03
0.06
0.09
0.12
regression
0.15
0.18
0.21
k
Fig. Al. Approximations of b(k) with m = 0.750 for k < 0.151. Approximation by binomial series yields (nn,pn) = (0*127,0.125), whereas regression analysis yields (nn,pn) = (0.182,0:146).
Multi-phaseflow modeling of air sparging order terms of eqns (AZ!) and (A3) become of the same order of magnitude. We get ns = m h/(4m+ 2 “=4m+l
1)
7
2(1 - “),(m-:l)/(4m+l), nD= 4m+l “=
2(1 -m) 4m+l ’
O
0
WI
333
Note, that using the upper bound for S, i.e. S < S,, see (4.2), or equivalently k < l/N,, it is possible to find other power law functions by regression analysis, which approximate D(k) for 0 5 k 5 l/N,. Figure Al shows the function D(k) with approximations for m = 0.750, which must be valid for k 5 1/IVr = 0.15 1. Because the approximation by binomial series shows the correct limiting behaviour of b(k) for k JO, we use this approximation for the derivative of the analytical solution.