On the onset of convection in a water-saturated porous box: Effect of conducting walls

IN H E A T A S D M A S S T R A N S F E R 0094-4548/78/ii01-0371502.00/0 Vol. 5, pp. 371-378, 1978 © P e r g a m o n P r e s s Ltd. Printed in GreatBritain

ON THE ONSET OF CONVECTION IN A WATER-SATURATED POROUS BOX: EFFECT OF CONDUCTING WALLS *

Robert P. Lowell 1.

2.

1

and Chuen-Tien Shyu

2

School o f Geophysical Sciences, Georgia I n s t i t u t e o f Technology Atlanta, Georgia 30332 School of Oceanography, Oregon State University, Corvallis, Oregon 97331

(Communicated by J.P. Hartnett and W.J. M/nkowycz)

ABSTRACT The c r i t i c a l Rayleigh number f o r the onset o f convection in a w a t e r - s a t u r a t e d porous box, heated from below, has been d e t e r mined by the G a l e r k i n method and by an approximate a n a l y t i c a l technique. The present problem d i f f e r s from previous ones which have been published in t h a t we consider the e f f e c t o f conducting v e r t i c a l boundaries on the c r i t i c a l number and on the c e l l p a t t e r n at the c r i t i c a l number. Moreover, we emphasize f a u l t or f r a c t u r e z o n e - l i k e box geometries; t h a t is , geometries in which one box dimension is much s h o r t e r than the o t h e r and much s h o r t e r than the h e i g h t . The r e s u l t s i n d i c a t e that f o r f a u l t / f r a c t u r e zone geometries, in which the long v e r t i c a l w a l l s are conducting and the short end w a l l s are i n s u l a t e d , the c r i t i c a l Raylelgh number is roughly four orders o f magnitude g r e a t e r than the c r i t i c a l number f o r a h o r i z o n t a l porous slab. The f l o w at the onset o f convection takes the form o f a r o l l w i t h i t s axis along the long h o r i z o n t a l dimension o f the box. There i s , however, l i t t l e d i f f e r e n c e between the c r i t i c a l numbers f o r two and three dimensional c e l l p a t t e r n s . These r e u s l t s i n d i c a t e t h a t convection may occur in n a t u r a l l y o c c u r r i n g f a u l t s and f r a c t u r e zones in the e a r t h ' s crust o n l y i f the p e r m e a b i l i t y is of the order o f Darcies. In n a t u r a l systems, the Rayleigh number would probably not exceed the c r i t i c a l number g r e a t l y , and the f l o w may tend to be f u l l y three dimensional. *This manuscript is submitted for publication with the understanding that the United States Government is authorized to reproduce and distribute reprints for governmental purposes. 371

372

R.P. Lowell and C.-T. Shyu

Vol. 5, No. 6

Introduction The onset o f convection in a closed r e c t a n g u l a r c o n t a i n e r o f w a t e r s a t u ra t e d porous m a t e r i a l w i t h uniform p r o p e r t i e s has been t r e a t e d by Beck (1), and the case o f a v a r i a b l e v i s c o s i t y f l u i d Zebib and Kassoy (2).

has been t r e a t e d by

These authors have considered v ar ious box geometries,

but they have assumed t h a t the v e r t i c a l w a l l s o f the box are i n s u l a t e d . In many p r a c t i c a l s i t u a t i o n s , however, such as in geothermal systems c o n t r o l l e d by f a u l t and f r a c t u r e zones, the assumption o f i n s u l a t e d v e r t i c a l w a l l s is not r e a l i s t i c .

We wish to consider the e f f e c t o f

conducting w a l l s on the c r i t i c a l at the c r i t i c a l

Rayleigh number and on the c e l l p a t t e r n

Rayleigh number.

Moreover, we w i l l

f r a c t u r e c o n t r o l l e d geometry; t h a t i s , we w i l l

emphasize a f a u l t or

emphasize porous c o n t a i n e r

