Numerical simulation of two-phase flow in a geocentrifuge

225

N u m e r i c a l simulation of t w o - p h a s e flow in a geocentrifuge B. Ataie-Ashtianil'3 S.M. Hassanizadeh 2, O. Oung 1, F.A. Weststrate l, A. Bezuijen 1 1Environmental Department, GeoDelft, Delft, The Netherlands 2Centre for Technical Geosciences, and Hydrology and Ecology Section, Delft University of Technology, The Netherlands 3Department of Civil Engineering, Sharif University of Technology, Tehran, Iran In this paper, results of DNAPL spreading experiments carried out in a geocentrifuge are analysed. The experiments are performed in the GeoDelft geocentrifuge facility in the framework of the EU-funded project NECER (Network of European Centrifuges for Environmental Research). The experiments constitute the intrusion of a known amount of PCE into a water saturated soil sample. They are carried out under 12g, 15g, and 30g. The objectives of these experiments were to investigate the potential of geocentrifuge tests for the study of environmental problems, to investigate the effect of acceleration on the spreading of DNAPL, and to determine whether pc-s curves of soils measured in a geocentrifuge are valid under 1-g conditions. Based on the experimental data and numerical simulations, it can be concluded that the capillary pressure-saturation curves measured in a geocentrifuge are valid for practical applications of multi-phase flow at 1-g. Furthermore, we have demonstrated the importance of numerical simulations in the understanding of measured data, through trying different scenarios and matching numerical and experimental results to describe the observed data and their usefulness in the design of a complex experiment on multiphase flow. 1. INTRODUCTION Geocentrifuges have been used for a long time to simulate a high gravitational acceleration by applying a centripetal force to a sample. Geocentrifuges have the ability to increase the body force on the sample and therefore decrease the required testing time. Mitchell (1998) has provided an extensive review of the use of centrifuge for geoenvironmental practice and education. One of the largest geocentrifuge facilities is located in GeoDelft, The Netherlands. The arm of the centrifuge is 6 m long and an acceleration of 300g can be reached. Characteristics of the GeoDelft geocentrifuge have been described by Nelissen (1991). Recently, a series of experiments have been carried out in this geocentrifuge in the framework of the EU-project NECER (Network of European Centrifuges for Environmental Geotechnic Research). To the authors' knowledge, this is the first time that a geocentrifuge has been used for obtaining Pc-S curves for a DNAPL-water system. The purpose of this paper is to present our efforts in modelling the experimental results and to show how the modelling work is used to investigate relevant physical processes and to identify possible shortcomings of the experimental design. It appeared that the simulation of the geocentrifuge experiments was not a straightforward task. This was because essential data on boundary conditions were not available and measurements of pressure transducers needed to be properly interpreted before they could be used. Thus, in this work we had to

226 make educated guesses and expert judgment in setting up our simulation work. The simulation procedure and scenarios developed have helped to gain insight in the physical processes involved in geocentrifuge experiments and have made it possible to interpret the results properly. Based on these simulations, problems and shortcomings of the current setup are explored and suggestions for its improvement are made. Also based on these experimental results and with the use of numerical simulation, we can conclude that measurement of constitutive relationships for two-phase at a higher g-level can be used for studying of NAPL transport at 1-g acceleration.

