Liquid Geotherinal Systems

GEOLOGICAL SURVEY PROFESSIONAL PAPER 1044-0

Numerical Modeling of 

Liquid Geothermal Systems

By M. L. Sorey

GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

GEOLOGICAL SURVEY PROFESSIONAL PAPER 1044-D

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON: 1978
UNITED STATES DEPARTMENT OF THE INTERIOR 

CECIL D. ANDRUS, Secretary

GEOLOGICAL SURVEY 

H. William Menard, Director

Sorey, M. L.

Numerical modeling of liquid geothermal systems.

(Geohydrology of geothermal systems)

(Geological Survey Professional Paper 1044-D)

Bibliography: p. 23-25

1. Geothermal resources-Mathematical models. I. Title. II. Series. III. Series: United States.

Geological Survey. Professional Paper 1044-D 

GB1199.5S67 551.4'9 78-606183

For sale by the Superintendent of Documents, U.S. Government Printing Office

Washington, D.C. 20402

Stock Number 024-001-03112-4
CONTENTS

	Page 

Abstract. _____________________________ Dl

Introduction ____________________________ 1

Previous investigations ___________________ 2

Acknowledgments ______________________ 3

Numerical model__________________________ 3

Partial differential equations _______________ 3

Flow equation ______________________ 3

Equation of state ____________________ 4

Energy equation ____________________ 4

Numerical solution _____________________ 5

General procedure.___________________ 5

Numerical formulation of flow equation _______ 6

Numerical formulation of energy equation ______ 7

Iterative solution ____________________ 7

Coupling of flow and energy equations __________ 8

Studies of natural convection _________________ D8

Circulatory convection ___________________ 9

Box models __________________________ 9

General considerations _________________ 9

Numerical solutions __________________ 10

Temperature-dependent parameters ________ 12

Effect of density variation with pressure _____ 13

Realistic models _______________________ 14

Conducting side walls _________________ 14

Two- and three-layer models _____________ 15

Hot-spring models _______________________ 16

Summary and conclusions ___________________ 22

References ____________________________ 23

ILLUSTRATIONS

Page 

FIGURE 1. Interlacing of temperature and pressure steps _________________________________________ D6

2. Cellular convection in an extensive horizontal layer_______________________________________ 8

3. Box model for cellular convection with heating from below __________________________________ 9

4. Convective rolls in box with Jix =2hz ______________________________________________ 10

5. Velocity pattern on lower surface of cube for three-dimensional convection ________________________ 10

6. Influence of grid spacing on steady state Nu, N = number of nodes in x and y directions ________________ 10

7. Steady-state temperature distribution in box model with aspect ratio = 1.0, Ra = 100, and uniform (symmetric) initial

temperature distribution ___________________________________________________ 11

8. Steady-state temperature distribution in box model with aspect ratio = 1.0, Ra = 100, and nonuniform (asymmetric)

initial temperature distribution._______________________________________________ 11

9. The relationship between Ra and Nu determined numerically for cellular convection in box model ___________ 11

10. Nusselt number with temperature-dependent fluid properties as a function of cold-side Rayleigh number (eq. 34, dashed 

lines) and mean Rayleigh number (eq. 35, solid lines) for temperature differences T, -T0 between 40°C and 300°C. 

Critical Rayleigh numbers at Nu = 1.0. Heavy solid line is for properties fndependent of T (see fig. 9) ______ 12

11. Temperature and velocity profiles at midplane (x = hx/2) for convection cell with Ram = 100 and Tl -T0 = 100°C. 

Temperature curve A is for constant fluid properties; temperature curve B and velocity curve C are for temperaturedependent fluid properties __________________________________________________ 13

12. Box model to test effects of density variation with pressure _________________________________ 13

13. Physical model for cellular convection with conducting side walls ______________________________ 14

14. Temperature distributions in convection model with conducting side walls for w/L = 2.1 and y/L= 0.5 ________ 14

15. Isotherms for 2-layer model with Ra = 100. Isotherms for 3-layer model with Ra = 100 _________________ 15

16. Isolated cylindrical conduit hot spring model __________________________________________ 16

17. Fault plane hot spring model____________________________________________________ 17

18. Geometry of hot spring models __________________________________________________ 17

19. Relationships between dimensionless flow rates (Mc,Mp) and dimensionless temperature drop (1 <j>sp) due to conductive heat loss in cylindrical and fault plane hot spring models ______________________________ 18

20. Steady-state temperature distribution in fault plane hot-spring model with discharge = 108kg/d, W = 10 m, D = 1 Km,

andflo = 0(fc = 0) ________________________________________________ 19

21. Steady-state temperature distribution in fault plane hot-spring model with discharge = 108kg/d, W = 10 m, D = 1 Km,

andflo = 0(k = 0), andff - 10-6 cal/sec °Ccm2 at land surface _________________________ 19

22. Steady-state temperature distribution in fault plane hot-spring model with discharge = 106kg/d, W = 10 m, D - 1 Km,

andflo = 20 (solid lines), Ra = 0 (dashed in figure 20) ________________________________ 20

23. Steady-state temperature distribution in fault plane hot-spring model with discharge = 108kg/d, W = 10 m, D = 1 Km,

andflo = 200 _________________________________________________ 21
IV CONTENTS

TABLES

TABLE 1. Temperature-dependent properties of pure water ________________________________________ D12

2. Results of computer runs with density a function of pressure _________________________________ 14

3. Heat transfer results for convection in vertical channel with conducting side walls _____________________ 15

4. Effect of multiple layers on heat transfer in permeable layer with circulatory convection _________________ 15

5. Spring temperature, Tsp, total conductive heat loss, Q, and heat flux at land surface near spring, q, for selected values of

coefficient of surface heat transfer H ____________________________________________ 20

NOMENCLATURE

Symbol

C 

c

cs 

D

g 

H

Kd

Km 

k

L

Mc

Mp

Nu 

P

Q

Ra

Description

Surface area connecting nodes n

and m

Compressibility coefficient 

Specific heat at constant volume

of fluid

Specific heat at constant pressure 

of fluid

Specific heat of solid phase 

Length of fault plane conduit

Distance between nodes n and m

Distance between node m and interface with node n 

Gravitational acceleration vector 

Coefficient of surface heat transfer

Thermal conductivity of solidfluid matrix

Thermal conductivity tensor for 

solid-fluid matrix

Coefficient of thermal dispersion

Total thermal dispersion tensor 

Intrinsic permeability 

Depth from land surface to source 

reservoir in hot spring models

Dimensionless mass flow rate in 

cylindrical hot spring model

Dimensionless mass flow rate in 

fault plane hot spring model

Nusselt number (dimensionless 

heat transfer coefficient)

Fluid pressure

Lateral conductive heat loss from 

spring conduit

Rayleigh number (natural convection parameter)

Dimensions

t2/L2 

L2/t2T

L2/t2T

L2/t2T

L 

L

L

L/t2 

M/t2TL

ML/t3T 

ML/t3T

ML/t3T 

ML/t3T

L2 

L

M/Lt2 ML2/t3

Symbol Description

Rac Critical Rayleigh number for onset of convection

rc Radius of cylindrical hot spring

conduit 

T Temperature

Tb Temperature at base of hot spring 

models

Ts Temperature at land surface in

hot spring models 

t_ Time

vd Fluid flux vector or Darcian velocity

W Half-width of fault plane conduit 

a Vertical compressibility of porous

medium

ft Thermal expansivity of fluid 

7 Coefficient in second order equation of state

0 Interpolation factor in implicit/- 

explicit numerical scheme

®sp Dimensionless temperature of

hot spring at land surface 

K Fluid compressibility 

n Fluid dynamic viscosity 

p Fluid density 

pc Heat capacity of solid-fluid matrix

PO Fluid density at reference state

Ps Density of solid phase 

T Stress deviator tensor

Tn Time constant or stability limit

node n 

<f> Porosity

Dimensions

T 

T

t

L/t

L

Lt2/M

T'1

rp-2

LtYM 

M/Lt M/L3

M/Lt2T M/L3

M/L3

M/Lt2

t
CONTENTS

CONVERSION OF UNITS

Metric units are used throughout this report. Thermal parameters are reported in more familiar "working" units rather than in 

the now-standard SI (Systeme Internationale) units. The following 

table lists metric and equivalent U.S. Customary units, and "working" units and SI units for the thermal parameters.

Multiply metric units By

centimeters (cm) 

meters (m) 

kilometers (km)

Length

0.3937 inches (in) 

3.281 feet (ft) 

.6214 miles (mi)

Area

square centimeters (cm2 ) 0.155 square inches (in2 ) 

square meters (m2 ) 10.76 square feet (ft2 )

Volume 

cubic centimeters (cm3 ) 0.061 cubic inches (in3 )

Mass

gram (g) 3.528 X 10'2 ounces (oz) 

kilogram (kg) 2.205 pounds (Ib) 

Temperature: degrees Celsius to degrees Fahrenheit 

°F = 9/5°C +32

THERMAL PARAMETERS

Multiply "working " units By To obtain SI units

Thermal Conductivity

calories per centimeter 0.4187 watts per meter 

second • degree Celsius degree Kelvin 