geometries in which one h o r i z o n t a l dimension is much s h o r t e r than the o t h e r and much s h o r t e r than the h e i g h t .

We w i l l

assume t h a t the long

w a l l s are conducting but t h a t the s h o r t , end- w alls are i n s u l a t e d . will

then discuss the p o s s i b i l i t y

f a u l t s or f r a c t u r e s

We

o f convection t a k i n g place in n a t u r a l

in the earthJs c r u s t .

Davis (3) and Catton (4) have t r e a t e d the s i m i l a r problem o f convection in a r e c t a n g u l a r c o n t a i n e r having conducting w a l l s and being filled

w i t h viscous f l u i d . Analysi.s Consider a c o n t a i n e r o f w a t e r - s a t u r a t e d porous m a t e r i a l o f h e i g h t L

and w i t h h o r i z o n t a l dimensions L 1 and L2.

Let a uniform thermal g r a d i e n t ,

, be a p p l i e d to the m a t e r i a l such t h a t the temperature at the upper surface is T = T (TO >Ts).

and the temperature at the lower surface is T = T s o The i n i t i a l s t a t e is assumed to be one o f pure thermal

conduction such t h a t i n i t i a l l y v = O, T = Ts +

Bz,

dP*/dz = Pfg

where ~ is the velocity, P* the pressure, Pf the density of the fluid and g the acceleration due to gravity.

The z axis assumed to be positive

downwards. The condition for the onset of convection is found by the principle of exchange of stabilities.

That is, infinitesimal perturbation of the

initial state leads to a marginally stable state of steady thermal convection.

The linearized, non-dimensional, marginal stability equations,

usina the Boussinesq anproximation are (5)

VOI. 5, NO. 6

~ I O N

I%.IA WAT~-SATURATED

P O E U S BOX

V.q : 0 + R½0k + VP : 0

373

(I) (2)

v20-R½qz = 0 (3) + where q = (qx,qy,qz) is the perturbation v e l o c i t y , 0 is the perturbation temperature, P is the perturbation pressure, and k is a u n i t vector in the z direction.

R is the dimensionless Rayleigh number given by R = pfs~SKgL2/vX

w h e r e s is t h e s p e c i f i c coefficient, k

its

ature

v its

thermal

heat of

kiematic

conductivity.

the fluid,

viscosity,

~ its

thermal expansion

K the permeability

These e q u a t i o n s

of

t h e r o c k and

a r e t o be s o l v e d w i t h

temper-

conditions, a e / a x = 0 a t x = 0,H 1 e = 0

(4) (S)

a t y = 0,H^

z

0,I z

where H I and H 2 are dimensionless horizontal

lengths.

Assuming impermeable

boundaries, the velocity conditions are .

n =

0

(6)

^

on a l l boundaries, where n is the unit vector, normal to the wall. Equations

(I) through (3) with conditions

(4) through (6) constitute an

eigenvalue problem for the Rayleigh parameter, the smallest eigenvalue being the critical Rayleigh number for the onset of convection. The set of equations and boundary conditions above are not separable, and consequently the normal solution procedures break down. One method of determining the critical Rayleigh number for problems of this type is the Galerkin method.

To use this method, we first eliminate

and P from the perturbation equations to obtain a single equation for e. We obta in V40 + R~Th20 = 0 where Vh2 is the horizontal Laplacian.

(7)

We then select a set of t r i a l

functions 0. which satisfy the temperature and the v e l o c i t y boundary J

conditions.

We substitute N

O=

Z

i=1

C.O. I

I

into (7) and apply the Galerkin c r i t e r i o n .

That is,

374

R.P. LOwell and C.-T. Shyu

Vol. 5, No. 6

N

! I C [ V2 (V20) + RVh2ei] .OkdV = 0

i I

i

(9)

i

where the integral (9) is carried out over the volume of the box, and 8 k, is one component of the set of trial functions

8.. I