2. EXPERIMENTAL SET-UP AND RESULTS Here, a brief description of the experimental set-up is given; further details can be found in Oung et al. (2000). Schematic presentation of the column set-up is shown in Figure 1. In these experiments, the displacement of water by PCE and its subsequent displacement by water in a column are studied. It is a 10-cm cylinder consisting of three parts: the inflow reservoir, a soil compartment and the PCE collection reservoir. Displaced water from the PCE collection reservoir may flow over to the water collection reservoir (see Figure 1). The whole system is placed in the geocentrifuge container. As the geocentrifuge reaches its designated maximum speed, the container tilts horizontally until the column assumes an almost horizontal position and the desired g-level is reached at the soil sample surface (position of PPT1 transducer in Figure 1). Instrumentation of the set-up consisted of a combination of hydrophobic Pore Pressure Transducers (PPT) and Tethra Probe Frequency Domain sensors (TP) for measuring water saturation. There are seven PPTs and three TPs. The positions of the transducers, measured from the bottom of the DNAPL collection reservoir, are given in Table 1. The water level in the inflow reservoir was kept constant at 59.3 cm. The Zeijen sand was used in these series of experiments. The sand column length is 18 cm. Initially, one liter of Perchloroethylene (PCE) is injected in the inflow reservoir. Pertinent sand properties, including permeability and Brooks-Corey (1964) parameters are listed in Table 1. The parameters were obtained by static measurement of capillary-pressuresaturation curves under quasi-static conditions (Oung and Bezuijen 2001). In this series of experiments, three different acceleration levels, namely 12g, 15g, and 30g, were applied. As all systems were below the water level, hydrostatic water pressure was subtracted from the raw signals of all PPTs. All pressure transducers were tested and calibrated before installation. However, there is always some shift in the readings of the pore pressure gauges. This is caused by forces on the pore pressure gauge due to their mounting and to the increased g-level during spinning up. It is therefore common practice to try to start an experiment with a well-defined pressure distribution. This was also done in these experiments. The hydrostatic pressure before the release of the DNAPL was used as a starting point. Thus, the correction of PPT readings for hydrostatic pressure also accounts for pressure shifts due to the mounting and the increased g-level. As pressure transducers were hydrophobic; i.e. they measured the non-wetting phase pressure, Pnw, once mobile PCE reached them. Before that, they measured water pressure. After the passage of PCE, however, they are in contact with a mix of water and residual PCE and therefore, that reading does not represent any of the pressures.

227

.•

PLUNGER

lATERTABLE

INJECTION

~'11 I II

JET BREAKER PPT2

~

PPT3 m ~ SOILCHAMBER ~ PPT4 " ~

RESERVOIR

THETAPROBES (TP 1, TP 2.TP3)

~ t

m6 q ~ MODELSAND DNAPL

i WATER

Figure 1. Sketch of the test column used in centrifuge experiments. Table 1. Sand and fluid properties

Units

Zeijensand

[m2]

3.0x10 "12

Porosity, n

[-]

0.375

Entry pressure, Pd

[m]

0.56

Pore size distribution index, k

[-]

6.3

Residual saturation, Sr

[-]

0.134

[kg/m3]

1623

PCE viscosity

[Pa.s]

3.15x10 4

PCE-water interfacial tension

[N/m]

44.4x10 "3

Property

Permeability, k

PCE density

Obviously, in this experiment the column undergoes a stage of water drainage when PCE invades and a stage of water imbibition as PCE leaves the sand. The hysteresis effects have not been included in the numerical simulations because the purpose of these simulations is to capture the general trend of process rather than matching all details.

228 3. N U M E R I C A L SIMULATIONS The two-phase flow simulations are carried out using the simulation package STOMP (Subsurface Transport Over Multiple Phases). The simulations discussed in this paper were conducted with a non-hysteretic two-phase version of the water-oil mode of STOMP (White and Oostrom, 1997). The multi-phase flow in the test column is conceptualized to be onedimensional. The case of Zeijen sand at 15g is described first. The procedure for other cases is similar. Properties of the Zeijen sand, including permeability and Brooks-Corey parameters as well as fluid properties are listed in Table 1. Relative permeability-saturation relations were computed based on the Burdine (1953) pore-size distribution model. The DNAPL pressure at the top of column is based on the PPTI reading. At the same time, no-flux boundary condition is imposed for water. When PPTI readings indicate that there is no PCE left (at time t = 2190 s), the top boundary condition is changed to a specified water pressure for water and no-flux boundary condition for PCE. The water pressure at this time is calculated based on the water level of the inflow reservoir. The water and PCE pressures at the bottom boundary are assumed to be equal (i.e. zero capillary pressure) and are based on the free water surface in the flow line containing valve V2. Initial condition is calculated based on the level of water and PCE in the inflow reservoirs. a) b)