(cal cm-1 s- 1 °C- 1 ) (Wm- 1 °K-0

Heat Transfer Coefficient

calories per square 41.87 watts per square

meter degree 

Kelvin

centimeter second 

degree Celsius 

(cal cm-2 - 1 °C- 1 )

calories per second 

(cal s-1 )

calories (cal)

calories per gram 

• degree Celsius

Heat Discharge

4.187 watts (W)

Energy 

4.187 joules (J)

Specific Heat

4.187 X 103 joules per kilogram 

• degree Kelvin

GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS

By M. L. SOREY

ABSTRACT

A mathematical model describing the physical behavior of hotwater geothermal systems is presented. The model consists of a set 

of coupled partial differential equations for heat and mass transfer 

in porous media and an equation of state relating fluid density to 

temperature and pressure. The equations are solved numerically 

using an integrated finite difference method which can treat arbitrary nodal configurations in one, two, or three dimensions.

The model is used to analyze cellular convection in permeable 

rock layers heated from below. Results for cases with constant fluid 

and rock properties are in good agreement with numerical and experimental results from other authors. Considering variations in 

fluid viscosity and thermal expansivity with temperature results in 

substantial differences in the values of the critical Rayleigh number for the onset of convection and the rate of vertical heat transfer 

compared with constant-parameter cases. For example, for a temperature difference of 300° C across a horizontal layer, the critical 

Rayleigh number based on parameters evaluated at the cold-side 

temperature is 2, whereas the corresponding constant-parameter 

value is 47T2 . Numerical simulations of more realistic models for cellular convection show that for laterally-bound reservoirs, conduction of heat across the side walls can lower the vertical heat 

transfer rate through the reservoir by 40 percent. For unbounded 

reservoirs, simulating heat conduction in impermeable layers 

above and below the convecting layer lowered the rate of vertical 

heat transfer by 50 percent.

Heat and mass transfer associated with hot-spring systems was 

analyzed for cylindrical (isolated conduit) and plane symmetric 

(fault-plane) models. Results of numerical simulations yield nondimensional relationships between spring discharge and temperature drop in the spring due to lateral conductive heat loss from the 

spring conduit. For the same spring discharge, the fault-plane 

model yields considerably greater conductive heat loss and temperature drop at springs than the cylindrical conduit model. 

Steady-state temperature distributions in the rock surrounding 

the conduit show the effects of this lateral heat flow out to distances of one conduit depth. The influence of convective motions 

in this region was investigated as a function of rock permeability. 

Finally, the time required for thermal equilibrium following hotspring development can be estimated as L2/2a, where L is the conduit depth and a is the thermal diffusivity of the saturated rock.

INTRODUCTION

Geothermal systems are receiving increasing attention as an alternative source of energy; consequently, 

there is growing interest in understanding their nature and behavior. One approach to this problem is to

attempt to deduce the physical behavior of such systems using a mathematical model. Such a model consists of a set of equations that describe the processes 

occurring within the system and the solution to these 

equations subject to conditions that prevail at a particular site.

The model approach has two important applications: the geothermal system under natural conditions 

before being disturbed by man, and the geothermal 

system during exploitation. Natural geothermal systems have been investigated by a great many workers. 

The main thrust of such studies has been to understand how geothermal systems can form and persist 

under the various geological constraints that can exist 

within the earth's crust. Mathematical modeling of 

the natural conditions of heat and fluid flow in a 

geothermal area can be used to describe certain basic 

phenomenon such as fluid convection induced by density imbalances and the associated convective heat 

transfer. Modeling also affords an economical alternative to an extensive drilling program in that various 

system properties such as reservoir volume and permeability can be inferred by using models to simulate 

known conditions and make parametric analyses.

The mathematical model can also be applied to the 

problem of evaluating the behavior of a geothermal 

system during exploitation. The main purpose, of 

course, is to estimate the quantity of recoverable 

energy and the rate at which mass and energy may be 

extracted. As information from drill holes and initial 

production characteristics becomes available, the 

model can be refined and used to simulate optimal development procedures.

Mathematical modeling of fluid and energy transfer 

in geothermal systems is a logical extension of similar 

techniques which have been successfully applied to 

the analysis of fluid flow in porous and fractured media. In geothermal areas we may need to consider the 

effects of density gradients on flow, heat transfer by 

conductive, convective, and dispersive mechanisms, 

and the simultaneous transfer of two fluid phases. BeD-l
D-2 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

cause of the complex nature of the equations and 

boundary conditions, numerical solutions are usually 

required, although analytical results are available for 

simplified cases. The development of general numerical models capable of treating all the complexities 

noted above is an area of considerable current interest. The purpose of this study is the development of a 

numerical model for single-phase fluid and heat transfer in three-dimensional porous media and the application of the model to problems involving natural 

convection and hot spring discharge in liquid geothermal systems before development.

PREVIOUS INVESTIGATIONS

No rigorous definition of what constitutes a geothermal system is really possible. In a general sense, 

we are concerned with areas of above-normal temperature and heat flow in which circulation of warm to 

hot fluids may be a significant process. Two broad 

types of geothermal systems involving water, that is 

hydrothermal systems, are recognized: liquid-dominated and vapor-dominated. The general characteristics of vapor-dominated systems in which steam is the 

continuous, pressure-controlling phase are discussed 

by White, Muffler and Truesdell (1971), Elder (1966), 

and James (1968). Most of the world's known hydrothermal systems are liquid-dominated, with water as 

the controlling pressure phase. The following discussion pertains mainly to liquid-dominated systems.

Mathematical modeling of hydrothermal systems 

has for a long time centered on the problem of convective heat transfer in a homogeneous porous layer 

heated from below. Early work by Lapwood (1948) 

and others followed the pattern of previous analyses 

of Benard convection in a layer of viscous fluid. This 

work was directed primarily at developing criteria for 

the onset of convective currents in a horizontal and 

laterally infinite permeable layer when the vertical 

temperature gradient exceeds a certain critical value. 

Beck (1972) extended the analysis of dynamic stability to a layer with lateral boundaries, that is—a box, 

and he developed criteria for the preferred shape and 

number of convection cells which form.

Many studies have been made of the increased rate 

of vertical heat transfer through a freely convecting 

horizontal porous layer. Analytical results have been 

reported by Wooding (1957, 1963), Donaldson (1962, 

1968, 1970), and others; numerical and experimental 

results are given by Hoist and Aziz (1972a, b), Combarnous and Bories (1973), and Home and O'Sullivan 

(1974). Some work has also been done on the problem 

of an inclined porous layer bounded by isothermal

surfaces as summarized by Combarnous and Bories 

(1973). They show that since the temperature gradient and gravity are no longer collineal, the fluid is constantly moving regardless of the temperature 

gradient. For the most part, previous investigations of 

the natural convection problem have treated one-layer systems with uniform boundary conditions and 

constant parameters. In the section on "Studies of 

Natural Convection," results using the numerical code 

developed in this study are compared with those previous investigations, and several more realistic models 

for natural convection are analyzed.

Another physical model used to describe hydrothermal systems is the pipe system, in which fluid is assumed to be channeled through zones of higher 

permeability in the rocks. Einarsson (1942) and Bodvarsson (1961) discuss the nature of thermal areas in 

Iceland in terms of pipe systems involving circulation 

of water to depths of 2-3 km and discharge in hot 

spring areas. Elder (1966) analyzed hydrothermal 

areas in Iceland and New Zealand using zero dimensional (lumped parameter), one-dimensional (pipe 

system), and multidimensional (homogeneous porous 

media) models to quantify the general features of heat 

and mass transfer. White (1957, 1961) used pipe systems to explain the chemical composition of waters associated with hydrothermal areas. A useful convective 

model analyzed by Donaldson (1968,1970) consists of 

a permeable channel or reservoir at the base connecting recharge and discharge columns at the sides. Geometrical and physical parameters of the model, 

including the heat available at depth, can be adjusted 

to simulate the gross features of several types of hot 

water systems. This model is analyzed numerically (in 

the section on "Studies of Natural Convection") to 

evaluate effects of conducting side walls on circulatory 

convection in the discharge channel. Additional studies of natural convection in homogeneous porous media and pipe models are reviewed by Witherspoon, 

Newman, Sorey, and Lippman (1975).

Developments of the appropriate mathematical 

equations for heat and mass transfer in single-fluid 

porous media where Darcy's law is valid have been 

made by Bear (1972), and Fernandez (1972), and Mercer (1973). Combarnous and Bories (1973) and Green 

(1963) considered the adequacy of using a single equation for energy transfer in the solid-fluid matrix; this 

equation requires assuming thermal equilibrium at 

points of contact between the two phases. They also 

analyzed the relative magnitude of dispersive and 

conductive heat fluxes as discussed in the section on 

"Numerical Model". Equations for two-phase flow in 

porous media are presented by Bear (1972), Mercer, 

Faust, and Finder (1974), and Witherspoon, Newman,
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-3

Sorey, and Lippman (1975). The common approach is 

to assume Darcian expressions for the flow of each 

phase and to combine them into a single flow equation 

under the assumption of equal pressure in both 

phases. The corresponding energy equation accounts 

for the effects of phase change, and the two-phase 

model can be applied to steam-water geothermal systems. Numerical solution of the coupled partial differential equations for mass and energy transfer by 

integral (Galerkin), finite element techniques is discussed by Mercer (1973) and Faust and Mercer (1975). 

Finite difference schemes were utilized by Hoist and 

Aziz (1972a, b) and Home and O'Sullivan (1974). A 

particular difficulty in numerical solution of transport 

equations is the generation of numerical dispersion error which can produce instabilities and spreading or 

smearing of temperature fronts. Lantz (1971) discusses the sources of these errors and methods for 

minimizing their effects.

Applications of numerical models to geothermal 

reservoir simulation are rather limited to date. Mercer 

(1973) attempted to simulate the Wairakei, New Zealand, thermal system under both steady state and exploitation conditions. At this stage his results are 

limited by the possible inadequacy of a two-dimensional areal model to simulate conditions in this system and the initiation of two-phase flow in the 

reservoir following several years of fluid production. 

Solutions to hypothetical two-phase reservoir problems are discussed by Lasseter, Witherspoon, and 

Lippman (1975) and Faust and Mercer (1975). Numerical models for steam stimulation in petroleum 

reservoirs have utilized similar equations and solution 

procedures as those applied to the geothermal problem (Weinstein and others, 1974; Coats and others, 

1973).

The equations and numerical scheme used in this 

study are described in the section on "Numerical 

model". An integrated finite difference method is 

used, and solutions to the heat and fluid flow equations are interlaced in time. The model is applied to 

the problem of circulatory convection in permeable 

layers to demonstrate its capabilities and to extend 

the analysis to more realistic geometries and boundary conditions. Applications of numerical modeling to 

the transfer processes associated with hot spring systems are also presented.

ACKNOWLEDGMENTS

An earlier version of this paper constituted a doctoral dissertation in the Department of Civil Engineering, University of California, Berkeley. The study 

was done under the general direction of Paul A. Witherspoon, Jr., and Ralph Greif of the University of California and Arthur H. Lachenbruch of the U.S. 

Geological Survey. I am also indebted to I. G. Donaldson of the Department of Scientific and Industrial Research, Lower Hutt, New Zealand, and to T. J. 

Lasseter, Marcello Lippman, and T. N. Narasimhan, 

colleagues at the University of California, for valuable 

discussions and suggestions. T. J. Lasseter also 

assisted in developing the numerical procedures used 

in this study.

NUMERICAL MODEL

Macroscopic description of mass and energy transport in porous media has been developed from the 

microscopic level using principles of statistical mechanics. Bear (1972) and Fernandez (1972) present 

detailed discussions of the assumptions involved and 

the resulting equations. In this study we utilize the results of these authors, making certain simplifications 

to yield a set of two partial differential equations for 

the flow of water and heat and an equation of state 

relating fluid density to temperature and pressure.

PARTIAL DIFFERENTIAL EQUATIONS

FLOW EQUATION

To describe mathematically the flow processes in 

porous media, we must consider various properties 

such as pressure and temperature as averages over 

volume elements which are large compared with dimensions of individual pores and solid particles, but 

small compared with the dimensions of the exterior 

boundaries of the medium. This introduces the 

assumption of what Bear (1972, p. 19) terms a representative elementary volume or REV and allows the 

actual porous medium to be treated as a continuum in 

which average values of dependent variables can be 

assigned to mathematical points. Thus we can write 

an equation for the conservation of mass as

(1) dt

where p = fluid density, 0 = porosity, t = time, and 

urf= fluid flux vector. The actual fluid velocity, averaged over the appropriate REV, would be equal to 

urf/0. The momentum equation for flow in porous media is

VP+V • T-pg = 0 (2)

where P — fluid pressure, T = stress deviator tensor, 

and g = gravitational acceleration vector. Equation 2 

is obtained from the momentum equation for a homo-
D-4 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

geneous fluid system (as given by Bird and others, 1960, p. 79) by neglecting inertial forces

which are small compared with viscous forces (V • T ) 

for laminar flow in porous media. For isotropic porous 

media, the empirical observations of Darcy (1856) and 

the similarity with Poiseuille-type flow in small-diameter tubes leads to the relation

V-7 =-£-i (3)

where n = fluid viscosity and k = intrinsic permeability. Combining equations 2 and 3 yields

(4)

or Darcy's law for laminar flow in isotropic porous media. For anisotropic media, Bear (1972, p. 105) and 

Raats and Klute (1968, p. 540) developed a similar expression for Vd in terms of a permeability tensor k. In 

this development, k will be treated as a scalar which 

can vary spacially, because in the geothermal systems 

we wish to model, little is known regarding the tensorial properties of JJ, and it is likely that large-scale inhomogeneities related to fractures and channeling will 

dominate the convective system.

Combining the continuity and momentum equations yields the flow equation

(5)

For slightly compressible fluids such as water, we wish 

to solve the flow equation in terms of pressure rather 

than of density. Hence, we expand the right-hand side 

of ( 5 ) as

d</>p _ 

dt = p

dP

dP

dP 1 dt

dp

__

dP

dP dp dT_ 

dt J

(6)

where x = fluid compressibility, 0 = fluid thermal expansivity, and a — vertical compressibility of the porous medium. Following Remson, Hornberger, and 

Molz (1971) we assume that the velocity of solid particles is negligible and deformation of the solid matrix is 

significant only in the vertical direction. For the problems considered in this study, the thermal expansivity 

term </>p/3 dT/dt can be shown to be small compared 

with the compressibility terms (pa + </>PK) dP/dt owing 

to the large difference in response times between the 

thermal and hydraulic systems, that is dT/dt « 

dP/dt Thus, we will neglect the expansivity term in

equation 6 and write the final form of our flow equation as

V- I+/D — 

M

dP

dt (7)

where C = p (</>K + a) = compressibility coefficient.

EQUATION OF STATE

We wish to consider density as a function of temperature and pressure. Following Fernandez (1972, p. 

54) we can write

(P-P) (8)

where p0 = fluid density at some reference state 

(T0 P0 ). Expanding the exponentials in series form and 

retaining terms of first order yields

p = p0 [l-0Cr-T0 )+K(P-P0 )]- (9)