Substitution of the

particular trial functions into (9) gives rise to a matrix equation of the form MAT = RMB~

(IO)

and the condition for a non-trivial solution is that the determinant of the coefficients be zero.

That is: MA -1M B - R -11

= 0

Equation (II) represents an eigenvalue problem for the Rayleigh number R; and the minimum eigenvalue found, R = Rc, represents the critical condition for the onset of convection.

The results for the critical

Rayleigh number, for different box geometries are given in Table I. TABLE I Critical Rayleigh Numbers 2-D and 3-D Convection in a Porous Box With Two Vertical Conducting Walls

H2~

2-D (i=O)

O.01

0.05 O.1

446831.6 17485.16 4327.77

HI~

3-D I

5

I0

446835.7

446832.1

446832.1

17488.8

17485.4

17485.3

4331.35

20

446832.1 17485.26

5O 446~32. I 17485.2

4327.92

4327 81

4327.75

4327.78

212.68

212.69

0.5

212.75

216.34

212.90

212 70

1.0

83.22

86.07

83.48

83 36

83.34

83.32

1.5

60.50

61.73

60.58

6O 52

60.50

60.50

40.66

40.66

3

43.57

40.85

40.86

4O 61

5

41.01

39.61

39.62

39 62

39.60

39.60

I0

39.87

39.48

39.48

39 49 39 48

39.49

39.49

39.48

39.48

5O

39.48

39.48

39.48

Vol. 5, NO. 6

CONVEC'gICN IN A ~ATER-SATURATED POROUS BOX

375

Discussion T a b l e 1 shows s e v e r a l with

interesting

two c o n d u c t i n g w a l l s .

dimensions,

the wall

First,

conditions

it

features of convection is c l e a r

if

the c o n d u c t i n g w a l l s

for

are far apart

f o r an i n f i n i t e

(H25 2 ) ,

between the i n s u l a t e d w a l l s

number

slab.

Secondly,

t h r e e d i m e n s i o n a l m o t i o n is

o v e r two d i m e n s i o n a l a t t h e o n s e t o f c o n v e c t i o n ;

separation

in a box

l a r g e box

become i m m a t e r i a l and t h e c r i t i c a l

approaches t h e v a ] u e 4 ~ 2 g i v e n by Lapwood (5)

preferred

that

increases

(HI+~),

and as t h e the three-

d i m e n s i o n a l v a l u e s o f R approaches t h e t w o - d i m e n s i o n a l v a l u e s . This c s u g g e s t s t h a t in n a t u r a l systems the m o t i o n s may tend t o be f u l l y t h r e e dimensional

rather

than r o l l - l i k e

suggest that at finite solutions

in c h a r a c t e r .

amplitude (i.e.,

Moreover, these results

R >R c) the

may be u n s t a b l e and t h e m o t i o n a t f i n i t e

three-dimenslonal.

Lastly,

two-dimensional a m p l i t u d e may become

in the case o f a f a u l t / f r a c t u r e

zone g e o m e t r y

(H 2 <<1, Hi 3 1 ) , T a b l e 1 shows t h a t Rc i s s e v e r a l o r d e r s o f m a g n i t u d e greater

than t h e v a l u e 4~ 2 f o r an i n f i n i t e

t a k e s t h e form o f a r o l l d i m e n s i o n o f t h e box. showed t h a t , para]lel

for

with This

its

axis parallel

is c o n t r a r y

insulated walls,

slab,

and t h e m o t i o n a t R = R c t o t h e long h o r i z o n t a l

t o the r e s u l t s

the r o l l s

their

axes

t o t h e s h o r t s i d e o f t h e box.

The r e s u l t s

in T a b l e 1 f o r a f a u l t / f r a c t u r e

zone geometry a r e in

r e a s o n a b l e agreement w i t h an a p p r o x i m a t e a n a l y t i c a l Lowell

(6).

velocity

The a n a l y t i c a l