~ 1.4 i.2

~90

>l

r

70

i

50 30

9

.,

,

1000

,

i . .ureA2L

2000 3000 Time, s

--TP

A2

4000

c)

~0.6 0.4 ~0.2

~

o

5000

d) 15 t

!"

ll~w

~

,

I 0

, 200

1200

2200 time, s

3200

200

! 200

2200 time, s

3200

Figure 2. 15g experimental and numerical results of, (a) water saturation at TP2, (c) accumulated volume of PCE and water inflow and outflow, (c) PCE pressure at PPT3, and (d) PCE pressure at PPT5. The soil column is discretized into a one-dimensional grid with a constant nodal spacing of 2.5 mm. Time step varies automatically from 1 s to 10 s. Comparisons of simulation results and experimental data are shown in Figure 2. There is a reasonable agreement between measured and calculated water saturation values. In Figure 2b, calculated and measured values of PCE outflow volumes (NWout and EXP.NW'out, respectively), water

229 outflow volumes (Woutand EXP.Wout) as well as calculated values of inflow water and PCE volumes (Wi, and NWi,, respectively) are shown. Here, there are considerable discrepancies between calculated and measured outflow volumes. Finally, pressure curves at PPT3 and PPT5 are shown in Figures 2c and 2d, respectively. While the result of simulation for PPT3 (at least the trend) is somewhat acceptable, the calculated pressure curve at PPT5 (near the bottom boundary)is completely off and nowhere near the measured curve. Also, the total PCE influx volume was one liter whereas the calculated value is 1.40; thus 40% more. As it was shown in Figure 2, the agreement of numerical results with experimental ones for Pnw deteriorated at positions closer to the bottom of the sand column. This trend pointed to the possibility that the boundary condition may be inconsistent with the actual situation. The bottom boundary condition imposed here was that fluids might flow out freely into the collection reservoir, at zero capillary pressure. The water and PCE pressures at the bottom boundary were assumed to be equal (i.e. zero capillary pressure) and were based on the free water surface in the flow line containing valve V2. However, head measurements in the DNAPL collection reservoir and in the flow lines behind valve V2 show that there is a considerable pressure drop in the flow line. This might be due to flow resistance caused by regulation valve V2 and other frictional effects. In order to include such effects, another set of simulations were performed as described below. a) b)

"~ ra~

-t E~ TM

70 7

> ! ~'0.8 I

s0

~0.4 1-

9~ , . 2

90

o

/"

~061 ..... ,I ~..... TP N~A2

,

30 0

............ . . . . . . . . . .

1000

2000

3000

4000

L

'//-

-- ~ ' - - ~ - - - ' ? ~ - - i l~lWin ' , . . . . I ~'--. . . . . . .., - : ,

......... r

//~Wp,,l

I

.,f~-',-

',

.....

I_/-.L,,

50O0

0

ua..". ~ '

Ey.p.

- --

"J;=

i

,, 200~ime, 3000

l O00

i 4000

50(0

Time, s

c)

d) 15

~o

......

.~10

'

j 0

0 200

1200

2200 time, s

3200

'

'

"

........... i

200

'

I

1200

,

2200 time. s

i

3200

Figure 3. 15g experimental and modified numerical results of (a) water saturation at TP2, (b) accumulated volume of PCE and water inflow and outflow, (c) PCE pressure at PPT3, and (d) PCE pressure at PPT5. The resistance to flow is simulated by adding a virtual layer of porous medium to bottom of the sample with a lower intrinsic permeability and a lower entry pressure. The lower entry pressure is chosen because this layer is not supposed to be a barrier to the flow of D N A P L . For the same reason, relative permeabilities of the bottom layer for both phases (water and PCE) are assumed to equal 1 at all times (not a function of saturation).

230

a)

b) -... _

. 2

~

NWin,,/

1.6

g

.,..-I"-

= i.6 "6 > ~.! r162.2

~.2

~

~0.8