For geothermal applications where (T-T0 ) may be 

large, a more accurate form for the temperature dependence was suggested by Wooding (1957, p. 274) as 

follows:

W (10)

Evaluating (p0 ,T0 ) at 25° C and setting 0 = 3.17 X 10'4 

°C'1 and 7 = 2.56 x 10" 6 °C-2 , yields p-p0) to an accuracy of 2 percent in the range of 25°C to 300° C. Between 4°C and 25° C a first order expression can be 

used.

In comparing the magnitude of the effects of temperature and pressure on the calculation of fluid density, it is seen that, for liquid geothermal systems, the 

pressure dependence can usually be neglected. For example, assuming (P-P0 ) = 100 kg/cm2 (equivalent to 

about 1,000 meters of head), K = 4.5 X 10'5 cm2/kg 

gives (P-PO )PO = K (P-P0 ) = 0.0045. However, for 

(T-T0 ) = 100°C and 0 = 5 X 10'4 °C'1 , (p-pfl)/P0 = 

0 (T0 -T )= 0.05, an order of magnitude greater. Thus, 

in the numerical results that follow, equation 10 was 

utilized to compute fluid density except as discussed 

in the section on "Box Models". It should be noted, 

however, that in solving the flow equation 7 a portion 

of the density change with time is always derived from 

the fluid compressibility term (frpKdP/dt. In effect, 

then, density changes with pressure contribute to 

transient pressure changes, but calculated values of 

density are not updated as pressures change, only as 

temperature changes.

ENERGY EQUATION

The energy equation in porous media can be written
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-5

following Mercer (1973, p. 7) as

-0) pscs] -pvdc-VT (11)

where ps = density of solid phase, cs = specific heat of 

solid phase, Km = total thermal dispersion tensor, and 

c = specific heat at constant volume of fluid phase. 

This relationship is obtained by combining the corresponding energy equations for the solid and fluid 

phases, as developed by Bird, Stewart, Lightfoot and 

others (1960, p. 314), with terms accounting for energy 

increases by viscous dissipation and compression neglected. To obtain (11) several assumptions are made. 

First, we assume that the fluid and solid are in thermal equilibrium at each point of contact so that a single temperature, T, can be assigned to the saturated medium in a REV. This assumption has been 

shown to be valid for most porous media problems 

(Combarnous and Bories, 1973; Hoist, 1970; and 

Green, 1963). However, a single energy equation may 

not be adequate under certain circumstances, such as 

forced convection at high flow rates through porous 

material with relatively low thermal conductivity 

(Curry, 1970; Combarnous and Bories, 1973). One can 

also envision difficulties in treating ground-water flow 

through fractured rock as porous media flow necessitating the use of an effective REV which is not small 

compared with the overall size of the system.

The flux of heat represented byXm-VT is a composite of heat conduction through the solid-fluid matrix 

and the spreading effects of dispersion in the fluid 

phase. Bear (1972, p.650) gives the components of the 

total thermal dispersion tensor as

Km =Kc+<j>Kd (12)

where Kc = thermal conductivity tensor for solidfluid matrix and K^ = coefficient of thermal dispersion. Bear states that heat transfer by dispersion is 

analogous to mass transfer by hydrodynamic dispersion and results from the distribution of local velocities caused by the presence of grains and molecular 

diffusion. Mercer (1973, p. 10) develops the functional 

form of the thermal dispersion tensor as

(13)

where e^i = dispersivity of the medium, a fourth rank 

tensor, and vkt vl = components of velocity in the k 

and / directions, and v = magnitude of the velocity. 

Green (1963) has shown analytically and experimentally that for relatively uniform porous media, heat 

transfer due to dispersion is small compared with conduction for values of velocity normally encountered in 

ground water systems. Unfortunately, it is difficult to 

apply these laboratory results to field situations where 

considerable variability in permeability may increase 

the dispersive effect. Indeed, it is possible that where 

dispersion is important in heterogeneous systems, the 

assumption of thermal equilibrium between phases 

may also break down and a more complicated mathematical description for heat transfer will be required. 

With this in mind, and considering the absence of information as to the magnitudes of the tensorial components of (12), it seems reasonable to treat Km as a 

scalar Km, which could represent an enhanced conductivity to account for dispersion. For example, Henry and Hilleke (1972, p. 33) used a value of Km which 

was five times the static thermal conductivity to 

match the temperature distribution in a coastal salt 

water-fresh water aquifer. Thus, the energy equation 

11 becomes

= V'KmVT-pvdc-VT (14)

where ~pc = [$pc + (1-0) pscs] = heat capacity of 

solid-fluid matrix.

NUMERICAL SOLUTION

The flow equation 7 is linear only if we neglect the 

pressure dependence of p and n; however, both p and \i 

must be treated as functions of temperature, and thus 

the pressure and velocity fields depend on temperature. The energy equation 14 is nonlinear because of p 

(T) and c (T), and the temperature field depends on 

the velocity field. Thus, we must solve a coupled nonlinear set of equations. Further, in the geothermal systems we wish to study, properties such as permeability 

and thermal conductivity will vary spacially and the 

boundary geometry may be complex. For these reasons numerical solutions are required, and a code 

which can treat multidimensional problems with variable geometries would be desirable. A computer program for slightly compressible heat and fluid flow was 

therefore developed.

GENERAL PROCEDURE

The basic procedures in the program were derived 

from the work of Edwards (1969) and Lasseter (1974). 

Before writing finite difference approximations, we 

integrate the equations over the volume of a finite element and apply the divergence theorem to reduce the 

flux terms to integrals over the area of the element. 

Thus, for an equation of the form

at (15)
D-6 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

integration yields,

vol v " area * ' 

To solve this equation numerically, we discretize the 

continuum into a finite series of elements, each with a 

specified volume and a surface area between it and 

surrounding elements. We must then assume that the 

volume-integrated factors on the left-hand side are 

constant over each element, and the surface-integrated terms are reasonably constant over each of the connecting areas.

Consider an element "n" and two adjacent elements 

ml and m2 as shown below

m\,n n,m\ mi,n

Writing equation 16 in difference form for this system 

yields

AXn

= " " lVAM>mi

+ Ar

Voln an "£* = An,mi

m

vn,rm rnjm 

>n,m (17)

To evaluate the flux at each interface, we assume that 

the average values of 0 and X over the element apply at 

a particular point or node within the element. The 

gradient (VX)n,m is approximated as

(18)

In general the distances dn>m and dm>n are measured 

normal to the interface An>m.

To preserve continuity of flux and potential across 

the interface, the coefficient (3n>m must be evaluated as 

the harmonic mean of the values /8m for node m and @n 

for node n. Thus

_ Pn,m = a Pm^n.m

(19)

where Dn,m = dn,m + dm>n- Equation 17 can now be 

written

Voln an AXn "'* (20)

The system of algebraic equations generated by 

writing equation 20 for each node in the system is 

solved by iteration. This general procedure puts no inherent restrictions on the number of nodes connected 

to node "n" or on the shape of the connecting surfaces. 

Thus it can be applied, as noted previously, to multidimensional problems with simple or complex element geometries. Narasimhan (1975) refers to this 

scheme as the IFD (Integrated Finite Difference 

Method) and gives a detailed conceptual comparison 

of the IFD method with the finite element method.

NUMERICAL FORMULATION OF FLOW 

EQUATION

The flow equation is

(7)

After volume integration, the difference equation can 

be written

AD

Cn Voln-^

-pkA_\ (Pm-Pn)

M >*>m Dm>n (21)

,( P2kAr,g\

T\ I

/n,m

where i\n>m is the direction cosine of the angle between 

the outward normal of "n" towards "m" and the gravitational acceleration vector g. Evaluating (k/n) n>m as 

the harmonic mean

(A) = V „ Jn,m

M

(22)

and weighting the density by

i

i.n Pn

Pn,m ~ D

(23)

n,m

we can rewrite equation 21 as 

C.VOLA5L, l/n m (Pm — P ) + w tn^\

( k\ where *n m = transductance = [- pn m \ ~~) n,m M

A ( k \

r\"m ] and U> n, m = \.PZn,m \ ~ )n,m An>m T\n,m #]• 

un,m M

We wish to solve the flow equation implicitly. Thus, 

we expand equation 24 into implicit and explicit parts 

as
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-7

VoL - Pn)

(25)

where P is the pressure at the old time and AP is the 

change in P during the time step A£. The interpolation 

factor 0 can be set to 1 for a fully implicit or backward 

differencing scheme, ¥2 for a Crank-Nicholson implicit scheme, or 0 for a fully explicit scheme. In solving 

the flow equation, 0 is set to 1. The explicit part of 

equation 25 is computed first and then the implicit 

part is solved by iteration.

NUMERICAL FORMULATION OF 

ENERGY EQUATION

The energy equation 14 is 

— dT (14)

Before integrating, we rewrite this equation in divergence form as

- V-pvdcT + T V-pvdc • (26)

The corresponding difference equation is

(PvdcA) n>m (Tn!m -Tn)]

where velocity is defined as positive into node n and 

Kn>m is the harmonic mean conductivity.

Expanding into implicit and explicit terms yields

n,m

(pvdcA) n>m (Tn>m -Tn)]

(28)

un,m

Temperature-dependent material properties are evaluated at some intermediate time between t and t + A£. 

The procedure by which this point is determined is 

discussed by Lasseter (1974, App. C). The value of 0 

used to solve the energy equation is adjusted between 

0.57 and 1.0 in order to damp out small oscillations in 

the solution.

The temperature at the interface Tn>m could be 

evaluated by a linear interpolation

n,m (29)

n,m

It can be shown, however, and anticipated physically, 

that the temperature variation in the presence of convection is not linear and that the interface temperature is weighted in favor of the upstream node 

temperature. Numerically, it is also found that instabilities develop unless we use this upstream weighting. 

Thus, we rewrite equation 29 in terms of the adjusted

distances d'n>m and d'm>n as

T 1 n,m =

" -L

D

(30)

n,m

The energy equation now becomes

(Vol - PvdcA)n>mTn]

(31)

in terms of transductances

,m

and

un,

n,m

dm,

m,n

The instability noted above may develop if W n , m 

and Sl2 n ,m have the same sign. To avoid this, we 

weight the interface temperature in favor of the upstream node temperature as follows. If (pud) n ,m> the 

mass flow rate from node "m"'to "n", is large and positive, the transductances may both be positive. In this 

case the weighting d'm ^ n /D m ^ n on Tm is reduced until 

Q 2n m is just negative. If (pvd} n>m is large and negative, 

the transductances will both be negative and the 

weighting d'n>m/Dn<m on Tm is reduced until Wn>m is 

positive. This computation is done internally by the 

computer program so that stability is always assured. 

The accuracy of this weighting technique was verified 

by Lasseter (1974) by comparing numerical results 