calculations

boundary c o n d i t i o n s

condition

result

d e r i v e d by

were p e r f o r m e d by r e l a x i n g

on t h e v e r t i c a l

walls.

the

I f one r e p l a c e s the

qy= 0 a t y = O, H2 w i t h ~qz/~Z = 0 a t y = O, H2, t h e e i g e n f u n c t i o n s

become s e p a r a b l e , and e q u a t i o n this

o f Beck (1) who

were a l i g n e d w i t h

case, w i t h H 1 + ~ ,

(7) can be s o l v e d by s t a n d a r d methods.

H2<< 1, L o w e l l

In

(6) o b t a i n s :

R = 4~ 2 /H 2 c 2 With H2 = 0 . 0 1 ,

Rc

4 x 105

w h i c h is a b o u t ten p e r c e n t s m a l l e r than

t h e v a l u e found by t h e G a l e r k i n method. e x p e c t e d in v i e w o f the f a c t restrictive

in t h e a n a l y t i c a l

The a n a l y t i c a l a roll

with

its

results

that

T h i s u n d e r e s t i m a t e o f R i s t o be c the boundary c o n d i t i o n s were l e s s

calculations

a l s o show t h a t

axis along the strike Geophysical

T a b l e 1 shows t h a t material

with a fault

than in t h e G a l e r k i n method.

t h e f l o w a t R = R t a k e s p l a c e as c of the fault.

Implications

in a c l o s e d c o n t a i n e r o f w a t e r - s a t u r a t e d

or fracture

zone g e o m e t r y , e . g . ,

porous

w i t h H2 = 0 . 0 1 ,

the

376

R.P. ~ i i

and C.-T.

Shyu

critical

R a y l e i g h number is n e a r l y 4.5 x 10 5 .

possible

in n a t u r a l

and a f a u l t

systems?

5, No. 6

Are such h i g h R a y l e i g h numbers

Assuming a t e m p e r a t u r e g r a d i e n t

8 = iO0°C/km

d e p t h o f 3 km, which a r e r e a s o n a b l e f o r a g e o t h e r m a l r e g i o n ,

we e s t i m a t e v a l u e s o f the p h y s i c a l

p a r a m e t e r s from S t r a u s and S c h u b e r t ( 7 ) :

pfs = I cal/cm 3 - °C,~ = 1.5 x 10-3/°C,v -

Vol.

°C - s e c .

e,v,l

these parameters

are average

into the Rayleigh

onset of convection

not unrealistic

separation

for

fault zone.

the p e r m e a b i l i t y

fractures

within

due t o a s e r i e s o f f l a t h, Bear (8, p.

= 2 x 10-3cm/sec, ~ = 5 x 10-3cal/cm

over the depth

range.

Substituting

number gives K =3 x 10-8cm 2, for the

in a 50 m wide

value,

from i n t e r c o n n e c t e d permeability

values

This

is quite a high,

in a f a u l t

but

zone m a i n ] y a r i s e s

t h e r o c k volume; and assuming a parallel

fractures

o f w i d t h d and

164) has shown t h a t K = d3/12h

This,

if

on t h i s

K = 3 x IO-8cm 2, h = 100 cm, then d = 3.3 x I0-2cm. s c a l e do n o t appear u n r e a s o n a b l e f o r a f a u l t

the earth's

crust.

in such systems, critical

that

the a c t u a l

gradient

we emphasize t h a t

a l o n g the w a l l s

of convection in t i m e ,

R a y l e i g h number l i e s

zone in

convection occurs rather

close to the

the a s s u m p t i o n o f a u n i f o r m t h e r m a l

is o n l y a p p r o p r i a t e

in a n a t u r a l

the g r a d i e n t

be m a i n t a i n e d .

fault

or fracture

insulated

in the a d j a c e n t

In t i m e t h e r e w i l l

a c r o s s the w a i l s ,

as a c o n d i t i o n zone.