0.8

~ 0.4

~ 0.4o,

,/1

1200

w~o..-" , , 1600

~

2000e,Tims 2400

':

o ,,.f,/ 1200

2800

d)

c)

,

1600

2000 e,Tim s 2400

2800

3~

30 25 ~20

g20

NI5

~15

1~

I !

t3

i

glo

0

,t" 5 0

0

-5

-5

1000

! 500

2000 time, s

2500

! 500

1000

3000

2000 time, s

2500

3000

e) 30

!,

30

.

.

.

.

.

.

.

.

i

25 ~20

vr5 Num~,'ical

_~!0

,,

i

......

gmo

0

0

-5

-5

1000

1500

2000 time, s

' Lj~,,

~5

2500

1000

3000

g)

1500

2000 time, s

2500

30O0

h)

~.90

~90 g,

.~_

70

70

t~

30

30 1000

i 500

2000 Time, s

2500

3000

1000

1500

2000 Time, s

2500

3000

Figure 4. 30g experimental and numerical results of (a) accumulated volume of PCE and water inflow and outflow, (b) accumulated volume of PCE and water inflow and outflow for modified simulation, (c) PCE pressure at PPT3, (d) PCE pressure at PPT3 for modified simulation, (e) PCE pressure at PPT5 and (f) PCE pressure at PPT5 for modified simulation, (g) water saturation at TP2, and (h) water saturation at TP2 for modified simulation. A large number of simulations, with various combination of intrinsic permeability and entry pressure were performed. The best result were obtained for the case that the bottom

231 porous medium is chosen to be 2.5 mm thick with a permeability 50 times smaller and an entry pressure 10 times smaller than the sample. The results are shown in Figure 3. It is evident that improved simulation results are obtained. The simulated and measured values of PCE outflow volumes (curves NWout and EXP.NWout, in Figure 3b, respectively), water outflow volumes (curves Wout and EXP.Wout, respectively) are quite close in this case. Also, the calculated volume of inflow PCE (curve NWin) approaches one liter, which is the amount that was infiltrated in the experiment. Obviously, this conceptualization of bottom boundary effects on the water and PCE movement is just a very rough approximation of reality. It is not clear how this resistance changes in time depending on the outflow of water and PCE. Commonly, the head loss in a pipe is assumed to be proportional to the square of velocity. The pressure loss cannot be exactly formulated and implemented in this problem. By adding a layer of porous medium to the bottom where Darcy's law applies, we effectively assume that the pressure loss in each phase is proportional to flux of that phase in the bottom of sample. Nevertheless, numerical results show that this conceptualization has been successful to a large extent. Based on these numerical analyses, for any further experiments, the bottom arrangement should be modified so as to have a controlled boundary condition for water and DNAPL flow. Also, to increase the accuracy of simulations, more controlled boundary condition at the top of sand column would be desirable. For example, maintaining a constant DNAPLlevel will lead to having only drainage and avoids hysteretic conditions. Accurate measurement of water pressure along the sand column is a necessity for any detailed study of phenomena like dynamic effects and effects of micro-heterogeneities on pc-s-k r. Therefore, pressure transducers for water pressure measurement should be installed. 4. G E O C E N T R I F U G E APPLICATION FOR MEASURMENT OF CONSTITUTIVE RELATIONSHIPS OF TWO-PHASE F L O W In spite of the shortcomings of the collected data in these experiments, we believe they can provide insight in the utility of geocentrifuge experiments for obtaining two-phase constitutive relationships. We try to reach such a conclusion by induction. We use the same Pc-S curve for simulating various experiments carried out at different g-levels. If we obtain satisfactory agreement with measurements, then we conclude that the measurement of pc-s curve at a higher g-level in a geocentrifuge can be used for 1-g simulations. Figures 4 present the results of 30g experiment and simulation without and with the effect of bottom resistance. As seen in this figure, with the same P~-S curve and imposed boundary condition, the general behavior of PCE movement has been reproduced in numerical simulations. The calculated saturation-time curves are similar to experimental results and also the amount of PCE entered from the top of sample, 1 liter, is re-produced well by numerical simulations. It should be noted that because of malfunctioning of PPT6 and PPT7 readings in these cases, the correct values of the water and PCE outflow from the experimental data could not be calculated. Numerical simulation results for outflow of water and PCE are reasonable. Also the pressure measurement is simulated reasonably well for all g-levels. Overall, the agreement between the numerical model and experimental data are reasonably good, considering the fact that in the numerical simulation secondary processes like hysteresis and non-wetting phase entrapment have not been included. We know that entrapment is happening in these experiment and about 15 % of PCE remains behind in the