with analytical solutions for forced convection 

through porous media.

ITERATIVE SOLUTION

For each time step, the flow equation 25 is solved in 

a fully implicit fashion with the interpolation factor 0 

set to 1. With the energy equation 31, 0 is adjusted 

between 0.57 and 1.0 and an implicit-explicit optimatization method is used. In the latter, a stability 

limit or time constant r is defined for each node as
D-8 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

Tn = (32)

where tt* n ,m signifies that we sum the conductive 

components of Q rtim over all connections and take the 

difference between inflow and outflow of the convective components. rn represents the largest stable time 

step that could be used to solve equation 31 explicitly 

for node n. Physically, it is the approximate time required for a node to react significantly to changes in 

temperature of nodes it is connected to. Thus, for each 

time step, At, the energy equation is solved explicitly 

for nodes with rn <At and implicitly for all other 

nodes.

The time step used to solve the energy equation is 

controlled mainly by an input parameter, TVARY, 

which is the desired maximum temperature change 

per time step. For accuracy in transient problems, 

TVARY is set to about 1 percent of the maximum 

temperature change expected in the system. For 

steady state solutions, TVARY can be increased to 

speed convergence. Control of the time step for the 

flow equation is discussed in the next section.

Within each time step, the set of equations 31 for 

each node is solved by using the iterative technique 

developed by Edwards (1969) and discussed by Lasseter (1974). Convergence is generally rapid, but the 

number of iterations necessary depends on the relative number and time constant of interconnected 

nodes, the relative transductances for these connections, and TVARY which is used in the convergence 

test. The program adjusts At to optimize on about 40 

iterations per time step.

COUPLING OF FLOW AND ENERGY EQUATIONS

Simultaneous solution of the flow and energy equations can be accomplished numerically by interlacing 

the solutions in time. Starting with initial temperature and pressure distributions at time t0 , we assume

that temperatures are reasonably constant over a 

short interval of time and solve for a new pressure distribution at time t 0 + At0/2. We then use this pressure 

distribution to determine the new temperature distribution at t0 + At. These temperatures are then used to 

compute pressures at t 0 -1-3At/2, and so on. This technique is illustrated in figure 1.

The actual scheme used is simplified due to the 

large difference between time constants of the two 

equations. The flow equation has time constants 

which are typically on the order of seconds while the 

energy equation has time constants on the order of 

tens of days. Thus after each temperature time step, 

the flow field will rapidly converge to a quasi-steady 

state solution. This is not a true steady state because 

temperatures are varying slowly with time and pressures will slowly but continually change. But we can 

solve for the quasi-steady state pressure distribution 

once we have computed a new temperature distribution without worrying about finding P and v<j at At/2. 

The difficulty with this technique is that many flow 

cycles per energy cycle are required to properly solve 

for the new pressure distribution because of the difference in time constants. A cycle is a step forward in 

time, whether it is a flow calculation or an energy calculation. For each temperature cycle, about 10 flow 

cycles are needed to compute the quasi-steady state 

solution. The time step for each flow cycle is adjusted 

to utilize about 40 iterations per cycle.

It should be noted that the equations presented and 

the problems solved in this paper do not involve 

source terms. With a mass source, such as a recharge 

or discharge well, the flow field would not necessarily 

equilibrate within each thermal cycle and the interlacing scheme shown in figure 1 would have to be followed.

STUDIES OF NATURAL CONVECTION

It is useful to identify two types of convective heat 

transfer in porous media. Forced convection results

SEQUENCE 

NUMBER 1

PRESSURE TIME STEPS 

3

t o +Ar/2 t0 +3A?/2 to +2A?

TEMPERATURE TIME STEPS 

FIGURE 1.—Interlacing of temperature and pressure steps.
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-9

from externally imposed hydraulic conditions; fluid 

velocity is independent of the temperature distribution. Natural or free convection arises from buoyancy 

effects caused by temperature differences. In geothermal systems, both types of flow occur and convection could be termed mixed. One example of this 

mixed convection is the circulation of hydrothermal 

fluids from recharge areas to hot spring discharge 

areas. The driving force for this circulation, which 

may extend to depths of several kilometers, could be 

derived from a combination of the difference in elevation of the recharge and discharge areas and the density difference between the cold downflowing water 

and hot upflowing water.

CIRCULATORY CONVECTION

Of particular interest in geothermal studies is the 

phenomenon of circular or cellular convection. Consider a laterally extensive layer of permeable material 

heated from below as in figure 2. For the idealized 

condition where the upper and lower surfaces of the 

permeable layer are at constant temperatures TQ and 

T}, the conditions for the existence of the convective 

pattern shown are related to the dimensionless Rayleigh number, defined as in Witherspoon, Newman, 

Sorey, and Lippman, (1975) as

_ figk (T, - T0 )L 

Ka ~ (Km/Poc)(»/p0 )

(33)

where L = thickness of the permeable layer, and the 

other terms are as previously defined. The Rayleigh 

number relates buoyancy forces to viscous forces; and 

for the system in figure 2, Ra must exceed a critical 

value of 4ir2 for convection to occur (Lapwood, 1948). 

For permeable layers with lateral boundaries, analytPermeable layer

A A A A A 

FIGURE 2.—Cellular convection in an extensive horizontal layer.

ical results by Sutton (1969) and Beck (1972) can be 

used to determine the critical Rayleigh number (Rac) 

and the preferred convective mode as discussed below. 

The preferred convective mode refers to the number 

and shape of convective cells which would form at values of the Rayleigh number above Rac.

In real systems, it is doubtful that a critical Rayleigh number exists because there would be horizontal 

temperature variations along the upper and lower 

bounding surfaces. Convection is then set up for any 

value of Ra>Q, provided the adiabatic temperature 