for

boundaries.

convection

f l o w may d e c l i n e

For t h i s

is

initiated.

As c o n v e c t i o n proceeds

be a d e c r e a s e in t h e l a t e r a l

porous c o n t a i n e r w i t h

and perhaps t o change i t s

cell

is reduced.

t i m e dependent a s p e c t s o f c o n v e c t i o n

in a

An e x a m i n a t i o n o f t h e s e p r o b l e m s

research.

r e s e a r c h was s u p p o r t e d by t h e U. S. G e o l o g i c a l Survey, Interior

In

may be somewhat more c o m p l i c a t e d

Acknowledgments

Department o f

pattern,

however, t h e c o n v e c t i v e

temperature gradient

insulated walls.

be a s u b j e c t o f f u t u r e

gradient

r e a s o n , one m i g h t e x p e c t c o n v e c t i o n

conducting walls

than in c o n t a i n e r s w i t h

not

tend t o have the appearance

At l a r g e t i m e s ,

s i n c e the d r i v i n g

summary, t h e f i n i t e - a m p l i t u d e ,

This

the o n s e t

impermeable r o c k m a t r i x w i l l

and hence the w a l l s w i l l

t o become somewhat more v i g o r o u s ,

will

if

a t v a l u e s o f R>Rc, and heat is conducted t h r o u g h the v e r t i c a l

boundaries,

after

or fracture

however, t h a t

number.

Finally,

of

One m i g h t e x p e c t ,

Fractures

under USGS Grant No.

14-O8-0001-G-365.

VOl. 5, ~b. 6

CONVECfION IN A ~TER-SATURATED POROUS BOX

Nomenclature C C° I

-

v e c t o r o f expansion c o e f f i c i e n t s o f t r i a l

solutions

component expansion c o e f f i c i e n t s o f t r i a l

solutions

d

fracture width

g

acceleration due to gravity

h

fracture separation

dimensionless horizontal box lengths hl,H 2 K permeability of rock -

k

unit vector in z direction

L

dimensional box height

LI,L 2 -

horizontal box dimensions

MA,MB h

matrices arising in Galerkin formulation

P

-

unit vector in z direction dimensionless perturbation pressure dimensional f l u i d pressure

+

q,qx,qy,qz-dimensionless perturbation velocity, components R

dimensionless Rayleigh number

R

critical

s

s p e c i f i c heat o f f l u i d

T

temperature

T o T

temperature at base o f box

C

Rayleigh number

ten~perature at surface

S

v

-

Darcy v e l o c i t y o f f l u i d

x,y,z

c a r t e s i a n c o o r d i n a t e s , z p o s i t i v e downwards

-

therma] expansion c o e f f i c i e n t

8

-

geothermal g r a d i e n t

0

-

dimensionless p e r t u r b a t i o n temperature

Oi,ek,X

and c a r t e s i a n

t r i a l temperature function in Galerkln formulation

-

thermal conductivity of saturated material

-

kinematic viscosity of f l u i d

pf -

f l u i d density

377

378

R.P. Lowell and C.-T. Shyu

Vol. 5, hb. 6

REFERENCES

I.

J. L. Beck, Phys. Fluids, 15, 1377 (1972).

2.

A. Zebib and D. R. Kassoy, Phys. Fluids, 20, 4 (1977).

3.

S. H. Davis, J. Fluid Mech. 3(3, 465 (1967).

4.

I. Catton, Jour. Heat Transfer,

5.

E. R. Lapwood, Proc. Cambridge Phil., Soc. 44, 508 (1948).

6.

R. P. Lowell, Semi-annual grant #14-08-0001-G-365,

9_.22,186 (1970).

technical

letter report no. I, USGS

10p (1977).

7.

J. M. Straus and G. Schubert, J. Geophys. Res., 82, 325 (1977).

8.

J. Bear, Dynamics of Fluids in Porous Media, Elsevier, New York, 764p (1972).