232 sample. Also the bottom boundary condition in these experiments is not very well controlled and it can not be exactly simulated. The effect of this boundary condition, however, has been roughly included in the model. All in all, based on this results it can be concluded that the measured Pc-S curve in geocentrifuge is valid for practical applications under 1-g. This means that the capillary pressure measured in the geocentrifuge is independent of the acceleration (g-level). Given the fact that the weight (centrifugal force) is directly proportional to N, we can easily establish that capillary pressure is independent of the g-level. 5. SUMMARY AND CONCLUSIONS A series of two-phase flow experiments has been performed in the GeoDelft geocentrifuge facility in order to study the process of DNAPL spreading in a water saturated porous medium. The experiments constitute the intrusion of a known amount of PCE into a water saturated soil sample. They are carried out under 12g, 15g, and 30g. In this work, a numerical model was employed for the interpretation of the experimental data. Numerical simulations indicated that some of the features observed in collected data are due to uncertainties in defining the bottom boundary condition. In an attempt to better simulate this boundary condition, a virtual layer is added to the bottom of the modeling domain. This layer has a lower permeability and a lower entry pressure than the sample. These results are in much better agreement with the experimental data. In spite of shortcomings in the experimental data, we have been able to simulate experimental data reasonably well, considering that in the numerical simulation secondary processes like hysteresis and non-wetting phase entrapment have not been included. All in all, by reproducing the experimental results of different g-levels using a unique constitutive relationship measured at l g, it can be concluded that the measured pc-s curve in geocentrifuge is valid for any practical application at lg. The common practice in almost all field-scale studies, or even laboratory studies is that these secondary processes are not included. Nevertheless, further investigation with modifications in the experimental set-up and procedure is required to reveal the significance of these secondary effects. 6. A C K N O W L E D G M E N T S This research has been carried out in the framework of the TRIAS Project "Upscaling micro-heterogeneities in two-phase flow in porous media" (Delft Cluster Project no. 5.3.1). REFERENCES 1. R.H. Brooks and A.T. Corey, 1964. Hydraulic properties of porous media. Hydrol. Pap. 3, Civ. Eng. Dep., Colo. State Univ., Fort Collins. 2. N.T. Burdine, 195. Relative permeability calculations from pore-size distribution data. Petroleum Transactions, 198:71-77. 3. R.J. Mitchell, 1998. The eleventh annual R.M. Hardy Keynote Address, 1997: Centrifugation in geoenvironmentalpractice and education. Canadian GeotechnicalJournal, 35: 630-640. 4. H.A.M. Nelissen, 1991. The Delft geotechnical centrifuge. "Centrifuge '91. H.Y. Ko and F.G. McLean, eds., Balkema, Rotterdam, The Netherlands, 35-42. 5. O. Oung and A. Bezuijen, 2001. Capillary pressure-saturation curves: Experimental determination, GeoDelft report no. DC-750301/350. 6. O. Oung, A. Bezuijen and F. Weststrate,2000. Centrifugeresearch on the transport of DNAPLs: 1-D flow of a DNAPL through water-saturated sand. IIIInt. Symp. Physical Modelling and Testing in environmental Geothenics, La Baule, France, 343-350. 7. M.D. White and M. Oostrom, 1997. STOMP, Subsurface Transport Over Multiple Phases. User's Guide. Pacific Northwest National Laboratory Report PNNL-11218