gradient dT/dz = gT(3/cp is exceeded (Landau and Lifshitz, 1959, p. 9). However, we would expect convective effects on the thermal and hydrologic regimes to 

be negligible except in geothermal areas with above 

normal heat flow and temperature gradients. Donaldson (1968) estimates that in geothermal areas the 

Rayleigh number is in the range of 500 to 5,000; and 

Home and O'Sullivan (1974) estimate that in the 

Wairakei, New Zealand geothermal reservoir Ra is 

near 5,000.

Previous investigations of cellular convection in porous media have treated either the infinite layer case 

or box models with insulated, impermeable vertical 

sides. Results from box model studies have application to both laterally extensive aquifers where the 

zone of influence of each of the convection cells forms 

a box with respect to adjacent cells and to structurally 

bounded aquifers or channels. The effects of assumptions of isothermal surfaces above and below the 

aquifer and insulated side walls are discussed in the 

section on "Realistic Models." In the next section, a 

general description of heat and fluid flow in box model 

systems will be given and numerical results from the 

program developed for this study will be compared 

with analytical, experimental, and numerical results 

from other authors.

BOX MODELS

GENERAL CONSIDERATIONS

The geometry and boundary conditions for the box 

model are shown in figure 3. All sides are impermeable 

and the vertical sides are insulated. Following the 

analysis by Beck (1972), the problem of the onset of 

convection is solved by writing steady state forms of 

equations 7 and 14 for static conditions (no convection) and for a small disturbance upon the static state, 

under assumptions of constant coefficients and constant fluid density except in the gravity term of 7. The 

two sets of equations are subtracted and the resulting 

set of equations nondimensionalized. Solutions have 

been obtained by the linear and energy methods, both 

yielding the same result for the box model problem
D-10 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

FIGURE 3.—Box model for cellular convection with 

heating from below.

(Beck, 1972). The critical Rayleigh number, which is 

the smallest eigenvalue of this eigenvalue problem, 

depends entirely on the aspect ratios hjhz and hy/hz. 

The minimum value of Rac is 4ir2 for the infinite layer 

case, and bounding the fluid tends to make it more 

stable, that is Rac >47r2 .

For values of Ra just above the critical value, several 

modes of convection are possible. The motion can be 

either 2- or 3-dimensional with single or multiple cells. 

In contrast with the corresponding problem of free 

convection of an enclosed homogeneous fluid, 2-dimensional motion is possible because a no-slip boundary condition at the walls is not required for the 

porous media system. A roll is defined as a cell with 

only two non-zero velocity components. It has the appearance of a cylinder with axis parallel to the zero 

velocity coordinate. Beck's (1972) solution to the 

minimum eigenvalue problem shows that rolls are 

preferred over 3-dimensional cells whenever the 

height hz is not the smallest dimension. When rolls do 

form, they are usually parallel to the shorter side, although the overriding rule is for the number of rolls 

and the direction of their axes to be such that each roll 

has the closest approximation to a square cross section 

as possible. For example, for hx = 2.0, hy = 0.75, and hz 

= 1.0, two rolls form parallel to the y-axis as shown in 

figure 4. Three-dimensional cells occur when hx , hy, 

and hz are nearly the same size, that is, a cube, and 

when the height hz is less than both lateral dimensions. For a cube, the motion resembles a toroid with 

vertical axis through the center of the box. The velocity pattern in the (x,y) plane is shown in figure 5.

One effect of this convective motion is to increase 

the heat transferred vertically through the system. 

The Nusselt number Nu is defined as the ratio of heat 

flow with convection to heat flow by conduction in the 

absence of convection. For Ra<Rac, Nu =1. In regard

Y

hx

FIGURE 4.—Convective rolls in box with hx = 2 hz. 

Z Y

X

FIGURE 5.—Velocity pattern on lower surface of 

cube for three-dimensional convection.

to the number and shape of the convection cells which 

form at Rayleigh numbers above Rac, the analysis of 

Platzman (1965) indicates that the solution exhibiting 

the maximum heat transfer is physically preferred. 

Thus, for a square box with Ra = 100, a single roll for 

which Nu = 2.60 would be preferred over a two-roll 

mode for which Nu = 2.20. Numerical and experimental studies by Combarnous and Bories (1973) and 

Hoist and Aziz (1972b) show that Nusselt numbers 

are not significantly different for two-and three-dimensional motions involving the same number of convection cells. Using the program code, comparison of 

plane symmetric, axisymmetric, and three-dimensional models with one convection cell at Ra = 100 

yielded differences in Nu at steady state of less than 

10 percent. In the discussions which follow, only twodimensional convective motions are considered.

NUMERICAL SOLUTIONS

To model the box problem using the program, a 10 

X 10 grid of nodes was normally used for 2-dimensional solutions and a 6 X 6 X 6 grid for 3-dimensional 

solutions. The effect of grid spacing on the Nusselt 

number for a typical 2-dimensional problem is illustrated in figure 6 where the error for the 10X10 grid is 

only 3 percent of the extrapolated value for A jc = Ay 

= 0.

For most problems, only the steady state solution 

was desired. Hence, starting with initial temperature
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS

Nu

2.6

2.4

2.2

2.0

.05 .10

-» 1/N

FIGURE 6.—Influence of grid spacing on steady state Nu, N = 

number of nodes in x and y directions.

and pressure conditions, the thermal time step was 

advanced rapidly using a large value of TVARY as discussed in the section on "Coupling of the Flow and 

Energy Equations"; steady state conditions were indicated when the hot and cold side Nusselt numbers became equal.

For a given aspect ratio, the choice of initial conditions affected the cellular pattern and the Nusselt 

number as steady state was approached. For example, 

using uniform initial temperature and pressure distributions in a square with Ra = 100 results in a 2-cell 

convective pattern with Nu = 2.20 as indicated in figure 7. If the initial temperature distribution is not 

symmetric, unicellular convection develops with Nu = 

2.60 as shown in figure 8. In addition, Home and 

O'Sullivan (1974) report that for a uniform initial 

temperature distribution in a square box, heating the 

lower boundary slowly instead of instantaneously, results in unicellular rather than multicellular motion. 

Thus, the convective pattern which may exist in a 

geothermal system probably depends, in part, on its 

thermal history.

The relationship between Ra and Nu determined 

numerically using the program is plotted in figure 9. 

These results are for a square, and the Nusselt number for each value of Ra corresponds to the cellular 

pattern which maximizes the heat transfer, which is 

the preferred cellular mode. As found numerically and 

experimentally by Combarnous and Bories (1973), the 

preferred cell size is a decreasing function of Rayleigh 

number, varying from an aspect ratio (width to 

height) of 1 at low values of Ra to an aspect ratio of 

0.25 at Ra = 1,000. Also shown in figure 9 are numerical results from Combarnous and Bories (1973), obtained from finite difference solutions of nondimensional forms of equations 7 and 14. These authors also present experimental data for several combinations of porous materials and fluids. The Ra-Nu 

relationship in figure 9 fits with the average experimental results; variations between the different experimental systems were attributed to finite heat 

transfer coefficients at solid-fluid contacts. However,

210 C

FIGURE 7.—Steady-state temperature distribution obtained 

from uniform (symmetric) initial temperature distribution.

FIGURE 8.—Steady-state temperature distribution obtained 

from nonuniform (asymmetric) initial temperature 

distribution.

for materials with thermal conductivity Km between 2 

and 5 mcal/(s, °C cm) (for example, natural earth materials), the numerical model assuming an infinite 

heat transfer coefficient was found to be adequate.
D-12 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

10

9

8

7

QC

LU 5 

00

LU 

CO 

CO

D

Combarnous and 

Bories(1973)

400

RAYLEIGH NUMBER

1000 2000

FIGURE 9.—The relationship between Ra and Nu determined 

numerically for cellular convection in box model.

The existence of oscillatory or fluctuating convective states has been investigated by Home and O'Sullivan (1974) and Caltagirone, Cloupean, and 

Combarnous (1971) for box models in which the lower 

boundary is uniformly and nonuniformly heated. 

Their work includes experimental and numerical observations of regular and irregular variations in the 

convection cell pattern and the resultant Nusselt 

number at Rayleigh numbers above about 300. Real 

time periods for these fluctuations at Ra near 1,000 

were approximately equal to the expression (0.006 

L2^/Km), which yields a range from 200 to 700 years 

for typical thermal properties in systems with 

L near 1 Km2 . Corresponding variations in Nusselt 

number of ± 15 percent were found. Numerical simulations using the program for the uniformly heated 

case indicate similar instabilities at a Rayleigh number of 1,000 for motion which is dominantly unicellular, although the Nusselt number variation was only 

±6 percent. Thus, from the standpoint of modeling 

natural convection in geothermal systems, the rather 

large periodicity and limited quantitative effect on 

heat flow associated with this oscillatory phenomenon 

suggest that its practical significance may be limited.

TEMPERATURE-DEPENDENT PARAMETERS

In previous studies of the critical Rayleigh number 

and the Ra-Nu relationship, the properties /u and c 

were treated as constants and the first-order equation 

of state (9) was used to relate fluid density to temperature. In the program the relationships n (T) and c (T) 

can be tabulated or suitable analytical expressions can

be solved at each thermal time step. Table 1 lists values of p, c, and p for pure water at temperatures covering the range found in hydrothermal systems. The 

viscosity variation with temperature is considerably 

larger than that for density or specific heat, and the 

variation suggests that results obtained with a constant viscosity assumption may be significantly in error. Actually, thermal expansivity 0, which appears in 

the equation (33) for Rayleigh number, increases considerably with temperature and further complicates 

the analysis.

TABLE 1. Temperature-dependent properties of pure water

[From Dorsey (1968)]

P

(atm)

P 

(g/cm3 ) (g/m-s) (cal/°C g)

10-----

60-----

160 -----

260-----

Q1 ft

1

1

50

50

inn

1.000

.983

.910

.788

692

1.307

.467

.174

.109

.090

1.000

.951

.833

.763

.741

The combined effects of the temperature dependence of p, c, IJL, and /3 on the analysis of cellular convection were evaluated for a box model with aspect 

ratio of 1.0. The Nu-Ra relationship was examined for 

several values of the vertical temperature difference 

(7VT,,), using an average temperature in each case of 

160°C. The second-order relation (eq. 10) for p(T was 

used and permeability k was adjusted to keep Ra constant as (T^-T,,) varied. For each run, two Rayleigh 

numbers were calculated. One could be termed the 

cold-side Rayleigh number Ra0 , 

where

Ran =

gkL(T-TQ ]

Km

(34)

and the subscript (0) indicates evaluation at T = T0 . A 

mean Rayleigh number Ram was also computed as

gkL(T-T0 )

Km

, 

L (35)

where the subscript (m) indicates evaluation at T = 

(TrT0 )/ 2 = 160°C, and ft, is an effective expansivity 

computed as

(36) * ~ Po (T-T0 )

For each (Tj—T0 ), eight-point tabulations of n(T) and 

c(T) were used with linear interpolation between 

points.

The results are shown in figure 10. The heavy solid 

curve is the same as the curve in figure 9 for the constant parameter case. We note that for a given system, 

computing the Rayleigh number based on n, c, and #
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-13

1

10 20 40 60 100 1000

RAYLEIGH NUMBER

FIGURE 10. — Nusselt number with temperature-dependent fluid 

properties as a function of cold-side Rayleigh number (eq. 34, 

dashed lines) and mean Rayleigh number (eq. 35, solid lines) 

for temperature differences Tj-To between 40°C and 300°C. 

Critical Rayleigh numbers at Nu = 1.0. Heavy solid line is for 

properties independent of T (see fig. 9).

evaluated at T0 yields a value of Nu which increases 

with (TV—T0 ) and is greater than the constant parameter result for (Tt—T0) greater than about 40°C. In 

contrast, the Nusselt number corresponding to Ram is 

less than the constant parameter result and decreases 

with(Tt— T0 ).The position of the Nu-Ram curves may 

vary with the actual average temperature (Tj—T0 )/2 

due to nonlinearities in the n(T), c(T), p(T) relations, 

but the choice of (T,-T,)I1 = 160° C should be representative of geothermal reservoirs. The explanation 

for the lower heat transfer with variable fluid properties (Nu-Ram) may also be related to the nonlinearity 

in /i(T) which causes the centroid value of n for T0 < T 

< TI to be greater than n at (Tt—T0 )/2. Because Ra is 

inversely proportional to /*, the effective Ra should be 

less than Ram. A similar trend of decreasing heat 

transfer with increasing (Tt —T0 ) was reported by 

Hoist and Aziz (1972a) with n-heptane as the convecting fluid.

Results in figure 10 also show that the critical Rayleigh number is influenced by variable fluid properties. In terms of the cold-side Rayleigh number Ra0 

the onset of convection occurs near the theoretical value of 4?r2 for (Tj-To) < about 40°C but occurs at values of Ra < 4?r2 for (T,-T0 ) > 40°C. For (T,-T0 ) = 

300° C, the critical value of Ra0 is only 1.7. However, 

the critical value in terms of Ram is always greater 

than 47r2 , reaching a value of 61 for a temperature difference of 300° C. Finally it should be noted that variable fluid properties also cause the temperature and

-15

DARCIAN VELOCITY (m/d x 10+3 ) 

-10 -5 0 +5 +10 +15

130 150 170 190 

TEMPERATURE, IN DEGREES CELSIUS

210

FIGURE 11.—Temperature and velocity profiles at midplane (x 

= hx/2) for convection cell with Ram = 100 and T^-T^ = 

100°C. Temperature curve A is for constant fluid properties; temperature curve B and velocity curve C are for temperature-dependent fluid properties.

flow patterns to be asymmetric as shown in figure 11. 

The asymmetry reflects the fact that the hotter, less 

viscous fluid near the bottom moves more easily than 

the colder fluid above.

EFFECT OF DENSITY VARIATION WITH PRESSURE

The effects of allowing fluid density to vary with 

pressure according to equation 7 were evaluated for 

the system shown in figure 12. Two runs were made, 

one with p (T) and one with p = p (T,P), forCTj-T,,) = 

1°C and 100°C. For a temperature difference of 100°C 

across the box, the effect of temperature on density 

should dominate, as discussed in the section on "Partial Differential Equations." However, for a temperature difference of 1°C, the density increase from top 

to bottom due to the increase in hydrostatic pressure 

is 10 times the density decrease due to the temperature difference.

The results of the computer runs are listed in table 

2. Although the Nusselt number for run 4 with p = p 

(T,P) is slightly higher than the value for run 3 with p 

= p (T), the difference in heat transfer between the 

various runs were negligible. Similarly, temperature 

and velocity distributions were essentially the same 

for each run. Thus, for problems or circulatory convection in porous media, density can probably be considered a function of temperature only.
D-14 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

-1km

Tkrn

Tor, -r0 =i,ioo°c

j8=5xlO~40 C~! 

=5xl()- 5 cmVkg

FIGURE 12.—Box model to test effects of 

density variation with pressure.

TABLE 2.—Results of computer runs with density a function of

pressure

Run

1

2

3 _____

4

r —T <°C\ j — 1 o ^ *->)

100

1001 1

p

P(T)

P (T,P)

P(T) 

P (T,P)

Nu

2.60

2.60

2.60 

2.61

REALISTIC MODELS

Application of results from box model investigations are somewhat limited because the boundary conditions are not realistic. In particular, the 

assumptions of insulated sides and constant temperature at the top and bottom of the convecting layer may 

not hold in the field problem. It is possible, however, 

to use the numerical code to evaluate the effects of

these assumptions on temperature distributions and 

heat transfer rates.

CONDUCTING SIDE WALLS

Consider the model shown in figure 13 which is 

similar to one proposed by Donaldson (1968,1970) for 

simulation of convecting geothermal systems. The 

permeable channel, labeled B, could represent a fractured zone of upflowing liquid, separated from a recharge area with downflowing water by impermeable 

blocks, labeled A. A horizontal channel connecting the 

recharge and discharge areas is not shown at the bottom of this model because our interest here is on circulatory convection in the vertical channel. For this 

reason we also treat the upper and lower surfaces of 

Channel B as impermeable. A constant temperature 

T0 is imposed along the top and sides, and a linear 

temperature variation from T0 to Tl is imposed along 

the bottom of the blocks A. The bottom of channel B 

is maintained at T15 and the vertical sides of the channel are impermeable but thermally conductive.

For a particular case in which T^—T^ = 100° C, w/L 

= 2.1, and y/L = 0.5, and a Rayleigh number for the 

channel of 100, the temperature distribution shown in 

figure 14 would result. The temperature distribution 

for the case of conduction only (Ra = 0) is also shown 

and was used as a reference in computing the Nusselt 

number to account for the influence of the distance w 

and the temperature T0 imposed at the sides of the 

model. Heat transfer rates along various boundaries 

are listed in table 3.

The Nusselt number for the channel is 1.3, compared with a value of 2.2 determined for the corresponding box model with insulated sides as in figure 7.

bT-.T,-Tn 

dX w

FIGURE 13.—Physical model for cellular convection with conducting side walls.
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS

0

D-15

FIGURE 14.—Temperature distributions in convection model with conducting side walls for w/L = 2.1 and

ylL = 0.5.

TABLE 3.—Heat transfer results for convection in vertical channel with conducting side walls

Location Q for Ra = 0 

(joules/d X 10')

Top of channel 

Bottom of channel _ 

Top of model _ _ _ _ 

Bottom of model

4.38 

6.41 

24.4 

29.0

Q for Ra = 100 Nu 

(joules/d X 10')

5.76 

7.98 

26.3 

31.0

1.32 

1.24 

-1.08 

1.07

A simulation was also performed on the box model in 

figure 7 with a linear temperature variation imposed 

along the vertical sides. This simulation also allows 

heat conduction through the sides and for Ra = 100, 

Nu = 1.8.

These results indicate that analysis of circulatory 

convection in bounded reservoirs should provide for 

the retarding effect of heat conduction through the 

lateral boundaries. The basic physical model in figure 

13 can also be used to simulate the combination of 

convective throughflow from recharge to discharge 

areas and circulatory or secondary convection in the 

upflow channel (Donaldson, 1968,1970). Results produced by using the program are in agreement with 

the general features of Donaldson's combined 

throughflow and secondary convection results; these 

results indicate that the amount of heat transferred 

by the secondary flow decreases markedly as the 

throughflow increases. Numerical analysis of this basic model, including recharge channel (s) at the sides 

and a reservoir at the base, could utilize more complex 

boundary conditions and parameter variations to simulate conditions in a specific hydrothermal system.

TWO- AND THREE-LAYER MODELS

For the case of convection in a laterally extensive 

reservoir, we could remove the restrictive assumption 

of constant T along the upper and lower surfaces by 

considering a multi-layer system with less permeable 

rock above and below the reservoir. Convective flow in

the reservoir would not draw heat evenly from below 

and heat would not be conducted evenly across the top 

of the layer. Consequently, isotherms in the less permeable layers would be bent near the contacts with 

the permeable layer in contrast to the constant temperature contacts in the single-layer model.

Donaldson (1962) analyzed the magnitude of this 

effect for a two-layer system with a permeable layer 

underlain by an impermeable layer. Both zones had 

the same finite thickness and were of infinite lateral 

extent; the upper surface of the permeable layer and 

the lower surface of the impermeable layer were held 

at constant temperatures. Donaldson found that the 

value of Rac was lower in the double-layer case than in 

the single-layer case, but that the conducting layer 

tended to retard convection effects for values of Ra 

above the Rac .

Several runs using the program were made to compare the results with Donaldson's results and to 

extend the models to the three-layer case. The multilayer-models are illustrated in figure 15. An equivalent Rayleigh number was computed using the 

vertical temperature difference across the permeable 

layer which would exist with heat flow by conduction 

only and the thickness L of the permeable layer. In 

table 4, the steady state Nusselt numbers for each 

model for Ra = 100 are listed. The retardation effect 

noted by Donaldson is evident in the decrease in Nu 

from 2.60 for one layer to 1.19 for three layers; although for an aspect ratio of 1.0, the Nusselt number 

for the three-layer system may have been closer to 

that for the two-layer case.

TABLE 4.—Effect of multiple layers on heat transfer in permeable 

layer with circulatory convection

[S = cell aspect ratio]

Model

3 Layers --------

Ra

100

100 

100

s

1.0 

1.0 

.8

Nu

2.60 

1.34 

1.19
D-16 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

110 C

-185.0°-

Permeable

210°C

60°C

90-

-180-

Permeable

210 C

B

FIGURE 15.—A, Isotherms for two-layer model with Ra = 

100; B, Isotherms for three-layer model with Ta = 100.

The corresponding temperature distributions in the 

multilayered systems as plotted in figure 15 show the 

distortion of the isotherms in the impermeable layers 

near the contacts with the convecting layer. Near the 

constant temperature boundaries, the isotherms flatten out but are more closely spaced than if there were 

no convection. This effect raises the question of how 

deep to place the constant temperature lower boundary to minimize its influence on simulating the convective system. We can use an analysis by Birch (1967) 

of the temperature distribution beneath the sea floor 

to show that the isotherm distortion decays as 

e—irz/x', where z = distance from the convective 

boundary and x' is the distance from peak to trough 

for a sinusoidal isotherm at z = 0. For irz/x' = 2, the 

isotherm distortion would be less than 10 percent of 

the value at z = 0. Hence, distortion will be small at

the outer boundaries of the multilayer models if x'<z, 

that is if the thickness of the impermeable layer exceeds the peak to trough distance (the width of the 

convection cell). Since x' must be equal to or less than 

the thickness of the permeable layer L, z>L ensures 

sufficient thickness for the impermeable layer.

HOT SPRING MODELS

In many of the potential geothermal areas in California and Nevada, the convective and conductive 

heat flow in and around hot springs are the dominant 

features in the near-surface thermal regime and account for most of the heat discharge from these areas 

(Sorey and Lewis, 1975; Olmsted and others, 1975). 

Commonly, information from deep drilling as to the 

location and properties of the reservoir or reservoirs 

supplying the hot spring discharge does not exist. 

However, we would expect to find fluid and rock with 

relatively high heat content in and around the hot 

springs at.depths considerably shallower than the underlying reservoirs which may be at depths of 1-5 km 

(Olmsted and others, 1975). The shallow occurrence of 

hot fluid and rock suggests the possibility of extraction of energy from relatively shallow depths in these 

areas of natural discharge, either by capturing the hot 

spring discharge or the injection and withdrawal of 

fluid from another source. The feasibility of this approach has been demonstrated at the Casa Diablo Hot 

Springs in Long Valley, California (McNitt, 1963).

To describe the transfer processes associated with 

an upflowing hot spring system, the effects of conductive heat losses, mixing, subsurface discharge, and 

boiling need to be accounted for. Geochemical techniques, including those based on concentrations of 

silica and the ratio of concentrations of the cations 

Na, K, and Ca in the hot spring water can provide estimates of the source temperature at depth for the 

spring water (Fournier and others, 1974). In some 

cases, the amount of colder, fresher ground water mixing with the thermal water can be determined from 

the chemistry and enthalpy of the spring discharge. 

However, to quantify the amount of heat loss by conduction away from the spring conduit and the resultant temperature drop in the spring water, we need to 

use a numerical model. In addition, it is of interest to 

analyze the thermal and hydrologic regimes in the 

rock adjacent to the spring conduit under conditions 

of circulatory convection and in the absence of convective motions.

For the simplest model, we assume an isolated cylindrical conduit with vertical orientation surrounded 

by impermeable rock as in figure 16. The lower boundary is formed by the top of a reservoir with fluid at
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-17

FIGURE 16.—Isolated cylindrical conduit hot spring model.

temperature Tb and the upper boundary is the land 

surface at temperature Ts. The model is extended far 

enough laterally to justify an insulated boundary condition on the sides. Water enters the bottom of the 

conduit with temperature Tb at some constant mass 

flow rate pv^. As the fluid moves up the conduit, heat 

is conducted into the adjacent rock and the fluid temperature falls. The fluid emerges at the land surface at 

temperature Tsp and the conductive heat loss from 

the conduit is discharged at the land surface within an 

area of above normal heat flow. Note that the constant 

temperature Ts is imposed on the land surface up to 

the edge of the spring conduit but that the temperature of the spring discharge is not fixed. This condition is possible numerically by neglecting heat 

conduction across the spring surface.

Another possible model for heat transfer in a hot 

spring system involves water moving up over a more 

extensive fault plane and then converging near the 

surface as shown in figure 17. In a specific area, the 

applicability of the fault plane model over the cylindrical model can be judged from data on the thermal 

regime adjacent to the springs (Olmsted and others, 

1975) and from hydraulic considerations of upward 

velocities permitted by the available head and conduit 

permeability (Sorey and Lewis, 1975). It is important 

to make this distinction, because for the same total 

spring discharge, the plane model would have considerably greater area for conductive heat loss and lower 

fluid velocity than the cylindrical model, and therefore greater heat loss and temperature drop.

Mathematical analysis of these models is facilitated 

by combining the relevant parameters in energy balances for the upflowing fluid columns. Geometries 

and boundary conditions for each model are illustrated in figure 18. As an approximation, the fault

Hot Creek

FIGURE 17.—Fault plane hot spring model.

plane model is treated as 2-dimensional and plane 

symmetric, and the mass flux rate is computed as the 

total spring discharge divided by the area of upflow. 

For the cylindrical model, neglecting heat conduction 

in the vertical direction,

2irrcKri dT 

dr = pvdcirr* dT

dz (37)

which can be rearranged as

dz

dT

dr (38)

Substituting 6 = (T-Ts)/(Tb -Ts), r' = r/rc , and 

z' = z/L we obtain

Me dQ dQ

dz' df (39)

where Mc = pv^cr^ /2KmL = dimensionless mass flow 

rate in the cylindrical spring conduit. Integration of 

(39) from z' = 0 to z' = 1 yields

Mc (6ssp-l) = -- dz' = Qc (40)

where Qsp = (Tsp - Tb)/(Tb - Ts) = dimensionless 

spring temperature at land surface and Qc = dimensionless radial heat loss from the cylinder. A similar 

development for the fault plane model in figure 18 

yields

Mp(Qsp-l) = Qp (41)

where Mp = pvdcW2/KmL = dimensionless mass flow 

rate in fault plane conduit and Qp = dimensionless radial heat loss from plane.
D-18 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

Cylindrical model 1

/

Fau\t p\ane mode\ \

^

FIGURE 18.—Geometry of hot spring models.

For numerical analysis of heat transfer in these 

models, a grid of 150 nodes with grid spacing increasing with distance from the spring and depth below 

land surface was used. The fluid mass flux and temperature were prescribed at the base of the conduit 

and the same mass flux was removed at the top of the 

conduit. The initial temperature condition was the 

conduction-only solution with dT/dz = (Tb —TS )/L everywhere, and following the initiation of springflow, 

transient and steady state solutions were obtained for 

a range of mass flow rates. Results for cases where the 

rock surrounding the spring is impermeable are shown 

in figure 19, in terms of the relationships between Mc, 

Mp and (1 - 9sp).

For both models, as the flow rate goes to zero, 

(l-Qsp)=(Tb~Tsp)/(Tb -Ts) approaches 1.0, that is, 

the spring temperature approaches the land surface 

temperature. Unfortunately, the relationships expressed by equations 40 and 41 are not unique because 

Qc depends on L and rc and Qp depends on L and W. 

The solid curves shown in figure 19 were obtained for 

the parameters rc =6 m, W=10 m, and L=l km. Portions of additional curves for other parameter values 

as shown by the dashed lines indicate that variations 

in (l-9Sp) for a given dimensionless mass flow rate are 

within about 20 percent of the values computed from 

the solid curves for parameter values varying by a factor of 2.

For an example of the use of figure 19, consider 

spring flow in a cylindrical conduit of radius rc = 6 meters and depth L = 1 km and a fault plane of halfwidth W— 10 meters, length D = 1 km, and depth L = 

1km. For a total mass flow of 1.1 X 105 kg/day, with 

Tb -Ts = 180°C-10°C, Km = 2X 10'3 cal/(sec °C cm), 

and c = 1 cal/gm°C, the corresponding values of Mc 

and Mp are 1.00 and 0.031. From figure 19, (l-esp) = 

0.11 and 0.61 for the cylinder and plane, respectively. 

These values correspond to a temperature drop of

18° C and a spring surface temperature of 162° C for 

the cylindrical model and a temperature drop of 

104° C and spring surface temperature of 76° C for the 

plane model. The total lateral conductive heat loss 

from spring conduit can be computed for the cylinder 

as

Q = 2*KmL(Tb-Ts) Qc

and for the fault plane as

Q = 2KmL(Tb-Ts) (D/W)Qp

(42)

(43)

These equations yield, for the example above, heat 

losses of 2.3 X 104 cal/sec and 1.3 X 105 cal/sec for the 

cylinder and plane.

The steady state temperature distribution in the 

impermeable rock adjacent to the fault plane model 

considered above is shown in figure 20. The isotherm 

pattern for the cylindrical model is similar in cross 

section. Distortion of the isotherms due to the spring 

flow occurs in a zone within about 1 km of the fault, 

and the excess heat is conducted upwards toward the 

land surface. As anticipated, conductive heat losses 

from the conduit are considerably larger in the upper 

portion near the land surface than at depth. For this 

reason, the boundary conditions at the land surface 

should have a significant influence on spring temperatures near the surface. One complication to consider is 

that by imposing a constant temperature Ts at the 

land surface, a singularity in heat flow exists at the 

wall of the conduit. Numerically, the infinite temperature gradient is avoided by averaging over finite distances; but in actual hot spring areas, the land surface 

temperature increases considerably near the spring. 

Thus, a more realistic boundary condition would be 

the linear heat transfer or radiation boundary 

condition

Km ^-=HT,atz = L (44)
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-19

FIGURE 19.—Relationships between dimensionless flow rates (Mc , Mp) and dimensionless temperature drop (1—Qsp) due 

to conductive heat loss in cylindrical and fault plane hot spring models.

0

10°C 1 km 2 km

1 km 180 C

FIGURE 20.—Steady state temperature distribution in fault plane hot spring model with discharge = 105 kg/d, W = 10 m, D

1 Km, and Ra = 0 (k = 0).
D-20 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

where H = coefficient of surface heat transfer (Carslaw and Jaeger, 1959, p. 19). This radiation boundary 

condition reduces to the constant temperature boundary condition as H—»oo and to an insulated boundary 

condition as H —»0. An upper limit for H could be obtained from experimental data on turbulent air flow 

over flat plates. For example, Rohsenow and Choi 

(1961, p. 200) give

H = 1.5 x 10-8 (pci>)< (45)

which for air at 15° C yields the following relationship 

at the land surface between H and mean annual wind 

velocity v.

(cm/sec)

H

(cal/sec °C cm2 )

10

100

1000

4.5 X 10'6 

4.5 X 10"5 

4.5 X 10"4

Other factors complicate the choice of an effective 

value of H, such as intermittent air flow, evaporation, 

and the presence of a dry, low conductivity soil layer. 

In the latter case, we could relate a particular value of 

H to a low conductivity layer of thickness 1 = Km/H. 

For example, for Km = 0.5 X 10'3 cal/sec °C cm, 1 

ranges from 0.01 to 1.1 meters for H between 4.5 X 10-4 

and 4.5 X 10"6 cal/sec °C cm2 , respectively. In this 

study, the effect of using the radiation boundary condition (44) was examined for a range of H from IO-6 to 

10-4 cal/(sec °C cm2 ).

Results are tabulated in table 5 for the fault plane 

case considered previously with spring discharge 

equal to 1.1 X 105 kg/day.

TABLE 5.—Spring temperature, Tsp, total conductive heat 

loss, Q, and heat flux at land surface near spring, q, for selected values of coefficient of surface heat transfer H

H

(cal/sec °C cm2 )

io-4

io-5

in-6

TS n (°C)

7C Q

7fi Q

7ft Q

n-i Q

Q

(cal/sec)

1 on v 1O5

1 on v 1O5

i 97 y in5

1 -iq v If)5

q(x = 5m) 

(li cal/sec-cm2 )

QO

7970

43

For values of H greater than about 10~4 , the radiation boundary condition yields essentially the same 

results as the constant temperature boundary condition. However, for H < 10"5 , the spring temperature is 

significantly greater and the total conductive heat loss 

and land surface heat flux near the spring less than for 

the constant temperature case. The temperature distribution for H = 10"6 cal/sec °C cm2 is shown in figure 

21. Land surface temperatures and heat flows vary 

from 53° C and 43 HFU at a distance of 5 meters to 

13.5° C and a normal conductive heat flux of 3.5 HFU 

beyond about 1 km from the spring. Thus, if the effective values of heat transfer coefficient in field areas 

are below 10" 5 , the radiation boundary condition 

should be used in simulating the near surface thermal

1 km 2 km

1 km 180 C

FIGURE 21.—Steady state temperature distribution in fault plane hot-spring model with discharge = 105 kg/d, W = 10 m, D • 1 Km, and Ra = 0 (k = 0), and H = 10"6 cal/sec0 C cm2 at land surface.
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-21

regime. It should also be noted that we are ignoring 

the details of heat and mass transfer at the surface in 

the spring. Many factors such as the geometry of the 

spring orifice, presence or absence of lateral discharge, 

and evaporative cooling may cause considerable variations in actual temperatures near the surface.

It is of interest to consider the time required for the 

conductive thermal regime to reach equilibrium following the development of the spring system. From 

equation 32, the effective time constant for this 2-dimensional heat transfer problem could be approximated by the expression L2 ~pc/4Km, where L2 can be 

thought of as the area through which heat flows from 

the spring to the land surface. For example, from figure 20, L2 might represent an area 1 km deep by 1 km 

wide or 1 km2 . This means that the time constant in 

this case would be about 9 X 106 days or 25,000 years. 

The transient numerical simulations for this problem 

show that temperatures are within 0.1° C of equilibrium after a period of about 50,000 years, or twice the 

effective time constant. The thermal diffusivity in 

this problem was a relatively low 0.0032 cm2/sec; equilibrium times for larger diffusivities would be proportionately smaller.

If the restriction is removed that the rock adjacent 

to the spring is impermeable, ground-water circulation in this region will affect the thermal regime. Even 

if the rock permeability is low enough that the critical 

Rayleigh number of 4ir2 is not exceeded, the fluid

would still be unstable because of the high temperature boundary formed by the spring conduit. In modeling this problem numerically, the simplification is 

made that there is no hydraulic connection between 

the spring conduit and the adjacent ground-water system. To simplify the problem, a thin impermeable 

zone is included at the contact between the two regions; such a zone may exist in many field situations 

where chemical deposition has sealed off the conduit 

from the surrounding rock.

Once again, consider the fault plane model analyzed 

previously, that is, discharge = 1.1 X 105 kg/day, 

Tb -Ts = 180°C - 10°C (fig. 20), but with a rock permeability of 1 millidarcy (10' 15 m2 ). In the absence of 

spring flow, the equivalent Rayleigh number for this 

case would be 20. With spring flow, the steady state 

temperature distribution shown in figure 22 would result. The existence of convection cells is not clearly 

delineated by the isotherm pattern although comparison with the 125°C and 150°C isotherms from figure 

20 show the shifts caused by upflowing fluid near the 

conduit and downflowing fluid to a distance of about 

1.5 km. The actual velocity distribution includes 3 

convection cells but with the magnitude of convection 

decreasing with distance from the spring. As a result 

of the convection, the temperature drop and total conductive heat loss from the spring are about 6 percent 

less than for the corresponding impermeable rock 

case. The 25° C isotherms in figure 22 indicate that

10°C 1 km 2 km

-175-

1 km 180 C

FIGURE 22.—Steady state temperature distribution in fault plane hot-spring model with discharge = 105 kg/d, W = 10 m, 

D = 1 Km, Ra = 20 (solid lines), Ra = 0 (dashed in fig. 20).
D-22 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

heat flow near the land surface is not appreciably affected by convection in the rock below. Similar results 

in terms of slight lowering of the conductive heat loss 

from the conduit and negligible changes in heat flow 

near the land surface were found for other spring discharge rates in systems with Rayleigh numbers below 

Rac although the distortion of the temperature distribution at depth increases somewhat for the higher 

spring discharge rates.

Now, if the rock permeability is increased to 10 millidarcies (10~14 mz ), the Rayleigh number would be 200; 

and for a discharge of 1.1 X 105 kg/d the resultant temperature distribution is plotted in figure 23. In this 

case the influence of cellular convection set up by the 

constant temperature boundaries at the top and bottom of the model dominates the thermal regime. The 

effect of the hot spring is evident out to a distance of 

about 0.5 km over which the 100°C and 125°C isotherms show reversals characteristic of cellular convection at large Ra. In this case, convection in the rock 

causes a more substantial lowering of the conductive 

heat loss from the spring conduit—0.99 X 105 cal/sec 

versus 1.29 X 105 for the Ra = 0 case—and a temperature drop of 84° C compared with 104° C for Ra = 0. 

The relative effects of the conductive heat loss from 

the conduit and the enhanced heat transfer caused by 

convection in the rock are indicated near the land surface by the more closely spaced isotherms between 0 

and 0.5 km from the spring compared with the spacing 

at corresponding intervals beyond 1 km.

These results represent preliminary attempts to 

simulate heat transfer and fluid flow associated with 

hot spring systems. To model specific field areas, it 

may be necessary to include nonvertical orientation of 

the spring conduit, subsurface discharge and recharge, and boiling near the land surface in the simulation. This work is being planned as sufficient data 

become available from field studies in California and 

Nevada. As noted previously, energy development in 

many of these geothermal areas may focus on localized 

regions of natural discharge. Thus a useful extension 

of the numerical modeling work described here will be 

the evaluation of various schemes for energy development.

SUMMARY AND CONCLUSIONS

We have discussed the development of a numerical 

code which can treat problems involving slightly compressible fluid and heat transfer in multidimensional 

porous media. Solutions to the appropriate partial differential equations are obtained by the integrated finite difference method which is essentially equivalent 

to balancing mass and energy over finite subregions or 

elements. The resultant system of finite difference 

equations is solved by an iterative procedure, and solutions to the fluid flow and energy equations are 

coupled by interlacing in time so that the temperature 

and velocity fields are interdependent. The useful 

concepts of fluid and thermal time constants as indi2 km

1 km 180 C

FIGURE 23.—Steady state temperature distribution in fault plane hot-spring model with discharge = 105 kg/d, W = 10 m,

D = 1 Km, and Ra = 200,
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-23

cators of nodal response times and numerical stability 

limits are an inherent part of the numerical scheme.

In applying the numerical model to the problems of 

circulatory convection in saturated porous media, the 

relevant aspects have been discussed as they pertain 

to geothermal systems, and we find that results from 

the program on the relationship between the Rayleigh 

number and the dimensionless heat transfer coefficient or Nusselt number are in good agreement with 

numerical and experimental results from other authors. Then the numerical model was used to extend 

these results to include the effects of temperature dependent parameters and density variations with pressure. Variations in fluid viscosity and thermal 

expansivity with temperature result in substantial 

differences in the values of the critical Rayleigh number for the onset of convection and the Rayleigh number-Nusselt number relationship compared with corresponding constant parameters results. However, 

consideration of fluid density as a function of pressure 

produced no noticeable effect on convective motion.

Numerical simulations of more realistic models for 

circulatory convection show that for laterally bounded 

reservoirs, conduction of heat across the vertical side 

walls results in significant lowering of the rate of vertical heat transfer through the reservoir. For a laterally 

extensive reservoir, consideration of impermeable or 

less permeable layers above and below tide convecting 

layer removes the restrictive assumption of constant 

temperature boundaries on the permeable layer and 

has the effect of lowering the value of the critical Rayleigh number Rac while retarding convective heat 

transfer at values of Ra above Rac.

Heat and mass transfer associated with hot spring 

systems was analyzed to determine the amount of 

heat lost by conduction to the rocks surrounding the 

spring conduit. An isolated cylindrical conduit model 

and a fault plane conduit model were considered, and 

the temperature drop in the hot spring water between 

the source reservoir and the surface due to conductive 

heat loss was determined numerically as a function of 

flow rate. The steady state temperature distribution 

for the case where the rock surrounding the spring is 

impermeable shows that heat loss from the spring distorts the normally horizontal position of the isotherms 

out to distances comparable to the depth of the spring 

conduit. Conductive heat flux at the land surface is 

high near the spring but near the normal or background level beyond one conduit depth. The time required for the conductive thermal regime to 

equilibrate following the development of hot spring 

activity can be approximated by the expression L2pc/ 

2Km where L is the depth to the source reservoir. For

unconsolidated sediments with low thermal conductivity, the equilibration time is about 50,000 years for 

a reservoir at 1 km.

The effects of fluid circulation in the rock surrounding the spring conduit were examined for systems in 

which the equivalent Rayleigh number (in the absence 

of hot spring activity) was both above and below the 

critical value of 47r2 . With Ra<4ir2 , circulatory convection is set up due to the presence of the hot spring, but 

it causes only slight effects on the thermal regime in 

the rock surrounding the spring conduit and on the 

conductive heat loss and temperature drop associated 

with the spring. For the case with Ra>4Trz, circulatory 

convection resulting from the temperature difference 

Tfc-Ts between the source reservoir and the land surface dominates the thermal and hydrologic regimes 

and significantly reduces the conductive heat loss and 

temperature drop for the spring.

The results of this investigation demonstrate the 

usefulness of numerical modeling to describe the natural conditions of heat transfer and fluid flow in 

geothermal areas. Given preliminary thermal, hydrologic, geophysical, and geochemical information, it is 

possible to construct and analyze simplified conceptual models of specific hydrothermal systems as a guide 

to further data collection in undeveloped areas. An 

application of this technique applied to the Long Valley caldera in California is described by Sorey, Lewis, 

and Olmsted (1978). As sufficient parametric and geometric data become available from deep test drilling 

to refine the model, it can then be used to evaluate 

rates and magnitudes of energy recovery and effects of 

development on existing water resources.

REFERENCES

Bear, J., 1972, Dynamics of fluids in porous media: New York,

American Elsevier Pub., 764 p. 

Beck, J. L., 1972, Convection in a box of porous material saturated

with fluid: Physics of Fluids, v. 15, no. 8, p. 1377-1383. 

Birch, F., 1967, Low values of oceanic heat flow: Jour. Geophys.

Research, v. 72, no. 8, p. 2261-2262. 

Bird, R. B., Stewart, W. E., and Lightfoot, E. N., 1960, Transport

Phenomenon: New York, Wiley, 780 p.

Bodvarsson, G., 1961, Physical characteristics of natural heat resources in Iceland: U.N. Conf. New Sources Energy, Rome,

1961, Proc., v. 2, p. 82-89. 

Caltagirone, J. P., Cloupean, M., and Combarnous, M. A., 1971,

Fluctuating natural convection in a horizontal porous layer [in

French]: Acad. Sci., [Paris], Comptes rendus, B, v. 273, no. 20,

p. 833-836. 

Carslaw, H. S., and Jaeger, J. C., 1959, Conduction of heat in solids:

London, Oxford Press, 510 p. 

Coats, K. H., George, W. D., Chu, C., and Marcum, B. E., 1973,

Three-dimensional simulation of steam flooding: Soc. Petrol.

Engr. Jour., v. 15, no. 6, p. 573-592.
D-24 GEOHYDROLOGY OF GEOTHERMAL SYSTEMS

Combarnous, M. A., and Bories, S. A., 1973, Thermal convection in 

saturated porous media: Groupe D'etude Sur Les Milieux Poreux, Institut Francais Du Petrole, Toulouse, Rapport G. E. 11.

Curry, D. M., 1970, Multi-dimensional analysis of heat and mass 

transfer in porous media: Houston, Texas, University of Houston, unpublished Ph. D. thesis, 197 p.

Darcy, H. P. G., 1856, Les Fontaines Publignes de la Villa de Dijon: 

Paris, Victor Dalmont, 647 p.

Donaldson, I. G., 1962, Temperature gradients in the upper layers 

of the earth's crust due to convective water flows: Jour. 

Geophys. Res., v. 67, no. 9, p. 3449-3459.

——— 1968, A possible model for hydrothermal systems and 

methods of studying such a model: Third Australian Conf. on 

Hydraulics and Fluid Mechanics, Nov. 25-29. Paper 2580. 

1970, The simulation of geothermal systems with a simple

convective model: Geothermics, Spec. Issue 2, U.N. Symposium on the Development and Utilization of Geothermal Resources, Pisa, Italy, v. 2, pt. 1, p. 649-654.

Dorsey, N. E., 1968, Properties of ordinary water-substance: New 

York, Hafner Pub. Co., 673 p.

Edwards, A. L., 1969, TRUMP: A computer program for transient 

and steady-state temperature distributions in multidimensional systems: Natl. Tech. Inf. Service, Springfield, Virginia, 

UCRL 14754, Rev. 2, 22151, 263 p.

Einarsson, T., 1942, The nature of springs of Iceland (in German): 

Rit. Visind. Es., 26, 91 p.

Elder, J. W., 1957, Some problems in fluid dynamics: Cambridge, 

England, Cambridge University, unpublished Ph. D. thesis.

___ 1966, Heat and mass transfer in the earth: Hydrothermal 

Systems, New Zealand, D.S.I.R. Bull. 169,115 p.

Faust, C. R. and Mercer, J. W., 1975, Mathematical modeling of 

geothermal systems [abs.]: Proc. of Second U.N. Symposium 

on Development and Use of Geothermal Resources, San Francisco, Calif., May 20-29, 1975, v. 3, p. 1635-1641.

Fernandez, R. T., 1972, Natural convection from cylinders buried 

in porous media: Berkeley, California, University of California, unpublished Ph. D. thesis, 300 p.

Fournier, R. O., White, D. E., and Truesdell, A. H., 1974, Geochemical Indicators of subsurface temperature—Pt. 1, Basic 

assumptions: Jour. Research U.S. Geol. Survey, v. 2, no. 3, p. 

259-262.

Green, D. W., 1963, Heat transfer with fluid flowing through porous media: Oklahoma, University of Oklahoma, unpublished 

Ph. D. thesis.

Henry, H. R. and Hilleke, J. B., 1972, Exploration of multiphase 

fluid flow in a saline aquifer system affected by geothermal 

heating: University of Alabama Bureau of Engr. Research 

Rept. No. 150-118,123 p.

Hoist, P. H., 1970, A theoretical and experimental investigation of 

natural convection in porous media: University Calgary, Alberta, Canada, Ph. D. thesis.

Hoist, P. H., and Aziz, K., 1972a, A theoretical and experimental 

study of natural convection in a confined porous medium: 

Canadian Jour. Chem. Eng., v. 50, p 232-241.

___ 1972b, Transient three-dimensional natural convection in 

confined porous media: Int. Jour. Heat and Mass Transfer, v. 

15, p. 73-90.

Home, R. N. and O'Sullivan, M. J., 1974, Oscillatory convection in 

a porous medium heated from below: Jour. Fluid Mech., v. 66, 

pt. 2, p. 339-352.

James, R., 1968, Wairakei and Larderello; geothermal power systems compared: New Zealand Jour. Sci. and Technology, v. 11, 

p. 706-719.

Landau, L. D., and Lifshitz, E. M., 1959, Fluid Mechanics: London, 

Pergamon Press, 536 p.

Lantz, R. B., 1971, Quantitative evaluation of numerical diffusion 

(truncation error): Soc. Petroleum Engineers Trans., v. 251, p. 

315-320.

Lapwood, E. R., 1948, Convection of a fluid in a porous media: 

Proc. Cambridge Phil. Soc., v. 44, p. 508-521.

Lasseter, T. J., 1974, Underground storage of liquified natural gas 

in cavities created by nuclear explosives: Berkeley, California, 

University of California, unpublished Ph. D. thesis, 273 p.

Lasseter, T. J., Witherspoon, P. A., and Lippman, M. J., 1975, Multiphase multidimensional simulation of geothermal reservoirs: 

Proc. of Second U.N. Symposium on Development and Use of 

Geothermal Resources, San Francisco, California, May 20-29, 

1975, v. 3, p. 1715-1723.

McNitt, J. R., 1963, Exploration and development of geothermal 

power in California: California Div. Mines and Geology, Spec. 

Rept. 75.

Mercer, J. W., 1973, Finite element approach to the modeling of 

hydrothermal systems: Urbana-Champaign, Illinois, University of Illinois, unpublished Ph. D. thesis, 115 p.

Mercer, J. W., Faust, C., and Pinder, G. F., 1974, Geothermal reservoir simulation: Proc. of National Science Foundation Conference on Research for the development of Geothermal Energy 

Resources, Pasadena, California, Sept. 23-25, 1974, p. 256- 

267.

Narasimhan, T. N., 1975, A unified numerical model for saturatedunsaturated groundwater flow: Berkeley, California, University of California, unpublished Ph. D. thesis, 244 p.

Olmsted, F. H., Glancy, P. A., Harrill, J. R., Rush, F. E., and Van 

Denburgh, A. S., 1975, Preliminary hydrologic appraisal of selected hydrothermal systems in northern and central Nevada: 

U.S. Geol. Survey Open-File Rept. 75-56, 267 p.

Platzman, G. W., 1965, The spectral dynamics of laminer convection: Jour. Fluid Mech., v. 23, p. 481-510.

Raats, P. A. C., and Klute, A., 1968, Development of equations describing transport of mass momentum and energy in soils, in 

Water in the unsaturated zone, Proc. Wageningen Symposium: Pub. No. 83, International Association of Scientific Hydrology, p. 535-546.

Remson, L, Hornberger, G. M., and Molz, F. J., 1971, Numerical 

methods in subsurface hydrology: New York, Wiley-Interscience, 389 p.

Rohsenow, W. N., and Choi, H. Y., 1961, Heat, mass and momentum transfer: Englewood Cliffs, New Jersey, Prentice-Hall.

Sorey, M. L., and Lewis, R. E., 1975, Convective heat flow from hot 

springs in the Long Valley caldera, Mono County, California: 

Jour. Geophys. Research, v. 81, no. 5, p. 785-791.

Sorey, M. L., Lewis, R. E., and Olmsted, F. H., 1978, The Hydrothermal System of Long Valley caldera, California: U.S. Geol. 

Survey Prof. Paper 1044-A, 60 p.

Spillette, A. G., and Nielson, R. L., 1968, Two-dimensional method 

for predicting hot waterflood recovery behavior: Jour. Petroleum Technology, v. 20, p. 627-638.

Sutton, F. M., 1969, Stability of flow in a model of a hydrothermal 

system: Victoria University, Wellington, New Zealand, unpublished M. Sc. thesis.

Weinstein, H. G., Wheeler, J. A., and Woods, E. G., 1974, Numeri-
NUMERICAL MODELING OF LIQUID GEOTHERMAL SYSTEMS D-25

cal model for steam stimulation in Improved oil recovery symposium, Apr. 22-24, 1974: Tulsa, Oklahoma, Soc. Petroleum 

Eng. preprint 4759, p.,253-270.

White, D. E., 1957, Thermal waters of volcanic origin: Geol. Soc. 

America Bull, v. 68, p. 1637-1658.

___ 1961, Preliminary evaluation of geothermal areas by geochemistry, geology, and shallow drilling, U.N. Conf. New 

Sources Energy, Rome, Italy, 1961, Proc. v. 2, p. 402-408.

White, D. E., Muffler, L. J. P., and Truesdell, A. H., 1971, Vapordominated hydrothermal systems compared with hot-water

systems: Econ. Geology, v. 66, no. 1, p. 75-97.

Witherspoon, P. A., Newman, S. P., Sorey, M. L., and Lippman, M. 

J., 1975, Modeling geothermal systems: Univ. California, Lawrence Berkeley Lab. Rept. 3263, 68 p.

Wooding, R. A., 1957, Steady state free thermal convection of liquid in a saturated permeable medium: Jour. Fluid Mech., v. 2, 

p. 274.

___ 1963, Convection in a saturated porous medium at large 

Rayleigh number or Peclet number: Jour.Fluid Mech., v. 15, p. 

527.