- Email: [email protected]
Heat: A fortran computer program for calculating 1-D conductive and advective heat transport in geological formations
Computers & Geosciences Vol. 19, No. 5, pp. 673~689, 1993
0098-3004/93 $6.00 + 0.00 Copyright c 1993 Pergamon Press Ltd
Printed in Great Britain. All rights reserved
HEAT: A FORTRAN COMPUTER PROGRAM FOR CALCULATING 1-D CONDUCTIVE AND ADVECTIVE HEAT TRANSPORT IN GEOLOGICAL FORMATIONS PATRICE DE CARITAT Department of Geology and Geophysics, University of Calgary, Alberta, Canada T2N IN4
(Received 6 April 1992; revised 9 July 1992; accepted 24 November 1992)
Abstract--A simple, accurate, numerical approximation of the one-dimensional equation of heat transport by conduction and advection is presented. The method is based on the iterative solution of an implicit, finite difference, Crank-Nicolson algorithm, featuring alternating differencing direction as a function of the thermal Peclet number, and is implemented in a documented FORTRAN code. The program allows users to evaluate how porewater flow may affect the thermal profile of permeable rock formations, under steady-state or unsteady-state boundary conditions, and with linear or nonlinear initial conditions. Several possible applications illustrate the potential usefulness of the program in first order, evaluative investigations. Key Wordr: Heat transport, Advection, Conduction, Finite difference, Numerical modeling.
INTRODUCTION Heat transport in saturated porous media, such as rock formations below the watertable, takes place as a result of conduction and, if the contained porewater is moving sufficiently fast, advection. The applications of heat transport to geological processes are numerous and important, for instance in metamorphic, ore deposit, hydrogeoiogical, thermal, and diagenetic studies. Geoscientists who are not inclined numerically, however, may be discouraged to tackle this issue because of the apparent complexity of the mathematical and numerical manipulations involved. The purpose of this communication is to provide these geoscientists with a simple, one-dimensional (l-D) computer simulation program, which will perform an accurate evaluation of conductive and advectire heat transport in porous media. HEAT TRANSPORT
IN P O R O U S
MEDIA
In a saturated permeable medium, such as a porous or fractured sedimentary rock, heat is transported by conduction ("diffusion") through the solid framework and within the fluid phase. If the contained fluid is moving through this permeable medium, advection ("convection") also contributes to heat transport. The relative importance of conductive versus advective heat transport is a function of the thermal diffusivity of the rock and fluid flow rate (Bickle and McKenzie, 1987). The principles of heat transport are presented only briefly herein, and a simple and accurate numerical solution is presented. Transport in the vertical direc-
tion will be specifically considered, although this solution applies to heat transport in any direction. As it would not be possible nor appropriate to present, within the context of this communication, a thorough review of the vast literature pertaining to heat transport, the reader is referred to Carslaw and Jaeger (1959) for an introduction to the physics of heat transport, and to Bickle and McKenzie (1987) for a presentation of geological applications, among many others. The conductive/advective heat transport equation, in its simplest form and in one-dimension, is given by: c?T O,t
k ?,2T_ pfQ -
= ~-
OT vf0 ~-,
p~c~ cz"
p~c~
cz
(I)
where T is the temperature (K), t the time (sec), k the thermal conductivity of the solid medium (W m 1K-J), Ps and pf are the solid and fluid specific density (kg m -3), respectively, cs and cf the solid and fluid specific heat capacity (J k g - ~ K - ~), respectively, vf is the fluid velocity (m sec-~), negative for upward flow, 4) the porosity of the solid medium, and z the direction in which heat transport is considered, here vertical (see Table 1 for a listing of symbols. suggested values and units). The following substitutions: ~," -
(2)
pfcf , P~Cs
(3)
W ~ ~t,~9,
(4)
M
673
k p~c~ ' ~ -
674
P. DE C^R,r^r Table 1. Symbols, suggested values and units (after Loosveld, 1989) symbol
meaning
value (range)
unit
ct" ca m
~ i I i c Imatc ~ c i t y of fluid ~aceific heat capacity of solid time-step
3800-4000 1000-1300 afleulamt
Liz
dcpth-st¢~
I~¢-~
J.kg-l.K-1 J.kg-l.K-1 see m
¢
r,om~ty
length of interest (for Pet) thermal diffusivity of solid
- 0.0-o.25 pre-selected
--
h ~'= k
- 10-7-10"6
m2.scc-1
2.0-3.0
W.m-1.K-1
Ps Cs
k M
thermal conductivity of solid
= Pfct
Ps cs Wh Pe, It"
r =
I¢ At
A z2
pf Ps
T t
v/ W = vf#
z
d2T
1, calculated
-
--
thermal Peclet ntmaber
calculated
--
stability factor
0.4 (fixed)
--
specific density of fluid specific denfity of solid temperature time fluid velocity effective fluid flux direction considemt
1000-1010 2700-3300 --prc-sel~"ted calcula~ --
kg.m-3 kg.m-3 K see
where x is termed the thermal diffusivity of the solid medium (m 2 see-t), and 14," the effective fluid flux (m sec- ~), result in the following simplified expression of Equation (I): dT
m
OT (5)
Ot = ~-ff~z2 - M W - ~ z "
This is the partial differential equation describing 1-D conductive/advective heat transport in porous media. It states that the rate of change of temperature at any point is equivalent to the sum of conductive heat transport (first term of the right-hand side) and advective heat transport (second term of the righthand side). In this form, it represents a simplification of the real-world situation, as it neglects internal heat generation, thermal expansion, compressibility, the temperature-dependency of heat capacity and thermal conductivity, etc. Furthermore, it is applied here to a homogeneous (nonlayered) 1-D medium. A time-dependent solution for this type of equation cannot be obtained analytically, except for the simplest of situations. However, a solution for the heat transport equation in the general transient (unsteadystate) situation can be approached satisfactorily using the Crank-Nicolson approximation (see Carnahan,
m.sec -I
m.sec -I m
Luther, and Wilkes, 1969; Croft and Liiley, 1977; Gerald and Wheatley, 1984). A computer program, HEAT, was designed to implement this procedure. In order to preserve the generality of this computer program, applicability to both slow and fast flowing systems was considered. In particular, provision was made to alternate the direction of the differencing scheme for the advection term as a function of the thermal Peclet number Pe, (Pe t = W h / x , where h is the length-scale of interest), as suggested by Spalding (1972). Thus, the final Crank-Nicolson approximation of Equation (5) is given by: i [TT÷ ~ _ T~] = x [T;':~ ~- 2T7 +~ + T;'_+~' at
~
"
+ T~+l - 2 T 7 + T'~ " Az 2 !I - MW
az 2
]
"+' ,,+l -o')T~+,-(20-I)T i +
4 (1 - tT)T'i' + 1-- (2o"2_Az-I)T7 + t;T'i'-I ],
Pet < -2 -2 < Pet < 2 2 < Pet
~valuc 0.0 0.5 1.0
(6)
..I
where a is the parameter that alternates the direction of the differencing scheme for the advection term as given in Table 2. In this formulation, the space and
Table 2. Conditions governing the direction of differencing scheme condition
~T,,+, i-l
y~
din~clion downmtam weighting midpoint weighting upstream weighting
Heat transport in geological formations time grids are represented by i (subscript) and n (superscript), respectively, and Az and At are the depth- (or space-) and time-steps, respectively. We now introduce r, the stability factor (or Fourier number): xAt r - Az"'
(7)
and #:
MWAt
U= - - , Az
675 t
.,,
at
IAz
l
IN
T~n+1 • ~7"~nlahi-1 qF
(8)
dhl qw
+ TT+l + r,
5+
-5
! -r
+T,'.' i
s'~l %J
'T.. n T.,n+1 ~i+1 ~i+1
and, rearranging Equation (6), we obtain the temperature at the ith depth-step and n + lth time-step:
r?+,=) ,+, 1_5-5' - 0 ) +
IF
En Tn+1
-5(2a
[-;
+
- 1) l+r+
]
(2~r-l) .
(9) Selecting r = 0.4 results in an unconditionally stable solution (Loosveld, 1989). The solution of Equation (5) is satisfactorily approximated by iteratively solving for the implicit, finite difference, Crank-Nicolson algorithm [Eq. (9)]. Figure 1 illustrates the grid or node system used here. Time increases from left to right, and depth increases downward. The t - z continuum under consideration is discretized using At and Az as time- and depth-steps, respectively. As can be seen in Figure i, the solution for the temperature at any node depends on the temperature at closest neighboring nodes at the same time-step, and on the temperature at the same and neighboring depth-step at the previous time-step. Such a configuration, which is referred to as implicit (Carnahan, Luther, and Wilkes, 1969), has the advantage of yielding stable results for any value of r, although small values for r ( ,~ 1) give more accurate results (Croft and Lilley, 1977; Gerald and Wheatley, 1984). Selected boundary conditions can be chosen to be the temperature at two fixed levels at the extremities of the space coordinate of interest (Dirichlet-type boundary conditions). A practical selection of boundary conditions would include the temperature at the surface (z = 0) and the temperature at the deepest level of interest (z = Zm,~), for instance the basement of a sedimentary basin, or at the base of the lithosphere. The selection of boundary conditions is a function of the scale of interest and the problem at hand. However, it is desirable to select the high temperature boundary condition to be as great as possible (greater depth, higher temperature), so as to avoid interference of this boundary within the (shallower) region of interest. In HEAT, the value of the
Figure 1. Discretization of time space (t :) continuum for finite difference approximation scheme. Approximation of temperature at any point (O) depends on temperature at five neighboring nodes (O) in t-: space, this being an implicit numerical scheme. two boundary conditions can be specified as a function of time, so that transient analysis can be performed. The initial condition may be taken to be the linear temperature gradient (conductive geothermal gradient) between these two levels. The program allows for a nonequilibrium (nonlinear) initial temperature profile to be specified. Figure 2 illustrates, in nondimensionalized terms, the effect of various conditions of upward, vertical porewater flow upon the equilibrium thermal profile (t = t~ ) in sedimentary basins, with respect to the linear conductive geotherm. APPLICATIONS To illustrate the purpose of this program, a brief description of a problem to which it was applied
0.0
I
,
J
0.2
,¢
0.4
06
o
08
,
?o,=+a5
• Pet=-17.6 1.0
'
0.0
' 0.2
'
' 0.4
"~'~ '
' 0.6
Temperature
'
' 0.8
'
-~1.0
Figure 2. Nondimensionalized representation of perturbations to conductive geotherm (11) as a result of upward, vertical porewater flow characterized by decreasing Pe, values (upward fluid flow requires negative sign).
676
P, DE CARITAT
follows. Baker and Caritat (1992) contended that ankerite cementation of the Aldebaran Sandstone, the main gas reservoir in the Denison Trough, eastern Australia, occurred as a result of vertical, upward, cross-formational porewater flow. This flow system developed as a thick sequence of overpressured continental mudrocks and coals underlying the Aldebaran Sandstone underwent overpressure relaxation, as a result of intensifying tectonic stress during the Late Triassic, immediately prior to basin inversion. The computer program H E A T was used to determine how much heat would have been advected under these circumstances as a function of time, and this constraint was used subsequently to model the temperature-sensitive fractionation of oxygen isotopes between the porewater and the precipitating ankerite (Caritat and Baker, 1992). Selected boundary conditions were the temperature at the surface and at the base of the sedimentary sequence, and initial condition was taken to be the conductive geotherm (25-3OC/km). Vertical flow velocities were estimated by considering the pressure differential and the hydraulic conductivity across the aquitard capping the overpressured section. Under the particular conditions of the Denison Trough, it was determined that ambient temperatures in the Aldebaran Sandstone were increased because of the upward advection of hotter water by 5-10' C relative to conductive temperatures of 120-140 C (depending on the location in the Denison Trough) during akerite cementation (Caritat and Baker, 1992). Other applications may include the conductive dissipation of a thermal anomaly caused by overthrusting (Fig. 3), the conductive/advective relaxation of a temperature spike introduced by sill emplacement or discharge of hot (hydrothermal) water in a permeable (p.e., fractured) horizion (Fig. 4), and the
0
200
400
Temperature 600 800 0
(°C) 200 400
600
2o
800 B
10 (conduction on~y) Imtantaneouscrusp,I ~ 20L ~zthl- c2koni n gbythrul C i n km ~~g~
o.,
o
0
100
Temperature (°C) 200 300 120 140 160 180 200
:up~va flow (.i0.o m.u¢,-1)~ 2 ~nltlal"hot"layer (200"C1,,t ~ 4.75 to 5.25 krn dq~h 3 ~z.2sOm \
~ 4
5 6
7
8 9 10
...........
Figure 4. A--Dissipation of thermal spike (200 C) in initial geotherm, mostly as result of conduction. This is followed by bowing-upward of geotherm through advection by upward fluid flow; B--details (geotherms are shown at 6 , 1, 2, 4, 6, 10, 50kyr, and t~).
propagation of a transient thermal anomaly in a tilted aquifer (Fig. 5). THE COMPUTER PROGRAM HEAT, a F O R T R A N computer code written tosolve Equation (9), is presented in Appendix 1. The code is annotated abundantly to allow the user to follow its logic, and is available from the author upon submission of a 3½in, floppy diskette. The H E A T source code is available in either U N I X or DOS versions. The code was run successfully on MIPS RC2030 and SUN 3/80 and 670-MP SPARC Temperature (°C) 0 10203040506070 10 20 , 0 A lo 1
downwardflowalongtilled 201-aqulfer(lO'°m.lec "t) "1 2 ~ 30~" iboundarytBmpemtureSO~C| ii at time. 15kyr(tranoienl)~
0.,
30
40
50 B
44
9
,t liT'" 80
60
'
Figure 3. A---conductive dissipation of thermal anomaly caused by instantaneous crustal thickening by thrusting; B~details (geotherms are shown at t,j, 0.4, 1, 2, 4, 10, 20, 40 Myr. and t, ).
10 ~
Figure 5. A-Evolution of temperature spike (50 C) introduced at t = 15 kyr (transient analysis) in tilted aquifer: B q e t a i l s (geotherms are shown at t~, 15, 30, 44, 59 kyr, 1.l, 4.3 Myr, and t, ).
Heat transport in geological formations
0.8 E 0.6
|
~0.4 r0.9 0.2
[i 0,I z' 0.5 z'=0.9
I
to5
I
-20
-10 0 10 20 Thermal Peclel number Figure 6. Comparison of equilibrium temperatures (t~) calculated by HEAT (curves), and exact temperatures (symbols) calculated analytically (after Spalding, 1972) for dimensionless distances (z') of 0.1, 0.5, and 0.9 for a wide range of Per ( - 2 0 to +20). workstations, and 386 and 486 MS-DOS personal computers during testing. A comparison of results obtained with HEAT and the exact (analytical) solution (as given by Spalding, 1972) is shown in Figure 6, in terms of the equilibrium temperature (time = o0) at dimensionless distances of 0.1, 0.5, and 0.9 as a function of fluid flow (Pe, from - 2 0 to +20). It can be seen that HEAT provides a reasonable estimate of temperature through this Pe, range. The input file contains the information on boundary and initial conditions, physical properties of the rock medium and porewater (see Table 1), grid size in the z direction and degree of accuracy required for the iterative scheme. Examples of input and output files are presented in Appendices 2 and 3, respectively. By default, a linear initial condition and constant boundary conditions (steady-state situation) are assumed at first. After entering or reading the input parameters, a prompt allowing the user to alter some initial condition temperatures will appear. To do this, simply enter the depth-steps and the new temperature values as couples of free-formatted numbers. Similarly, to change any value of the boundary conditions and thus perform a transient thermal analysis, three free-formatted numbers must be entered. These represent the first and last time-step of a given time interval, and the new boundary temperature value for that time interval. As important discontinuities in the boundary condition may strain the Crank-Nicolson scheme (Smith, 1985), it is advisable to keep these functions as smooth as possible. Before calculating the temperature profile at a new time-step n + 1, the program assigns it the temperature profile from the previous time-step, n. This is used as a seed for iteration, which then is performed
677
over and over again until the requested accuracy level between two consecutive iterations is obtained at every depth-step. Only then does the program proceed to the next time-step. For an accuracy of 0.0001 C, it was determined that less than l0 iterations were required in most test runs. After having reached the preselected accuracy level at any one time-step, the computer performs a stability check to determine whether the temperature profile is in a transient phase, or whether a thermal equilibrium has been reached. The stability check, which can be customized according to specific needs, currently consists of summing the differences in temperature for all the nodes between the current and the previous time-step. If that sum is smaller than a specified value (0.1~C for instance, see Appendix 1), the system is assumed to have reached steady-state, and a ' Y ' (yes) appears at the bottom of that column in response to the question "equil. '~ in the output file (Appendix 3). Otherwise, an ' N ' (no) is printed in the output file. If a transient system is under investigation, it is possible to reach thermal equilibrium, as defined, and thereafter to return to a nonequilibrium situation when one or both boundary conditions change. The size of the depth-step, Az, is specified in the input file (Appendix 2). Because the stability factor, r, is fixed at 0.4, changing A: will change At, the size of the time-step, via Equation (7). For each specific application, the user will need to change depth-step size, and time-step size, in order to match the particular scale of interest. The computer performs up to 9000 time-step calculations, but results of the preequilibrium thermal evolution are displayed only, together with the final, equilibrium temperature profile. The user then can decide to see details of the temperature profile evolution (TPE) by selecting the appropriate n-value (time-step). If equilibrium is not reached by the end of the 9000 time-steps, the user may decide to increase the depth-step size in order to increase the time-step size. Alternatively, the array declaration in the program can be altered, but increasing the array sizes will increase memory requirements and compilation and execution time. The PC-version of HEAT only performs up to 3000 time-step calculations to reduce memory requirements. The output file written by the program (Appendix 3) only presents fourteen time-steps in order to permit easy presentation in one single table. If the user elected to see certain parts of the TPE in detail, these also will appear in the output file. CONCLUSIONS HEAT is a simple, accurate, numerical approximation of the I-D conductive/advective heat transport equation. The method, based on the iterative solution of an implicit, finite difference, Crank-Nicolson algorithm, is implemented in a
678
P. DE CARITAT
F O R T R A N code presented herein. The p r o g r a m makes it possible to evaluate, in a straightforward m a n n e r , p e r t u r b a t i o n s to the thermal profile of geological formations, as caused by c o n d u c t i o n and advection by fluid flow. N o n l i n e a r initial condition and transient b o u n d a r y conditions can be specified.
Acknowledgments--This program was written, following a suggestion by Ron Spencer, while the author was the recipient of a University of Calgary Post-Doctoral Fellowship. lan Hutcheon provided moral support for this study as well as access to workstations. The original code was improved on the basis of discussions with Bernie Boudreau. Terry Gordon's assistance with compiling HEAT for the DOS environment is much appreciated. Joern Springer and an anonymous reviewer are acknowledged for their comments. REFERENCES Baker, J. C., and Caritat, P. de, 1992, Postdepositional history of the Permian sequence in the Denison Trough, eastern Australia: Am. Assoc. Petroleum Geologists Bull., v. 76, no. 8, p. 1224 1249. Bickle, M. J., and McKenzie, D., 1987, The transport of heat and matter by fluids during metamorphism: Contrib. Mineral. and Petrol., v. 95, no. 3, p. 384 392.
Caritat, P. de, and Baker, J. C., 1992, Oxygen-isotope evidence for upward, cross-formational porewater flow in a sedimentary basin near maximum burial: Sed. Geology, v. 78, no. 3/4, p. 155-164. Carahan, B., Luther, H. A., and Wilkes, J. O., 1969, Applied numerical methods: John Wiley & Sons, New York, 604 p. Carslaw, H. S., and Jaeger, J. C., 1959, Conduction of heat in solids (2nd ed.): Oxford Univ. Press, New York, 510 p. Croft, D. R., and Lilley, D. G., 1977, Heat transfer calculations using finite difference equations: Applied Science Publ. Ltd., London, 283 p. Gerald, C. F., and Wheatley, P. O., 1984, Applied numerical analysis (3d ed.): Addison-Wesley, Reading, Massachusetts, 579 p. Loosveld, R. J. H., 1989, The synchronism of crustal thickening and low-pressure facies metamorphism in the Mount lsa Inlier, Australia. 2. Fast convective thinning of mantle lithosphere during crustal thickening: Tectonophysics, v. 165, no. 1-4, p. 191-218. Smith, G. D., 1985, Numerical solution of partial differential equations: finite difference methods (3rd ed.): Oxford Univ. Press, Oxford, 337 p. Spalding, D. B., 1972, A novel finite difference formulation for differential expressions involving both first and second derivatives: Intern. Jour. for Num. Methods in Engineering, v. 4, no. 4, p. 551 559.
APPENDIX
1
HEAT Program Listing This program evaluates heat transport by advection and conduction, when one-dimensional fluid flow occurs in a rock sequence, using a finite difference, Crank-Nicolson approximation, solved here by iteration. COPYRIGHT date created:
PATRICE de CARITAT 1991,1992 910821 -- last revision:
921030
variable names: acc accuracy of the iteration scheme (C) cf specific heat capacity of fluid (J.kg-l.K-l) cs specific heat capacity of solid (J.kg-l.K-l) delt time grid size (s) delty time grid size (years) delz spatial grid size (m) dist distance (length, thickness) considered (m) gg geothermal gradient (C.m-l) kappa thermal diffusivity of solid [-thcond/(ros*cs)] (m2.s-l) lob lower boundary condition (surface temperature) (C) M cf*rof / cs*ros (dimensionless) mu M*W*delt/delz PeT thermal Peclet number (dimensionless) phi solid medium porosity (dimensionless) rof specific density of fluid (kg.m-3) ros specific density of solid (kg.m-3) rsf stability factor ('r' term = 0.4) sigma weighting parameter for f.d. scheme T(i,n) temperature (C) at time n & space i thcond thermal conductivity of solid (W.m-I.K-I) upb upper boundary condition (bottom temperature) (C) vf velocity of the fluid (m.s-l), negative if upward vf2 velocity of the fluid (m.a-l) w effective fluid flux [=vf*phi] (m.s-l) implicit integer parameter integer
none imaxd, nmax (imaxd=lO O, nmax=9 0 0 0 ) age(O:nmax),c,cmax,timestep(O:nmax),iflagl
Heat transport in geological formations & real & & & &
,imax,itn, nmx, nl,n2,n3,n4,ndet acc,cf, cs,delt,delty,delz,diff(nmax),dist,gg ,kappa,lob, upb, M, mu,PeT, phi,rage(0:nmax), rof, ros ,rsf, sigma,T(0:imaxd, 0:nmax) ,Tchk(0:imaxd, 0:nmax),thcond, Tlob, Tupb, Tinc ,vf,vf2,W
character*l answerl,answer2,answer3,answer4,equil(nmax) c c define constant used rsf - 0.4 iflagl = 0 c c open files open (unit-9, file-"heat.out",status-"unknown",err-lO000) rewind 9 c write (*,*) write (*,500) read (*,510) answer1 if ((answerl.eq.'y').or.(answerl.eq.'Y')) then open (unit-8, file-"heat.inp",status-"old",err-10000) rewind 8 goto 1 else open (unit-8,file-"heat.inp",status-"unknown",err-10000) rewind 8 endif c c enter input parameters from keyboard/screen write (*,520) read (*,*) lob write (%530) read (*,*) upb write (*,540) read (*,*) dist write (*,550) read (*,*) vf write (*,560) read (*,*) ros write (*,570) read (*,*) cs write (*,580) read (*,*) rof write (*,590) read (*,*) cf write (*,600) read (*,*) thcond write (*,610) read (*,*) phi write (*,620) read (*,*) delz write (*,630) read (*,*) acc c c write input parameters to "heat.inp" file (overwrites if exists} write (8,640) lob write (8,660) upb write (8,670) dist write (8,680) vf write (8,690) ros write (8,700) cs write (8,710) rof write (8,720) cf write (8,730} thcond write (8,740) phi write (8,750) delz write (8,760) acc goto 2 continue c reading parameters from input file ("heat.inp") & write to screen read (8,770) lob, upb,dist,vf, ros,cs, rof,cf, thcond, phi,delz, acc write (*,650) lob write (*,660) upb write (*,670) dist write (*,680) vf write (*,690) ros write (*,700) cs write (*,710) rof write (*,720) cf write (*,730) thcond
679
680
P. DE CARITAT
write write write
(*,740) phi (*,750) delz (*,760) acc
c 2 continue c c write input parameters to output file ("heat.out") wr~te (9,650) lob wrlte (9,660) upb wr~te (9,670) dist wr~te (9,680) vf wr~te (9,690) ros wrlte (9,700) cs wr~te (9,710) rof wrlte (9,720) cf write (9,730) thcond write (9,740) phi wrlte (9,750) delz write (9,760) acc c c calculate linear geothermal gradient gg=(upb-lob)/dist write (*,780) gg write (9,780) gg c c calculate M, the ratio of the products of specific c capacities by the specific densities of the fluid c to the solid (cf*rof/cs*ros) M=(cf*rof)/(cs*ros) write (*,790) M write (9,790) M c c calculate kappa, the thermal diffusivity of solid, c kappa = (thcond / (ros * cs)) kappa=(thcond/(ros*cs)) write (*,800) kappa write (9,800) kappa c c calculate W, the effective fluid flux, where c W = vf * phi W=vf*phi write (*,810) W write (9,810) W c c calculate fluid velocity in m/a vf2=vf*31557600 write (*,820) vf2 write (9,820) vf2 c c calculate thermal Peclet number PeT=W*dist/kappa write (*,830) PeT write (9,830) PeT c 3 continue c c calculate delt, the time grid size, where c delt = rsf * (delz)**2 / kappa
heat
where
delt=(rsf*(delz)**2)/kappa write (*,840) delt write (9,840) delt delty=delt/31557600 write (*,850) delty write (9,850) delty c c calculate mu, a reccurring ratio in the C-N scheme mu=M*W*delt/delz c c calculate sigma, the f.d. weighting parameter, which dictates c whether the advective term will be approximated using c upstream weighting f.d. scheme sigma-l.0 for PET>2 c midpoint weighting f.d. scheme sigma-0.5 for -2<=PET<=2 c downstream weighting f.d. scheme sigma-0.0 for PET<-2 sigma=0.5 if (PeT.gt.2.0) then sigma=l.0 endif if (PeT.it.-2.0) then sigma=0.0 endif imax=dist/delz
Heat transport in geological formations
if (imax.gt.imaxd) then write (*,860) imaxd write (*,1030) delz read (*,*) delz write (*,1040) write (9, 1040) goto 3 endi f c if ((imax*delz).ne.dist) write (*, 870) write (*,1030) delz read (*,*) delz write (*,1040) write (9,1040) goto 3 endi f
then
c c calculate initial condition do i=0, imax T (i, 0) =lob+ (i*delz*gg) end do c c calculate boundary conditions do n=l, nmax T (0, n) =lob T (imax, n) =upb end do c c allow for nonlinear initial condition write (*,2000) write (*,880) read (*,510) answer2 if ( (answer2.eq. 'n') .or. (answer2.eq. 'N') goto 5 write (*,890) write (9,890) do i=0,imax, 3 if (i+2.1e.imax) then write (*,900) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & , i+2, T (i+2, 0 ), (i+2) *del z write (9,900) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & ,i+2, T (i+2, 0) , (i+2) *delz elseif (i+l.eq.imax) then write (*,910) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz write (9,910) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz elseif (i.eq.imax) then write (*,920) i,T(i,0),i*delz write (9,920) i,T(i,0),i*delz endi f end do write (*,930) 4 continue read (*,*) i,Tinc if ((i.ne.-l) .and. (Tinc.ne.0)) then T (i, 0) =Tinc goto 4 endi f write (*,940) write (9,940) do i=0, imax, 3 if (i+2.1e.imax) then write (*,950) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & , i+2, T (i+2, 0) , (i+2) *delz write (9,950) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & ,i+2,T(i+2, 0), (i+2) *delz elseif (i+1.eq. lmax) then write (*,960) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz write (9,960) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz elseif (i.eq. imax) then write (*,970) i,T(i,0),i*delz write (9,970) i,T(i,0),i*delz endif end do c 5 continue c write (*, 2000) c c allow for transient boundary conditions write (*, 980) read (*,510) answer3 CAGEO19:5~E
681
682
P. DE CARITAT
6
c 7 c 8
c 9 c
if ((answer3.eq. 'n') .or. (answer3.eq.'N')) goto 9 write (*,990) lob, delty continue read (*,*) nl,n2,Tlob if ((nl.eq.0).and. (n2.eq.0) .and. (Tlob.eq.0)) then goto 7 else write (9,1000) nl,n2,Tlob do n=nl, n2 T (0, n) =Tlob end do endi f goto 6 continue write (*,1010) upb, delty continue read (*,*) n3,n4,Tupb if ((n3.eq.0) .and. (n4.eq.0) .and. (Tupb.eq.0)) goto 9 else write (9,1000) n3,n4,Tupb do n=n3, n4 T (imax, n) =Tupb end do endi f goto 8
then
continue
write (*,2000) c c Crank-Nicolson algorithm cmax=0 c c start the timestep loop in 'n' .................................. do n=0, nmax-i diff (n+l) =0.0 c=0 c c preset the temperature in column 'n+l' do i=l, imax-i T (i, n+l) =T (i, n) end do 10 continue c c set Tchk equivalent to previous value of T do i=l, imax-I Tchk (i, n+l) =T (i, n+l) end do c c calculate the new temperature value using the C-N method do i=l, imax-i T(i,n+l) = (T(i+l,n+l)*(rsf/2-mu*(1-sigma)/2)+ & T (i-l, n+l) * (rsf/2+mu*sigma/2)+ & T (i+l, n) * (rsf/2-mu* (1-sigma)/2) + & T(i,n) * (l-rsf-mu* (-l+2*sigma)/2) + & T (i-l, n) * (rsf/2+mu*sigma/2)) / & (l+rsf+mu* (-l+2*sigma)/2) end do c check that C-N temperature is within 'acc' of preset temperature c if not, reset the preset value as the C-N value, and recalculate. c record number of iteration performed do i=l, imax-i if (abs(Tchk(i,n+l)-T(i,n+l)).ge.acc) then C=C+I if (c.gt.cmax) then cmax=c endi f goto 10 endl f end do c c check the difference in temperature between timestep 'n+l' c and (previous) timetstep 'n'. Sum up the diff at any timestep do i=1, imax-1 cliff (n+l) =diff (n+l) +abs (T (i, n+l ) -T (i, n) ) end do if (diff(n+l).lt.O.1) then equil (n+l) ='Y' else
Heat transport in geological formations equil(n+l)='N' endif c c close the timestep loop in vnv .................................. end do c if (equil(nmax).eq.'N') then write (*, i020) read (*,510) answer4 if ((answer4.eq.'y') .or. (answer4.eq. 'Y')) then write (*,1030) delz read (*,*) delz write (*,1040) write (9,1040) goto 3 endi f endi f c write (*,1050) cmax write (9,1050) cmax c c write out the output in 14 columns (screen & output file) c nmx=nmax do n-nmax, 130, -i if ((equil(n) .eq.'N') .and. (equil(n+l) .eq.'Y')) then write (*,1060) n+l,int((n*delty/1000)+0.5) write (9,1060) n+l,int((n*delty/1000)+0.5) nmx-n+13 goto Ii else nmx=130 endif end do c ii continue c if (nmx.gt.nmax) then D_rax~D/sax endl f c do n~0t nmx timestep (n) =n age (n) =int ((n*delty/1000) +0.5) end do c if (age(nmx).ge.10000) then iflagl=l do n=0, nmx rage (n) =n*delty/1000000 end do endi f c write (*,1070) (timestep(n),n=0,nmx, int((nmx/13)+0.5)) write (9,1070) (timestep(n),n=0,nmx, int((nmx/13)+0.5)) if (iflagl.eq.0) then write (*,1080) (age(n),n=0,nmx, int((nmx/13)+0.5)) write (9,1080) (age(n),n=0,nmx, int((nmx/13)+0.5)) elseif (iflagl.eq.l) then write (*,1090) (rage(n),n=0,nmx, int((nmx/13)+0.5)) write (9,1090) (rage(n),n=0,nmx, int((nmx/13)+0.5)) endif c write (*,*) write (9,*) write (*,1100) write (9,1100) do i=0,imax write (*,1110) i, (int ( (T (i, n) +0.5) ), n=0, nmx, int ((nmx/13) +0.5) ), & int ((i*delz) +0.5) write (9,1110) i, (int ( (T (i,n) +0.5) ), n=0, nmx, int ((nmx/13) +0.5) ), & int ((i*delz) +0.5) end do (equil (n) ,n=int ((nmx/13) +0.5) ,nmx, write (*,1120) & int ((nmx/13) +0.5) ) write (9,1120) (equil (n) ,n=int ((nmx/13) +0.5) ,nmx, &
int ((nmx/13) +0.5) )
c write write
(*,2000) (9,2000)
683
684
P. DE CARITAT
c display parts of TPE tableau in detail write (*,1130) write (9,1130) c continue 12 c write (*,1140) read (*,*) ndet if (ndet.eq.-l) goto 9000 if (ndet.gt.nmax-13) then write (*,1150) nmax-12 goto 12 endif write (9,1140) write (9,*) ndet
(every timestep)
write (*,1070) (timestep(n),n=ndet,ndet+13) write (9,1070) (timestep(n),n=ndet,ndet+13) if (iflagl.eq.0) then write (*,I080) (age(n),n=ndet,ndet+13) write (9,1080) (age(n),n=ndet,ndet+13) elseif (iflagl.eq.l) then write (*,1090) (rage(n),n=ndet,ndet+13) write (9,1090) (rage(n),n=ndet,ndet+13) endif write (*,*) write (9,*) write (*,Ii00) write (9,1100) do i~0, imax write (*,1110) i, (int((T(i,n)+0.5)),n=ndet,ndet+13), & int ((i*delz) +0.5) write (9,1110) i, (int((T(i,n)+0.5)),n=ndet,ndet+13), & int((i*delz)+0.5) end do write (*,1120) (equil(n),n=ndet+l,ndet+13) write (9,1120) (equil(n),n=ndet+l,ndet+13) c write (*, 2000) write (9, 2000) goto 12 c 9000 c
continue
stop c c errors i0000 write write
(*,*) (9,*)
' <<< error occurred -- check code or input file >>>' ' <<< error occurred -- check code or input file >>>'
c close close
(unit=8) (unit=9)
c c formats 500 format (Ix, 79('*'),/,/,8x,' Program HEAT',I0x, &'Copyright 1991, 1992 Patrice de Caritat',/,/,Ix, 79('*'), &/,/," Do you want to read input from ", &"existing 'heat.inp' file (y/n)? ",$) format (al) 510 520 format (/,' Enter the appropriate parameters as prompted',/, & ' Lower boundary temperature (C): ',8x, ',$) ',8x, ',$) 530 format (' Upper boundary temperature (C): ',8x, ',$) 540 format (' Distance of interest (m): 550 format (' Fluid velocity (m.s-1): ',8x, ',$) 560 format (' Solid specific density (kg.m-3): ',8x, ',$) ',8x, ',$) 570 format (' Solid heat capacity (J.kg-l.K-l): ',8x, ',$) 580 format (' Fluid specific density (kg.m-3): ',8x, ',$) 590 format (' Fluid heat capacity (J.kg-l.K-1): 600 format (' Solid thermal conductivity (W.m-l.K-1):',Sx,' ',$) ',Sx,' ',$) 610 format (' Porosity: ',8x,' ',$) 620 format (' Grid size (space dimension) (m): ',8x, ' ',$) 630 format (' Accuracy (C) : ', 5x, f8.2) 640 format (' Lower boundary temperature (C): 650 format (/,' *** INPUT PARAMETERS',/, ',5x, f8.2) & ' Lower boundary temperature (C): ', 5x, f8.2) 660 format (' Upper boundary temperature (C): ', 5x, f8.2) 670 format (' Distance of interest (m): ', 8x, e10.3) 680 format (' Fluid velocity (re.s-l) : ', 5x, f8.2) 690 format (' Solid specific density (kg.m-3): ',5x, f8.2) 700 format (' Solid heat capacity (J.kg-l.K-1):
Heat transport in geological formations 710 720 730 740 750 760 770 780
format format format format format format format format
(' (' (' (' (' ('
685
Fluid specific density (kg.m-3): ',Sx, f8.2) Fluid heat capacity (J. kg-l.K-l): ',5x, f8.2) Solid thermal conductivity (W.m-l.K-l):',5x, f8.2) Porosity: ',5x, f8.2) Grid size (space dimension) (m): ',5x, f8.2) Accuracy (C): ',5x, f12.6)
(3(45x, f8.2,/),48x, elO.3,/,7(45x, f8.2,/),45x, f12.6)
(Ix,79(' '),/,/,' *** CALCULATED PARAMETERS',/, ' Linea~ geothermal gradient (C.m-l): ',6x, f10.4) format (' M (dimensionless): ',8x, f9.5) format (' Kappa (m2.s-l): ',10x, e9.3) format (' w (m.s-l): ',10x,e9.3) format (' Fluid velocity (m.a-l): ',10x,e9.3) format (' Thermal Peclet number (dimensionless):',10x, e9.3) format (' Delta time (s) : ',10x, e9.3) format (' Delta time (years): ',6x, f9.3) format (/,' Ratio dist/delz must be less than ',i5) format (/,' Ratio dist/delz must be an integer!!') format (/,' Do you want to modify initial (linear) ', &'condition (y/n)? ',$) format (/,' Present initial condition reads:',/, &3(Ix,' i',2('_'),'T(C)',5(' '),'z(m)', ('_'),3x)) format (3(i3,1x, f6.1,1x, fS.2,3x)) format (2(i3,1x, f6.1,1x, fS.2,3x)) format (i3,1x, f6.1,1x, fS.2) format (/,' Enter i and new T (-1,0 to finish)') format (/,' Modified initial condition reads:',/, &3(lx,' i',2(' ') 'T(C)',5(' '),'z(m)', (' '),3x)) format (2(i3,~x, f6.1,1x, f8.2,3x),i3,1x, f6.1,1x, f8.2) format (i3,1x, f6.1,1x, f8.2,3x, i3,1x, f6.1,1x, f8.2) format (i3,1x, f6.1,1x, f8.2) format (/,' Do you want to modify the boundary ', &'conditions (transient case) (y/n)? ',$) format (/,' Present lower boundary condition reads:', & f8.2,' C',/, & ' Current time-step is',f12.2,' years',/, &' Enter nl, n2 and new lower boundary T (0,0,0 to finish)') format (/,' Modified lower boundary condition:',/, &' from n = ',i3,' to n = ',i3,' T(iob) = ',f7.2) format (/,' Present upper boundary condition reads:',fS.2, & ' C',/, & ' Current time-step is',fl2.2,' years',/, &' Enter nl, n2 and new upper boundary T (0,0,0 to finish)') format (/,' Thermal equilibrium not attained',/, &' Do you want to increase delta z (y/n)? ',$) format (' Old delta z value (m):',22x, f8.2,/, & ' Enter new delta z value (m): ',15x,' ',$) format (/,' *** RECALCULATED PARAMETERS') format (/,/,' Max. # iterations performed =',Ix, i3) format (' Thermal equilibrium reached at timestep ',i4, &' (time ~ ',i6,' kyrs)') format (/,' n:',14(i4,1x)) format (' age:',14(i4,1x),' kyrs') format (' age:',14(f4.1,1x),' Myrs') format (' iJ',16x, &' TEMPERATURE PROFILE EVOLUTION (TPE)',I5x,'J z') format (ix,i3,'r',13(i4,1x),i4,'J',i5) format (/,' equil?',6x,13(al,4x)) format (/,' To see part of TPE tableau in detail') format (' Enter n-value of leftmost column (-i to exit): ',$) format (' Warning: n-value must be less than ',i5) format (ix,79(' '))
&
790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950
960 970 980 990
1000 1010
1020 1030 1040 1050 1060 1070 1080 1090 1100 1110 1120 1130 1140 1150 2000 c
end
APPENDIX 2 Input File Lower boundary temperature (C): Upper boundary temperature (C): Distance of interest (m): Fluid velocity (m.s-1): Solid specific density (kg.m-3): Solid heat capacity (J.kg-l.K-1): Fluid specific density (kg.m-3): Fluid heat capacity (J.kg-l.K-1): Solid thermal conductivity (W.m-l.K-1): Porosity: Grid size (space dimension) (m): Accuracy (C):
10.00 300.00 10000.00 -0.100E-08 270O.00 1300.00 i010.00
3800.00 3.00 0.I0 250.00 0.000100
686
P. DE CARITAT
APPENDIX 3 Output File *** INPUT PARAMETERS Lower boundary temperature (C): Upper boundary temperature (C): Distance of interest (m): Fluid velocity (m.s-1): Solid specific density (kg.m-3): Solid heat capacity (J.kg-l.K-1): Fluid specific density (kg.m-3): Fluid heat capacity (J.kg-l.K-1): Solid thermal conductivity (W.m-l.K-1): Porosity: Grid size (space dimension) (m): Accuracy (C):
10.00 300.00 i0000.00 -0.100E-08 2700.00 1300.00 1010.00 3800.00 3.00 0.10 250.00 0.000100
*** CALCULATED PARAMETERS Linear geothermal gradient (C.m-l): M (dimensionless): Kappa (m2.s-l): W (m.s-l): Fluid velocity (m.a-l): Thermal Peclet number (dimensionless): Delta time (s): Delta time (years):
0.0290 1.09345 0.855E-06 -.100E-09 -.316E-01 -.II7E+01 0.293E+II 926.877
Present initial condition reads: T(C) z(m) i T(C) z(m) --i-- 17.3 250.00 o lO.-0-~ o.o~ 3 31.8 750.00 4 39.0 i000.00 6 53.5 1500.00 7 60.8 1750.00 9 75.3 2250.00 i0 82.5 2500.00 13 104.3 3250.00 12 97.0 3000.00 16 126.0 4000.00 15 118.8 3750.00 19 147.8 4750.00 18 140.5 4500.00 21 162.3 5250.00 22 169.5 5500.00 24 184.0 6000.00 25 191.3 6250.00 28 213.0 7000.00 27 205.8 6750.00 30 227.5 7500.00 31 234.8 7750.00 33 249.3 8250.00 34 256.5 8500.00 37 278.3 9250.00 36 271.0 9000.00 39 292.8 9750.00 40 300.0 i0000.00
i 2 5 8 ii 14 17 20 23 26 29 32 35 38
T(C) 24.5 46.3 68.0 89.8 111.5 133.3 155.0 176.8 198.5 220.3 242.0 263.8 285.5
z(m) 500.00 1250.00 2000.00 2750.00 3500.00 4250.00 5000.00 5750.00 6500.00 7250.00 8000.00 8750.00 9500.00
Modified initial condition reads: i T(C) z(m) i T(C) z(m) 0 I0.0 0.0~ 1 17.3 250.00 3 31.8 750.00 4 39.0 I000.00 6 53.5 1500.00 7 60.8 1750.00 9 75.3 2250.00 I0 82.5 2500.00 12 97.0 3000.00 13 104.3 3250.00 15 118.8 3750.00 16 126.0 4000.00 19 200.0 4750 00 18 140.5 4500.00 21 200.0 5250.00 22 169.5 5500 00 25 191.3 6250 00 24 184.0 6000.00 27 205.8 6750.00 28 213.0 7000 00 30 227.5 7500.00 31 234.8 7750 00 33 249.3 8250.00 34 256.5 8500 00 36 271.0 9000.00 37 278.3 9250 00 39 292.8 9750.00 40 300.0 i0000.00
i 2 5 8 ii 14 17 20 23 26 29 32 35 38
T(C) 24.5 46.3 68.0 89.8 111.5 133.3 200.0 176.8 198.5 220.3 242.0 263.8 285.5
z(m) 500.00 1250.00 2000.00 2750.00 3500.00 4250.00 5000.00 5750.00 6500.00 7250.00 8000.00 8750.00 9500.00
i
Max. # iterations performed = 7 Thermal equilibrium reached at t i ~ s t e p
1271
(time ~
1177 kyrs)
n: age:
686 636
784 727
980 1078 1176 1274 908 999 1090 1181
i 0 1 2 3 4 5 6 7 8 9 I0
0 0
98 91
196 182
294 273
392 363
490 454
588 545
i0 17 25 32 39 46 54 61 68 75 83
i0 19 28 37 46 54 62 70 79 87 94
I0 20 30 39 49 58 67 75 84 92 100
TEMPERATURE PROFILE i0 I0 i0 i0 21 21 21 22 31 32 33 33 41 42 43 44 51 52 54 55 60 62 64 65 70 72 74 75 79 81 83 85 88 90 93 94 96 99 102 104 105 108 111 113
EVOLUTION i0 I0 22 22 33 34 45 45 56 56 66 67 76 77 86 87 96 97 105 106 114 116
882 818 I0 22 34 45 57 67 78 88 98 107 117
(TPE) i0 22 34 46 57 68 78 89 98 108 117
i0 22 34 46 57 68 79 89 99 109 118
10 22 34 46 57 68 79 89 99 109 118
I i01 221 351 461 581 691 791 901 100i 1091 1191
kyrs z 0 250 500 750 1000 1250 1500 1750 2000 2250 2500
~
~
~
~
~
~
~
~
~
~
o
ooo
o
o o ~ o o ~
~
~
0
~
o
~
o
~ ~
~ ~
~
~
0
~~
O
~
~
~
~
J
o
~
~
~
~
O
~
~
J
~
~
~
~
~
~
o
o
~o~
~
~
~
~
~
~
~ ~
O
~
~
~
~
~
O
~
~ 0
v
0
O
~
~
~~
~
oo
. . . .
~
~ o ~
~ ~o
0 ~ 0 ~ 0 ~ 0 ~ 0 ~ 0 ~ 0 ~ ~0~0~ ~0~0~ ~0~0 0~0~ oooooooooooooogooooo~ooooogoooo~ooooo~
~
0
~
o ~.o.~. .
O
o ~ ~ o ~ ~ o ~ ° ~ ~ ~ o ~ ~ o
~ o
o ~ ~ ~ ~ ~ ~ o ~ ~ o ~ ~ . ~ ~o~~
o ~
~
~
0
o
~
~
~
~ ~
O
o
0
o o o
rt
rt 0
0 0
0
I...-
I-.-, rP
~0
Z
Z
Z
Z
Z
Z
m:
~Z
Z
Z
Z
Z
H~ .J
o~o~o~o~o~o~o~o~g~o~o~g~o~o~ o o o o oooog~ o o o o o o o o o o o o o o
0 0 ~
0
5~
0
o
"0 0
:Z
0 0 ~ 0 0
0 0 0 0 0 0
~o~o~o~o~o~o~o~o~o~
~ ~ o ~ o ~ ~
~
0 0 0
~
o
O 0 0 N
0
0
~
oo
0
t,ol~ co~ r.
z
.
.
.
.
.
.
000
°
~
0 0 0 0 0
~
.
~
.
0 0 0 0 ~ 0
~
o~o~~o~~o~~~~~
.
~
.
~
~
.
~
~
~
.
.
o
.
~
.
.
~
o
~
.
0 0 0 0 0 ~ 0 0 0 0 0 0 0 0 0 0 0
.
~
o
~ O 0 0 0 0 0 N
~
~ ~ ~
.
~
c
~
~
..,,
°°
g
So
o
I"1
z
t..,I-. ,.0
>
.,..a
Heat transport in geological formations 231 241 251 261 271 284 291 301 311 321 331 34t 35f 36i 371 38J 391 401
188 194 200 206 212 218 224 231 238 245 252 259 266 273 280 287 293 300
equil?
188 194 200 206 212 218 225 231 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 218 225 231 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 218 225 232 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 219 225 232 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 219 225 232 239 245 252 259 266 273 280 287 294 300
188 194 200 206 213 219 225 232 239 246 253 260 266 273 280 287 294 300
N
N
N
N
N
N
213 219 226 232 239 246 253 260 267 274 280 287 294 300
188 194 201 207 213 219 226 232 239 246 253 260 267 274 280 287 294 300
188 194 201 207 213 219 226 233 239 246 253 260 267 274 280 287 294 300
188 194 201 207 213 220 226 233 239 246 253 260 267 274 280 287 294 300
188 195 201 207 213 220 226 233 239 246 253 260 267 274 280 287 294 300
N
N
N
N
N
N
48 44
49 45
50 46
51 47
52 48
53 49
TEMPERATURE PROFILE EVOLUTION (TPE) i0 i0 i0 i0 I0 I0 i0 I0 18 18 18 18 18 18 18 18 26 26 27 27 27 27 27 27 34 34 34 35 35 35 35 35 42 42 42 42 42 43 43 43 50 50 50 50 50 50 50 51 58 58 58 58 58 58 58 58 65 66 66 66 66 66 66 66 73 73 73 74 74 74 74 74 81 81 81 81 82 82 82 82 89 89 89 89 90 90 90 90 97 97 97 97 98 98 98 98 105 105 105 106 106 106 106 106 113 113 114 114 114 114 114 114 122 122 122 122 122 122 122 122 130 130 130 130 130 130 130 130 138 138 138 138 138 138 138 138 146 146 146 146 146 146 146 146 153 153 153 153 153 153 153 153 161 161 161 161 161 161 161 161 168 168 168 168 168 168 168 168 175 175 175 175 175 175 175 175 182 182 182 182 182 182 182 182 188 188 188 188 188 188 188 188 195 195 195 195 195 195 195 195 201 201 201 201 201 201 201 201 207 207 207 208 208 208 208 208 214 214 214 214 214 214 214 214 220 220 220 220 221 221 221 221 227 227 227 227 227 227 227 227 233 233 233 234 234 234 234 234 240 240 240 240 240 240 241 241 247 247 247 247 247 247 247 247 253 254 254 254 254 254 254 254 260 260 260 261 261 261 261 261 267 267 267 267 267 267 268 268 274 274 274 274 274 274 274 274 281 281 281 281 281 281 281 281 287 287 287 287 287 287 287 287 294 294 294 294 294 294 294 294 300 300 300 300 300 300 300 300
i0 1B 27 35 43 51 58 66 74 82 90 98 106 114 122 130 138 146 153 161 168 175 182 189 195 202 208 214 221 227 234 241 247 254 261 268 274 281 287 294 300
i0 18 27 35 43 51 59 66 74 82 90 98 106 114 122 130 138 146 153 161 168 175 182 189 195 202 208 215 221 228 234 241 248 254 261 268 274 281 287 294 300
N
N
Enter n-value o f l e f t m o s t 4O
column
n: age:
40 37
41 38
42 39
44 41
il OI I) 21 31 4t 5f 61 71 81 9[ 10r III 121 13i 141 151 161 17t 181 19f 201 211 221 231 241 25[ 261 27[ 281 29] 301 311 32i 331 341 351 361 37[ 381 39i 40[
i0 18 26 34 42 50 57 65 73 81 89 97 105 113 121 129 138 146 153 161 168 175 182 188 195 201 207 213 220 226 233 239 246 253 260 267 274 280 287 294 300
I0 18 26 34 42 50 57 65 73 81 89 97 105 113 121 130 138 146 153 161 168 175 182 188 195 201 207 213 220 226 233 240 246 253 260 267 274 281 287 294 300
I0 18 26 34 42 50 58 65 73 81 89 97 105 113 121 130 138 146 153 161 168 175 182 188 195 201 207 214 220 226 233 240 247 253 260 267 274 281 287 294 300
N
N
equil?
43 40
N
N
188 194 200
188 194 200
206
207
213 219 225 232 239 246 253 260 267 273 280 287 294 300 N
689 188 5750 1 9 5 ' 6000 2 0 1 6250 207 6500 2131 6750 220~ 7000 226 7250 233 7500 240 7750 246 8000 253 8250 260 8500 267 8750 274 9000 281 9250 287 9500 294 9750 300 I0000
(il to exit):
45 42
N
46 43
N
47 44
N
N
N
N
kyrs
I z i01 0 191 250 271 500 351 750 43J 1000 511 1250 591 1500 671 1750 741 2000 82( 2250 901 2500 98# 2750 106 3000 114 3250 122 3500 130 3750 138 4000 146 4250 153 4500 161 4750 168 5000 175 5250 182 5500 189 575O 195 6000 202 6250 208 6500 215 6750 221 7000 228 7250 234 7500 241 7750 248 8000 254 8250 261 8500 268 8750 274 9000 281 9250 288 9500 294 9750 300 i0000 N
0098-3004/93 $6.00 + 0.00 Copyright c 1993 Pergamon Press Ltd
Printed in Great Britain. All rights reserved
HEAT: A FORTRAN COMPUTER PROGRAM FOR CALCULATING 1-D CONDUCTIVE AND ADVECTIVE HEAT TRANSPORT IN GEOLOGICAL FORMATIONS PATRICE DE CARITAT Department of Geology and Geophysics, University of Calgary, Alberta, Canada T2N IN4
(Received 6 April 1992; revised 9 July 1992; accepted 24 November 1992)
Abstract--A simple, accurate, numerical approximation of the one-dimensional equation of heat transport by conduction and advection is presented. The method is based on the iterative solution of an implicit, finite difference, Crank-Nicolson algorithm, featuring alternating differencing direction as a function of the thermal Peclet number, and is implemented in a documented FORTRAN code. The program allows users to evaluate how porewater flow may affect the thermal profile of permeable rock formations, under steady-state or unsteady-state boundary conditions, and with linear or nonlinear initial conditions. Several possible applications illustrate the potential usefulness of the program in first order, evaluative investigations. Key Wordr: Heat transport, Advection, Conduction, Finite difference, Numerical modeling.
INTRODUCTION Heat transport in saturated porous media, such as rock formations below the watertable, takes place as a result of conduction and, if the contained porewater is moving sufficiently fast, advection. The applications of heat transport to geological processes are numerous and important, for instance in metamorphic, ore deposit, hydrogeoiogical, thermal, and diagenetic studies. Geoscientists who are not inclined numerically, however, may be discouraged to tackle this issue because of the apparent complexity of the mathematical and numerical manipulations involved. The purpose of this communication is to provide these geoscientists with a simple, one-dimensional (l-D) computer simulation program, which will perform an accurate evaluation of conductive and advectire heat transport in porous media. HEAT TRANSPORT
IN P O R O U S
MEDIA
In a saturated permeable medium, such as a porous or fractured sedimentary rock, heat is transported by conduction ("diffusion") through the solid framework and within the fluid phase. If the contained fluid is moving through this permeable medium, advection ("convection") also contributes to heat transport. The relative importance of conductive versus advective heat transport is a function of the thermal diffusivity of the rock and fluid flow rate (Bickle and McKenzie, 1987). The principles of heat transport are presented only briefly herein, and a simple and accurate numerical solution is presented. Transport in the vertical direc-
tion will be specifically considered, although this solution applies to heat transport in any direction. As it would not be possible nor appropriate to present, within the context of this communication, a thorough review of the vast literature pertaining to heat transport, the reader is referred to Carslaw and Jaeger (1959) for an introduction to the physics of heat transport, and to Bickle and McKenzie (1987) for a presentation of geological applications, among many others. The conductive/advective heat transport equation, in its simplest form and in one-dimension, is given by: c?T O,t
k ?,2T_ pfQ -
= ~-
OT vf0 ~-,
p~c~ cz"
p~c~
cz
(I)
where T is the temperature (K), t the time (sec), k the thermal conductivity of the solid medium (W m 1K-J), Ps and pf are the solid and fluid specific density (kg m -3), respectively, cs and cf the solid and fluid specific heat capacity (J k g - ~ K - ~), respectively, vf is the fluid velocity (m sec-~), negative for upward flow, 4) the porosity of the solid medium, and z the direction in which heat transport is considered, here vertical (see Table 1 for a listing of symbols. suggested values and units). The following substitutions: ~," -
(2)
pfcf , P~Cs
(3)
W ~ ~t,~9,
(4)
M
673
k p~c~ ' ~ -
674
P. DE C^R,r^r Table 1. Symbols, suggested values and units (after Loosveld, 1989) symbol
meaning
value (range)
unit
ct" ca m
~ i I i c Imatc ~ c i t y of fluid ~aceific heat capacity of solid time-step
3800-4000 1000-1300 afleulamt
Liz
dcpth-st¢~
I~¢-~
J.kg-l.K-1 J.kg-l.K-1 see m
¢
r,om~ty
length of interest (for Pet) thermal diffusivity of solid
- 0.0-o.25 pre-selected
--
h ~'= k
- 10-7-10"6
m2.scc-1
2.0-3.0
W.m-1.K-1
Ps Cs
k M
thermal conductivity of solid
= Pfct
Ps cs Wh Pe, It"
r =
I¢ At
A z2
pf Ps
T t
v/ W = vf#
z
d2T
1, calculated
-
--
thermal Peclet ntmaber
calculated
--
stability factor
0.4 (fixed)
--
specific density of fluid specific denfity of solid temperature time fluid velocity effective fluid flux direction considemt
1000-1010 2700-3300 --prc-sel~"ted calcula~ --
kg.m-3 kg.m-3 K see
where x is termed the thermal diffusivity of the solid medium (m 2 see-t), and 14," the effective fluid flux (m sec- ~), result in the following simplified expression of Equation (I): dT
m
OT (5)
Ot = ~-ff~z2 - M W - ~ z "
This is the partial differential equation describing 1-D conductive/advective heat transport in porous media. It states that the rate of change of temperature at any point is equivalent to the sum of conductive heat transport (first term of the right-hand side) and advective heat transport (second term of the righthand side). In this form, it represents a simplification of the real-world situation, as it neglects internal heat generation, thermal expansion, compressibility, the temperature-dependency of heat capacity and thermal conductivity, etc. Furthermore, it is applied here to a homogeneous (nonlayered) 1-D medium. A time-dependent solution for this type of equation cannot be obtained analytically, except for the simplest of situations. However, a solution for the heat transport equation in the general transient (unsteadystate) situation can be approached satisfactorily using the Crank-Nicolson approximation (see Carnahan,
m.sec -I
m.sec -I m
Luther, and Wilkes, 1969; Croft and Liiley, 1977; Gerald and Wheatley, 1984). A computer program, HEAT, was designed to implement this procedure. In order to preserve the generality of this computer program, applicability to both slow and fast flowing systems was considered. In particular, provision was made to alternate the direction of the differencing scheme for the advection term as a function of the thermal Peclet number Pe, (Pe t = W h / x , where h is the length-scale of interest), as suggested by Spalding (1972). Thus, the final Crank-Nicolson approximation of Equation (5) is given by: i [TT÷ ~ _ T~] = x [T;':~ ~- 2T7 +~ + T;'_+~' at
~
"
+ T~+l - 2 T 7 + T'~ " Az 2 !I - MW
az 2
]
"+' ,,+l -o')T~+,-(20-I)T i +
4 (1 - tT)T'i' + 1-- (2o"2_Az-I)T7 + t;T'i'-I ],
Pet < -2 -2 < Pet < 2 2 < Pet
~valuc 0.0 0.5 1.0
(6)
..I
where a is the parameter that alternates the direction of the differencing scheme for the advection term as given in Table 2. In this formulation, the space and
Table 2. Conditions governing the direction of differencing scheme condition
~T,,+, i-l
y~
din~clion downmtam weighting midpoint weighting upstream weighting
Heat transport in geological formations time grids are represented by i (subscript) and n (superscript), respectively, and Az and At are the depth- (or space-) and time-steps, respectively. We now introduce r, the stability factor (or Fourier number): xAt r - Az"'
(7)
and #:
MWAt
U= - - , Az
675 t
.,,
at
IAz
l
IN
T~n+1 • ~7"~nlahi-1 qF
(8)
dhl qw
+ TT+l + r,
5+
-5
! -r
+T,'.' i
s'~l %J
'T.. n T.,n+1 ~i+1 ~i+1
and, rearranging Equation (6), we obtain the temperature at the ith depth-step and n + lth time-step:
r?+,=) ,+, 1_5-5' - 0 ) +
IF
En Tn+1
-5(2a
[-;
+
- 1) l+r+
]
(2~r-l) .
(9) Selecting r = 0.4 results in an unconditionally stable solution (Loosveld, 1989). The solution of Equation (5) is satisfactorily approximated by iteratively solving for the implicit, finite difference, Crank-Nicolson algorithm [Eq. (9)]. Figure 1 illustrates the grid or node system used here. Time increases from left to right, and depth increases downward. The t - z continuum under consideration is discretized using At and Az as time- and depth-steps, respectively. As can be seen in Figure i, the solution for the temperature at any node depends on the temperature at closest neighboring nodes at the same time-step, and on the temperature at the same and neighboring depth-step at the previous time-step. Such a configuration, which is referred to as implicit (Carnahan, Luther, and Wilkes, 1969), has the advantage of yielding stable results for any value of r, although small values for r ( ,~ 1) give more accurate results (Croft and Lilley, 1977; Gerald and Wheatley, 1984). Selected boundary conditions can be chosen to be the temperature at two fixed levels at the extremities of the space coordinate of interest (Dirichlet-type boundary conditions). A practical selection of boundary conditions would include the temperature at the surface (z = 0) and the temperature at the deepest level of interest (z = Zm,~), for instance the basement of a sedimentary basin, or at the base of the lithosphere. The selection of boundary conditions is a function of the scale of interest and the problem at hand. However, it is desirable to select the high temperature boundary condition to be as great as possible (greater depth, higher temperature), so as to avoid interference of this boundary within the (shallower) region of interest. In HEAT, the value of the
Figure 1. Discretization of time space (t :) continuum for finite difference approximation scheme. Approximation of temperature at any point (O) depends on temperature at five neighboring nodes (O) in t-: space, this being an implicit numerical scheme. two boundary conditions can be specified as a function of time, so that transient analysis can be performed. The initial condition may be taken to be the linear temperature gradient (conductive geothermal gradient) between these two levels. The program allows for a nonequilibrium (nonlinear) initial temperature profile to be specified. Figure 2 illustrates, in nondimensionalized terms, the effect of various conditions of upward, vertical porewater flow upon the equilibrium thermal profile (t = t~ ) in sedimentary basins, with respect to the linear conductive geotherm. APPLICATIONS To illustrate the purpose of this program, a brief description of a problem to which it was applied
0.0
I
,
J
0.2
,¢
0.4
06
o
08
,
?o,=+a5
• Pet=-17.6 1.0
'
0.0
' 0.2
'
' 0.4
"~'~ '
' 0.6
Temperature
'
' 0.8
'
-~1.0
Figure 2. Nondimensionalized representation of perturbations to conductive geotherm (11) as a result of upward, vertical porewater flow characterized by decreasing Pe, values (upward fluid flow requires negative sign).
676
P, DE CARITAT
follows. Baker and Caritat (1992) contended that ankerite cementation of the Aldebaran Sandstone, the main gas reservoir in the Denison Trough, eastern Australia, occurred as a result of vertical, upward, cross-formational porewater flow. This flow system developed as a thick sequence of overpressured continental mudrocks and coals underlying the Aldebaran Sandstone underwent overpressure relaxation, as a result of intensifying tectonic stress during the Late Triassic, immediately prior to basin inversion. The computer program H E A T was used to determine how much heat would have been advected under these circumstances as a function of time, and this constraint was used subsequently to model the temperature-sensitive fractionation of oxygen isotopes between the porewater and the precipitating ankerite (Caritat and Baker, 1992). Selected boundary conditions were the temperature at the surface and at the base of the sedimentary sequence, and initial condition was taken to be the conductive geotherm (25-3OC/km). Vertical flow velocities were estimated by considering the pressure differential and the hydraulic conductivity across the aquitard capping the overpressured section. Under the particular conditions of the Denison Trough, it was determined that ambient temperatures in the Aldebaran Sandstone were increased because of the upward advection of hotter water by 5-10' C relative to conductive temperatures of 120-140 C (depending on the location in the Denison Trough) during akerite cementation (Caritat and Baker, 1992). Other applications may include the conductive dissipation of a thermal anomaly caused by overthrusting (Fig. 3), the conductive/advective relaxation of a temperature spike introduced by sill emplacement or discharge of hot (hydrothermal) water in a permeable (p.e., fractured) horizion (Fig. 4), and the
0
200
400
Temperature 600 800 0
(°C) 200 400
600
2o
800 B
10 (conduction on~y) Imtantaneouscrusp,I ~ 20L ~zthl- c2koni n gbythrul C i n km ~~g~
o.,
o
0
100
Temperature (°C) 200 300 120 140 160 180 200
:up~va flow (.i0.o m.u¢,-1)~ 2 ~nltlal"hot"layer (200"C1,,t ~ 4.75 to 5.25 krn dq~h 3 ~z.2sOm \
~ 4
5 6
7
8 9 10
...........
Figure 4. A--Dissipation of thermal spike (200 C) in initial geotherm, mostly as result of conduction. This is followed by bowing-upward of geotherm through advection by upward fluid flow; B--details (geotherms are shown at 6 , 1, 2, 4, 6, 10, 50kyr, and t~).
propagation of a transient thermal anomaly in a tilted aquifer (Fig. 5). THE COMPUTER PROGRAM HEAT, a F O R T R A N computer code written tosolve Equation (9), is presented in Appendix 1. The code is annotated abundantly to allow the user to follow its logic, and is available from the author upon submission of a 3½in, floppy diskette. The H E A T source code is available in either U N I X or DOS versions. The code was run successfully on MIPS RC2030 and SUN 3/80 and 670-MP SPARC Temperature (°C) 0 10203040506070 10 20 , 0 A lo 1
downwardflowalongtilled 201-aqulfer(lO'°m.lec "t) "1 2 ~ 30~" iboundarytBmpemtureSO~C| ii at time. 15kyr(tranoienl)~
0.,
30
40
50 B
44
9
,t liT'" 80
60
'
Figure 3. A---conductive dissipation of thermal anomaly caused by instantaneous crustal thickening by thrusting; B~details (geotherms are shown at t,j, 0.4, 1, 2, 4, 10, 20, 40 Myr. and t, ).
10 ~
Figure 5. A-Evolution of temperature spike (50 C) introduced at t = 15 kyr (transient analysis) in tilted aquifer: B q e t a i l s (geotherms are shown at t~, 15, 30, 44, 59 kyr, 1.l, 4.3 Myr, and t, ).
Heat transport in geological formations
0.8 E 0.6
|
~0.4 r0.9 0.2
[i 0,I z' 0.5 z'=0.9
I
to5
I
-20
-10 0 10 20 Thermal Peclel number Figure 6. Comparison of equilibrium temperatures (t~) calculated by HEAT (curves), and exact temperatures (symbols) calculated analytically (after Spalding, 1972) for dimensionless distances (z') of 0.1, 0.5, and 0.9 for a wide range of Per ( - 2 0 to +20). workstations, and 386 and 486 MS-DOS personal computers during testing. A comparison of results obtained with HEAT and the exact (analytical) solution (as given by Spalding, 1972) is shown in Figure 6, in terms of the equilibrium temperature (time = o0) at dimensionless distances of 0.1, 0.5, and 0.9 as a function of fluid flow (Pe, from - 2 0 to +20). It can be seen that HEAT provides a reasonable estimate of temperature through this Pe, range. The input file contains the information on boundary and initial conditions, physical properties of the rock medium and porewater (see Table 1), grid size in the z direction and degree of accuracy required for the iterative scheme. Examples of input and output files are presented in Appendices 2 and 3, respectively. By default, a linear initial condition and constant boundary conditions (steady-state situation) are assumed at first. After entering or reading the input parameters, a prompt allowing the user to alter some initial condition temperatures will appear. To do this, simply enter the depth-steps and the new temperature values as couples of free-formatted numbers. Similarly, to change any value of the boundary conditions and thus perform a transient thermal analysis, three free-formatted numbers must be entered. These represent the first and last time-step of a given time interval, and the new boundary temperature value for that time interval. As important discontinuities in the boundary condition may strain the Crank-Nicolson scheme (Smith, 1985), it is advisable to keep these functions as smooth as possible. Before calculating the temperature profile at a new time-step n + 1, the program assigns it the temperature profile from the previous time-step, n. This is used as a seed for iteration, which then is performed
677
over and over again until the requested accuracy level between two consecutive iterations is obtained at every depth-step. Only then does the program proceed to the next time-step. For an accuracy of 0.0001 C, it was determined that less than l0 iterations were required in most test runs. After having reached the preselected accuracy level at any one time-step, the computer performs a stability check to determine whether the temperature profile is in a transient phase, or whether a thermal equilibrium has been reached. The stability check, which can be customized according to specific needs, currently consists of summing the differences in temperature for all the nodes between the current and the previous time-step. If that sum is smaller than a specified value (0.1~C for instance, see Appendix 1), the system is assumed to have reached steady-state, and a ' Y ' (yes) appears at the bottom of that column in response to the question "equil. '~ in the output file (Appendix 3). Otherwise, an ' N ' (no) is printed in the output file. If a transient system is under investigation, it is possible to reach thermal equilibrium, as defined, and thereafter to return to a nonequilibrium situation when one or both boundary conditions change. The size of the depth-step, Az, is specified in the input file (Appendix 2). Because the stability factor, r, is fixed at 0.4, changing A: will change At, the size of the time-step, via Equation (7). For each specific application, the user will need to change depth-step size, and time-step size, in order to match the particular scale of interest. The computer performs up to 9000 time-step calculations, but results of the preequilibrium thermal evolution are displayed only, together with the final, equilibrium temperature profile. The user then can decide to see details of the temperature profile evolution (TPE) by selecting the appropriate n-value (time-step). If equilibrium is not reached by the end of the 9000 time-steps, the user may decide to increase the depth-step size in order to increase the time-step size. Alternatively, the array declaration in the program can be altered, but increasing the array sizes will increase memory requirements and compilation and execution time. The PC-version of HEAT only performs up to 3000 time-step calculations to reduce memory requirements. The output file written by the program (Appendix 3) only presents fourteen time-steps in order to permit easy presentation in one single table. If the user elected to see certain parts of the TPE in detail, these also will appear in the output file. CONCLUSIONS HEAT is a simple, accurate, numerical approximation of the I-D conductive/advective heat transport equation. The method, based on the iterative solution of an implicit, finite difference, Crank-Nicolson algorithm, is implemented in a
678
P. DE CARITAT
F O R T R A N code presented herein. The p r o g r a m makes it possible to evaluate, in a straightforward m a n n e r , p e r t u r b a t i o n s to the thermal profile of geological formations, as caused by c o n d u c t i o n and advection by fluid flow. N o n l i n e a r initial condition and transient b o u n d a r y conditions can be specified.
Acknowledgments--This program was written, following a suggestion by Ron Spencer, while the author was the recipient of a University of Calgary Post-Doctoral Fellowship. lan Hutcheon provided moral support for this study as well as access to workstations. The original code was improved on the basis of discussions with Bernie Boudreau. Terry Gordon's assistance with compiling HEAT for the DOS environment is much appreciated. Joern Springer and an anonymous reviewer are acknowledged for their comments. REFERENCES Baker, J. C., and Caritat, P. de, 1992, Postdepositional history of the Permian sequence in the Denison Trough, eastern Australia: Am. Assoc. Petroleum Geologists Bull., v. 76, no. 8, p. 1224 1249. Bickle, M. J., and McKenzie, D., 1987, The transport of heat and matter by fluids during metamorphism: Contrib. Mineral. and Petrol., v. 95, no. 3, p. 384 392.
Caritat, P. de, and Baker, J. C., 1992, Oxygen-isotope evidence for upward, cross-formational porewater flow in a sedimentary basin near maximum burial: Sed. Geology, v. 78, no. 3/4, p. 155-164. Carahan, B., Luther, H. A., and Wilkes, J. O., 1969, Applied numerical methods: John Wiley & Sons, New York, 604 p. Carslaw, H. S., and Jaeger, J. C., 1959, Conduction of heat in solids (2nd ed.): Oxford Univ. Press, New York, 510 p. Croft, D. R., and Lilley, D. G., 1977, Heat transfer calculations using finite difference equations: Applied Science Publ. Ltd., London, 283 p. Gerald, C. F., and Wheatley, P. O., 1984, Applied numerical analysis (3d ed.): Addison-Wesley, Reading, Massachusetts, 579 p. Loosveld, R. J. H., 1989, The synchronism of crustal thickening and low-pressure facies metamorphism in the Mount lsa Inlier, Australia. 2. Fast convective thinning of mantle lithosphere during crustal thickening: Tectonophysics, v. 165, no. 1-4, p. 191-218. Smith, G. D., 1985, Numerical solution of partial differential equations: finite difference methods (3rd ed.): Oxford Univ. Press, Oxford, 337 p. Spalding, D. B., 1972, A novel finite difference formulation for differential expressions involving both first and second derivatives: Intern. Jour. for Num. Methods in Engineering, v. 4, no. 4, p. 551 559.
APPENDIX
1
HEAT Program Listing This program evaluates heat transport by advection and conduction, when one-dimensional fluid flow occurs in a rock sequence, using a finite difference, Crank-Nicolson approximation, solved here by iteration. COPYRIGHT date created:
PATRICE de CARITAT 1991,1992 910821 -- last revision:
921030
variable names: acc accuracy of the iteration scheme (C) cf specific heat capacity of fluid (J.kg-l.K-l) cs specific heat capacity of solid (J.kg-l.K-l) delt time grid size (s) delty time grid size (years) delz spatial grid size (m) dist distance (length, thickness) considered (m) gg geothermal gradient (C.m-l) kappa thermal diffusivity of solid [-thcond/(ros*cs)] (m2.s-l) lob lower boundary condition (surface temperature) (C) M cf*rof / cs*ros (dimensionless) mu M*W*delt/delz PeT thermal Peclet number (dimensionless) phi solid medium porosity (dimensionless) rof specific density of fluid (kg.m-3) ros specific density of solid (kg.m-3) rsf stability factor ('r' term = 0.4) sigma weighting parameter for f.d. scheme T(i,n) temperature (C) at time n & space i thcond thermal conductivity of solid (W.m-I.K-I) upb upper boundary condition (bottom temperature) (C) vf velocity of the fluid (m.s-l), negative if upward vf2 velocity of the fluid (m.a-l) w effective fluid flux [=vf*phi] (m.s-l) implicit integer parameter integer
none imaxd, nmax (imaxd=lO O, nmax=9 0 0 0 ) age(O:nmax),c,cmax,timestep(O:nmax),iflagl
Heat transport in geological formations & real & & & &
,imax,itn, nmx, nl,n2,n3,n4,ndet acc,cf, cs,delt,delty,delz,diff(nmax),dist,gg ,kappa,lob, upb, M, mu,PeT, phi,rage(0:nmax), rof, ros ,rsf, sigma,T(0:imaxd, 0:nmax) ,Tchk(0:imaxd, 0:nmax),thcond, Tlob, Tupb, Tinc ,vf,vf2,W
character*l answerl,answer2,answer3,answer4,equil(nmax) c c define constant used rsf - 0.4 iflagl = 0 c c open files open (unit-9, file-"heat.out",status-"unknown",err-lO000) rewind 9 c write (*,*) write (*,500) read (*,510) answer1 if ((answerl.eq.'y').or.(answerl.eq.'Y')) then open (unit-8, file-"heat.inp",status-"old",err-10000) rewind 8 goto 1 else open (unit-8,file-"heat.inp",status-"unknown",err-10000) rewind 8 endif c c enter input parameters from keyboard/screen write (*,520) read (*,*) lob write (%530) read (*,*) upb write (*,540) read (*,*) dist write (*,550) read (*,*) vf write (*,560) read (*,*) ros write (*,570) read (*,*) cs write (*,580) read (*,*) rof write (*,590) read (*,*) cf write (*,600) read (*,*) thcond write (*,610) read (*,*) phi write (*,620) read (*,*) delz write (*,630) read (*,*) acc c c write input parameters to "heat.inp" file (overwrites if exists} write (8,640) lob write (8,660) upb write (8,670) dist write (8,680) vf write (8,690) ros write (8,700) cs write (8,710) rof write (8,720) cf write (8,730} thcond write (8,740) phi write (8,750) delz write (8,760) acc goto 2 continue c reading parameters from input file ("heat.inp") & write to screen read (8,770) lob, upb,dist,vf, ros,cs, rof,cf, thcond, phi,delz, acc write (*,650) lob write (*,660) upb write (*,670) dist write (*,680) vf write (*,690) ros write (*,700) cs write (*,710) rof write (*,720) cf write (*,730) thcond
679
680
P. DE CARITAT
write write write
(*,740) phi (*,750) delz (*,760) acc
c 2 continue c c write input parameters to output file ("heat.out") wr~te (9,650) lob wrlte (9,660) upb wr~te (9,670) dist wr~te (9,680) vf wr~te (9,690) ros wrlte (9,700) cs wr~te (9,710) rof wrlte (9,720) cf write (9,730) thcond write (9,740) phi wrlte (9,750) delz write (9,760) acc c c calculate linear geothermal gradient gg=(upb-lob)/dist write (*,780) gg write (9,780) gg c c calculate M, the ratio of the products of specific c capacities by the specific densities of the fluid c to the solid (cf*rof/cs*ros) M=(cf*rof)/(cs*ros) write (*,790) M write (9,790) M c c calculate kappa, the thermal diffusivity of solid, c kappa = (thcond / (ros * cs)) kappa=(thcond/(ros*cs)) write (*,800) kappa write (9,800) kappa c c calculate W, the effective fluid flux, where c W = vf * phi W=vf*phi write (*,810) W write (9,810) W c c calculate fluid velocity in m/a vf2=vf*31557600 write (*,820) vf2 write (9,820) vf2 c c calculate thermal Peclet number PeT=W*dist/kappa write (*,830) PeT write (9,830) PeT c 3 continue c c calculate delt, the time grid size, where c delt = rsf * (delz)**2 / kappa
heat
where
delt=(rsf*(delz)**2)/kappa write (*,840) delt write (9,840) delt delty=delt/31557600 write (*,850) delty write (9,850) delty c c calculate mu, a reccurring ratio in the C-N scheme mu=M*W*delt/delz c c calculate sigma, the f.d. weighting parameter, which dictates c whether the advective term will be approximated using c upstream weighting f.d. scheme sigma-l.0 for PET>2 c midpoint weighting f.d. scheme sigma-0.5 for -2<=PET<=2 c downstream weighting f.d. scheme sigma-0.0 for PET<-2 sigma=0.5 if (PeT.gt.2.0) then sigma=l.0 endif if (PeT.it.-2.0) then sigma=0.0 endif imax=dist/delz
Heat transport in geological formations
if (imax.gt.imaxd) then write (*,860) imaxd write (*,1030) delz read (*,*) delz write (*,1040) write (9, 1040) goto 3 endi f c if ((imax*delz).ne.dist) write (*, 870) write (*,1030) delz read (*,*) delz write (*,1040) write (9,1040) goto 3 endi f
then
c c calculate initial condition do i=0, imax T (i, 0) =lob+ (i*delz*gg) end do c c calculate boundary conditions do n=l, nmax T (0, n) =lob T (imax, n) =upb end do c c allow for nonlinear initial condition write (*,2000) write (*,880) read (*,510) answer2 if ( (answer2.eq. 'n') .or. (answer2.eq. 'N') goto 5 write (*,890) write (9,890) do i=0,imax, 3 if (i+2.1e.imax) then write (*,900) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & , i+2, T (i+2, 0 ), (i+2) *del z write (9,900) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & ,i+2, T (i+2, 0) , (i+2) *delz elseif (i+l.eq.imax) then write (*,910) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz write (9,910) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz elseif (i.eq.imax) then write (*,920) i,T(i,0),i*delz write (9,920) i,T(i,0),i*delz endi f end do write (*,930) 4 continue read (*,*) i,Tinc if ((i.ne.-l) .and. (Tinc.ne.0)) then T (i, 0) =Tinc goto 4 endi f write (*,940) write (9,940) do i=0, imax, 3 if (i+2.1e.imax) then write (*,950) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & , i+2, T (i+2, 0) , (i+2) *delz write (9,950) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz & ,i+2,T(i+2, 0), (i+2) *delz elseif (i+1.eq. lmax) then write (*,960) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz write (9,960) i,T(i,0),i*delz,i+l,T(i+l,0), (i+l)*delz elseif (i.eq. imax) then write (*,970) i,T(i,0),i*delz write (9,970) i,T(i,0),i*delz endif end do c 5 continue c write (*, 2000) c c allow for transient boundary conditions write (*, 980) read (*,510) answer3 CAGEO19:5~E
681
682
P. DE CARITAT
6
c 7 c 8
c 9 c
if ((answer3.eq. 'n') .or. (answer3.eq.'N')) goto 9 write (*,990) lob, delty continue read (*,*) nl,n2,Tlob if ((nl.eq.0).and. (n2.eq.0) .and. (Tlob.eq.0)) then goto 7 else write (9,1000) nl,n2,Tlob do n=nl, n2 T (0, n) =Tlob end do endi f goto 6 continue write (*,1010) upb, delty continue read (*,*) n3,n4,Tupb if ((n3.eq.0) .and. (n4.eq.0) .and. (Tupb.eq.0)) goto 9 else write (9,1000) n3,n4,Tupb do n=n3, n4 T (imax, n) =Tupb end do endi f goto 8
then
continue
write (*,2000) c c Crank-Nicolson algorithm cmax=0 c c start the timestep loop in 'n' .................................. do n=0, nmax-i diff (n+l) =0.0 c=0 c c preset the temperature in column 'n+l' do i=l, imax-i T (i, n+l) =T (i, n) end do 10 continue c c set Tchk equivalent to previous value of T do i=l, imax-I Tchk (i, n+l) =T (i, n+l) end do c c calculate the new temperature value using the C-N method do i=l, imax-i T(i,n+l) = (T(i+l,n+l)*(rsf/2-mu*(1-sigma)/2)+ & T (i-l, n+l) * (rsf/2+mu*sigma/2)+ & T (i+l, n) * (rsf/2-mu* (1-sigma)/2) + & T(i,n) * (l-rsf-mu* (-l+2*sigma)/2) + & T (i-l, n) * (rsf/2+mu*sigma/2)) / & (l+rsf+mu* (-l+2*sigma)/2) end do c check that C-N temperature is within 'acc' of preset temperature c if not, reset the preset value as the C-N value, and recalculate. c record number of iteration performed do i=l, imax-i if (abs(Tchk(i,n+l)-T(i,n+l)).ge.acc) then C=C+I if (c.gt.cmax) then cmax=c endi f goto 10 endl f end do c c check the difference in temperature between timestep 'n+l' c and (previous) timetstep 'n'. Sum up the diff at any timestep do i=1, imax-1 cliff (n+l) =diff (n+l) +abs (T (i, n+l ) -T (i, n) ) end do if (diff(n+l).lt.O.1) then equil (n+l) ='Y' else
Heat transport in geological formations equil(n+l)='N' endif c c close the timestep loop in vnv .................................. end do c if (equil(nmax).eq.'N') then write (*, i020) read (*,510) answer4 if ((answer4.eq.'y') .or. (answer4.eq. 'Y')) then write (*,1030) delz read (*,*) delz write (*,1040) write (9,1040) goto 3 endi f endi f c write (*,1050) cmax write (9,1050) cmax c c write out the output in 14 columns (screen & output file) c nmx=nmax do n-nmax, 130, -i if ((equil(n) .eq.'N') .and. (equil(n+l) .eq.'Y')) then write (*,1060) n+l,int((n*delty/1000)+0.5) write (9,1060) n+l,int((n*delty/1000)+0.5) nmx-n+13 goto Ii else nmx=130 endif end do c ii continue c if (nmx.gt.nmax) then D_rax~D/sax endl f c do n~0t nmx timestep (n) =n age (n) =int ((n*delty/1000) +0.5) end do c if (age(nmx).ge.10000) then iflagl=l do n=0, nmx rage (n) =n*delty/1000000 end do endi f c write (*,1070) (timestep(n),n=0,nmx, int((nmx/13)+0.5)) write (9,1070) (timestep(n),n=0,nmx, int((nmx/13)+0.5)) if (iflagl.eq.0) then write (*,1080) (age(n),n=0,nmx, int((nmx/13)+0.5)) write (9,1080) (age(n),n=0,nmx, int((nmx/13)+0.5)) elseif (iflagl.eq.l) then write (*,1090) (rage(n),n=0,nmx, int((nmx/13)+0.5)) write (9,1090) (rage(n),n=0,nmx, int((nmx/13)+0.5)) endif c write (*,*) write (9,*) write (*,1100) write (9,1100) do i=0,imax write (*,1110) i, (int ( (T (i, n) +0.5) ), n=0, nmx, int ((nmx/13) +0.5) ), & int ((i*delz) +0.5) write (9,1110) i, (int ( (T (i,n) +0.5) ), n=0, nmx, int ((nmx/13) +0.5) ), & int ((i*delz) +0.5) end do (equil (n) ,n=int ((nmx/13) +0.5) ,nmx, write (*,1120) & int ((nmx/13) +0.5) ) write (9,1120) (equil (n) ,n=int ((nmx/13) +0.5) ,nmx, &
int ((nmx/13) +0.5) )
c write write
(*,2000) (9,2000)
683
684
P. DE CARITAT
c display parts of TPE tableau in detail write (*,1130) write (9,1130) c continue 12 c write (*,1140) read (*,*) ndet if (ndet.eq.-l) goto 9000 if (ndet.gt.nmax-13) then write (*,1150) nmax-12 goto 12 endif write (9,1140) write (9,*) ndet
(every timestep)
write (*,1070) (timestep(n),n=ndet,ndet+13) write (9,1070) (timestep(n),n=ndet,ndet+13) if (iflagl.eq.0) then write (*,I080) (age(n),n=ndet,ndet+13) write (9,1080) (age(n),n=ndet,ndet+13) elseif (iflagl.eq.l) then write (*,1090) (rage(n),n=ndet,ndet+13) write (9,1090) (rage(n),n=ndet,ndet+13) endif write (*,*) write (9,*) write (*,Ii00) write (9,1100) do i~0, imax write (*,1110) i, (int((T(i,n)+0.5)),n=ndet,ndet+13), & int ((i*delz) +0.5) write (9,1110) i, (int((T(i,n)+0.5)),n=ndet,ndet+13), & int((i*delz)+0.5) end do write (*,1120) (equil(n),n=ndet+l,ndet+13) write (9,1120) (equil(n),n=ndet+l,ndet+13) c write (*, 2000) write (9, 2000) goto 12 c 9000 c
continue
stop c c errors i0000 write write
(*,*) (9,*)
' <<< error occurred -- check code or input file >>>' ' <<< error occurred -- check code or input file >>>'
c close close
(unit=8) (unit=9)
c c formats 500 format (Ix, 79('*'),/,/,8x,' Program HEAT',I0x, &'Copyright 1991, 1992 Patrice de Caritat',/,/,Ix, 79('*'), &/,/," Do you want to read input from ", &"existing 'heat.inp' file (y/n)? ",$) format (al) 510 520 format (/,' Enter the appropriate parameters as prompted',/, & ' Lower boundary temperature (C): ',8x, ',$) ',8x, ',$) 530 format (' Upper boundary temperature (C): ',8x, ',$) 540 format (' Distance of interest (m): 550 format (' Fluid velocity (m.s-1): ',8x, ',$) 560 format (' Solid specific density (kg.m-3): ',8x, ',$) ',8x, ',$) 570 format (' Solid heat capacity (J.kg-l.K-l): ',8x, ',$) 580 format (' Fluid specific density (kg.m-3): ',8x, ',$) 590 format (' Fluid heat capacity (J.kg-l.K-1): 600 format (' Solid thermal conductivity (W.m-l.K-1):',Sx,' ',$) ',Sx,' ',$) 610 format (' Porosity: ',8x,' ',$) 620 format (' Grid size (space dimension) (m): ',8x, ' ',$) 630 format (' Accuracy (C) : ', 5x, f8.2) 640 format (' Lower boundary temperature (C): 650 format (/,' *** INPUT PARAMETERS',/, ',5x, f8.2) & ' Lower boundary temperature (C): ', 5x, f8.2) 660 format (' Upper boundary temperature (C): ', 5x, f8.2) 670 format (' Distance of interest (m): ', 8x, e10.3) 680 format (' Fluid velocity (re.s-l) : ', 5x, f8.2) 690 format (' Solid specific density (kg.m-3): ',5x, f8.2) 700 format (' Solid heat capacity (J.kg-l.K-1):
Heat transport in geological formations 710 720 730 740 750 760 770 780
format format format format format format format format
(' (' (' (' (' ('
685
Fluid specific density (kg.m-3): ',Sx, f8.2) Fluid heat capacity (J. kg-l.K-l): ',5x, f8.2) Solid thermal conductivity (W.m-l.K-l):',5x, f8.2) Porosity: ',5x, f8.2) Grid size (space dimension) (m): ',5x, f8.2) Accuracy (C): ',5x, f12.6)
(3(45x, f8.2,/),48x, elO.3,/,7(45x, f8.2,/),45x, f12.6)
(Ix,79(' '),/,/,' *** CALCULATED PARAMETERS',/, ' Linea~ geothermal gradient (C.m-l): ',6x, f10.4) format (' M (dimensionless): ',8x, f9.5) format (' Kappa (m2.s-l): ',10x, e9.3) format (' w (m.s-l): ',10x,e9.3) format (' Fluid velocity (m.a-l): ',10x,e9.3) format (' Thermal Peclet number (dimensionless):',10x, e9.3) format (' Delta time (s) : ',10x, e9.3) format (' Delta time (years): ',6x, f9.3) format (/,' Ratio dist/delz must be less than ',i5) format (/,' Ratio dist/delz must be an integer!!') format (/,' Do you want to modify initial (linear) ', &'condition (y/n)? ',$) format (/,' Present initial condition reads:',/, &3(Ix,' i',2('_'),'T(C)',5(' '),'z(m)', ('_'),3x)) format (3(i3,1x, f6.1,1x, fS.2,3x)) format (2(i3,1x, f6.1,1x, fS.2,3x)) format (i3,1x, f6.1,1x, fS.2) format (/,' Enter i and new T (-1,0 to finish)') format (/,' Modified initial condition reads:',/, &3(lx,' i',2(' ') 'T(C)',5(' '),'z(m)', (' '),3x)) format (2(i3,~x, f6.1,1x, f8.2,3x),i3,1x, f6.1,1x, f8.2) format (i3,1x, f6.1,1x, f8.2,3x, i3,1x, f6.1,1x, f8.2) format (i3,1x, f6.1,1x, f8.2) format (/,' Do you want to modify the boundary ', &'conditions (transient case) (y/n)? ',$) format (/,' Present lower boundary condition reads:', & f8.2,' C',/, & ' Current time-step is',f12.2,' years',/, &' Enter nl, n2 and new lower boundary T (0,0,0 to finish)') format (/,' Modified lower boundary condition:',/, &' from n = ',i3,' to n = ',i3,' T(iob) = ',f7.2) format (/,' Present upper boundary condition reads:',fS.2, & ' C',/, & ' Current time-step is',fl2.2,' years',/, &' Enter nl, n2 and new upper boundary T (0,0,0 to finish)') format (/,' Thermal equilibrium not attained',/, &' Do you want to increase delta z (y/n)? ',$) format (' Old delta z value (m):',22x, f8.2,/, & ' Enter new delta z value (m): ',15x,' ',$) format (/,' *** RECALCULATED PARAMETERS') format (/,/,' Max. # iterations performed =',Ix, i3) format (' Thermal equilibrium reached at timestep ',i4, &' (time ~ ',i6,' kyrs)') format (/,' n:',14(i4,1x)) format (' age:',14(i4,1x),' kyrs') format (' age:',14(f4.1,1x),' Myrs') format (' iJ',16x, &' TEMPERATURE PROFILE EVOLUTION (TPE)',I5x,'J z') format (ix,i3,'r',13(i4,1x),i4,'J',i5) format (/,' equil?',6x,13(al,4x)) format (/,' To see part of TPE tableau in detail') format (' Enter n-value of leftmost column (-i to exit): ',$) format (' Warning: n-value must be less than ',i5) format (ix,79(' '))
&
790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950
960 970 980 990
1000 1010
1020 1030 1040 1050 1060 1070 1080 1090 1100 1110 1120 1130 1140 1150 2000 c
end
APPENDIX 2 Input File Lower boundary temperature (C): Upper boundary temperature (C): Distance of interest (m): Fluid velocity (m.s-1): Solid specific density (kg.m-3): Solid heat capacity (J.kg-l.K-1): Fluid specific density (kg.m-3): Fluid heat capacity (J.kg-l.K-1): Solid thermal conductivity (W.m-l.K-1): Porosity: Grid size (space dimension) (m): Accuracy (C):
10.00 300.00 10000.00 -0.100E-08 270O.00 1300.00 i010.00
3800.00 3.00 0.I0 250.00 0.000100
686
P. DE CARITAT
APPENDIX 3 Output File *** INPUT PARAMETERS Lower boundary temperature (C): Upper boundary temperature (C): Distance of interest (m): Fluid velocity (m.s-1): Solid specific density (kg.m-3): Solid heat capacity (J.kg-l.K-1): Fluid specific density (kg.m-3): Fluid heat capacity (J.kg-l.K-1): Solid thermal conductivity (W.m-l.K-1): Porosity: Grid size (space dimension) (m): Accuracy (C):
10.00 300.00 i0000.00 -0.100E-08 2700.00 1300.00 1010.00 3800.00 3.00 0.10 250.00 0.000100
*** CALCULATED PARAMETERS Linear geothermal gradient (C.m-l): M (dimensionless): Kappa (m2.s-l): W (m.s-l): Fluid velocity (m.a-l): Thermal Peclet number (dimensionless): Delta time (s): Delta time (years):
0.0290 1.09345 0.855E-06 -.100E-09 -.316E-01 -.II7E+01 0.293E+II 926.877
Present initial condition reads: T(C) z(m) i T(C) z(m) --i-- 17.3 250.00 o lO.-0-~ o.o~ 3 31.8 750.00 4 39.0 i000.00 6 53.5 1500.00 7 60.8 1750.00 9 75.3 2250.00 i0 82.5 2500.00 13 104.3 3250.00 12 97.0 3000.00 16 126.0 4000.00 15 118.8 3750.00 19 147.8 4750.00 18 140.5 4500.00 21 162.3 5250.00 22 169.5 5500.00 24 184.0 6000.00 25 191.3 6250.00 28 213.0 7000.00 27 205.8 6750.00 30 227.5 7500.00 31 234.8 7750.00 33 249.3 8250.00 34 256.5 8500.00 37 278.3 9250.00 36 271.0 9000.00 39 292.8 9750.00 40 300.0 i0000.00
i 2 5 8 ii 14 17 20 23 26 29 32 35 38
T(C) 24.5 46.3 68.0 89.8 111.5 133.3 155.0 176.8 198.5 220.3 242.0 263.8 285.5
z(m) 500.00 1250.00 2000.00 2750.00 3500.00 4250.00 5000.00 5750.00 6500.00 7250.00 8000.00 8750.00 9500.00
Modified initial condition reads: i T(C) z(m) i T(C) z(m) 0 I0.0 0.0~ 1 17.3 250.00 3 31.8 750.00 4 39.0 I000.00 6 53.5 1500.00 7 60.8 1750.00 9 75.3 2250.00 I0 82.5 2500.00 12 97.0 3000.00 13 104.3 3250.00 15 118.8 3750.00 16 126.0 4000.00 19 200.0 4750 00 18 140.5 4500.00 21 200.0 5250.00 22 169.5 5500 00 25 191.3 6250 00 24 184.0 6000.00 27 205.8 6750.00 28 213.0 7000 00 30 227.5 7500.00 31 234.8 7750 00 33 249.3 8250.00 34 256.5 8500 00 36 271.0 9000.00 37 278.3 9250 00 39 292.8 9750.00 40 300.0 i0000.00
i 2 5 8 ii 14 17 20 23 26 29 32 35 38
T(C) 24.5 46.3 68.0 89.8 111.5 133.3 200.0 176.8 198.5 220.3 242.0 263.8 285.5
z(m) 500.00 1250.00 2000.00 2750.00 3500.00 4250.00 5000.00 5750.00 6500.00 7250.00 8000.00 8750.00 9500.00
i
Max. # iterations performed = 7 Thermal equilibrium reached at t i ~ s t e p
1271
(time ~
1177 kyrs)
n: age:
686 636
784 727
980 1078 1176 1274 908 999 1090 1181
i 0 1 2 3 4 5 6 7 8 9 I0
0 0
98 91
196 182
294 273
392 363
490 454
588 545
i0 17 25 32 39 46 54 61 68 75 83
i0 19 28 37 46 54 62 70 79 87 94
I0 20 30 39 49 58 67 75 84 92 100
TEMPERATURE PROFILE i0 I0 i0 i0 21 21 21 22 31 32 33 33 41 42 43 44 51 52 54 55 60 62 64 65 70 72 74 75 79 81 83 85 88 90 93 94 96 99 102 104 105 108 111 113
EVOLUTION i0 I0 22 22 33 34 45 45 56 56 66 67 76 77 86 87 96 97 105 106 114 116
882 818 I0 22 34 45 57 67 78 88 98 107 117
(TPE) i0 22 34 46 57 68 78 89 98 108 117
i0 22 34 46 57 68 79 89 99 109 118
10 22 34 46 57 68 79 89 99 109 118
I i01 221 351 461 581 691 791 901 100i 1091 1191
kyrs z 0 250 500 750 1000 1250 1500 1750 2000 2250 2500
~
~
~
~
~
~
~
~
~
~
o
ooo
o
o o ~ o o ~
~
~
0
~
o
~
o
~ ~
~ ~
~
~
0
~~
O
~
~
~
~
J
o
~
~
~
~
O
~
~
J
~
~
~
~
~
~
o
o
~o~
~
~
~
~
~
~
~ ~
O
~
~
~
~
~
O
~
~ 0
v
0
O
~
~
~~
~
oo
. . . .
~
~ o ~
~ ~o
0 ~ 0 ~ 0 ~ 0 ~ 0 ~ 0 ~ 0 ~ ~0~0~ ~0~0~ ~0~0 0~0~ oooooooooooooogooooo~ooooogoooo~ooooo~
~
0
~
o ~.o.~. .
O
o ~ ~ o ~ ~ o ~ ° ~ ~ ~ o ~ ~ o
~ o
o ~ ~ ~ ~ ~ ~ o ~ ~ o ~ ~ . ~ ~o~~
o ~
~
~
0
o
~
~
~
~ ~
O
o
0
o o o
rt
rt 0
0 0
0
I...-
I-.-, rP
~0
Z
Z
Z
Z
Z
Z
m:
~Z
Z
Z
Z
Z
H~ .J
o~o~o~o~o~o~o~o~g~o~o~g~o~o~ o o o o oooog~ o o o o o o o o o o o o o o
0 0 ~
0
5~
0
o
"0 0
:Z
0 0 ~ 0 0
0 0 0 0 0 0
~o~o~o~o~o~o~o~o~o~
~ ~ o ~ o ~ ~
~
0 0 0
~
o
O 0 0 N
0
0
~
oo
0
t,ol~ co~ r.
z
.
.
.
.
.
.
000
°
~
0 0 0 0 0
~
.
~
.
0 0 0 0 ~ 0
~
o~o~~o~~o~~~~~
.
~
.
~
~
.
~
~
~
.
.
o
.
~
.
.
~
o
~
.
0 0 0 0 0 ~ 0 0 0 0 0 0 0 0 0 0 0
.
~
o
~ O 0 0 0 0 0 N
~
~ ~ ~
.
~
c
~
~
..,,
°°
g
So
o
I"1
z
t..,I-. ,.0
>
.,..a
Heat transport in geological formations 231 241 251 261 271 284 291 301 311 321 331 34t 35f 36i 371 38J 391 401
188 194 200 206 212 218 224 231 238 245 252 259 266 273 280 287 293 300
equil?
188 194 200 206 212 218 225 231 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 218 225 231 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 218 225 232 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 219 225 232 238 245 252 259 266 273 280 287 294 300
188 194 200 206 212 219 225 232 239 245 252 259 266 273 280 287 294 300
188 194 200 206 213 219 225 232 239 246 253 260 266 273 280 287 294 300
N
N
N
N
N
N
213 219 226 232 239 246 253 260 267 274 280 287 294 300
188 194 201 207 213 219 226 232 239 246 253 260 267 274 280 287 294 300
188 194 201 207 213 219 226 233 239 246 253 260 267 274 280 287 294 300
188 194 201 207 213 220 226 233 239 246 253 260 267 274 280 287 294 300
188 195 201 207 213 220 226 233 239 246 253 260 267 274 280 287 294 300
N
N
N
N
N
N
48 44
49 45
50 46
51 47
52 48
53 49
TEMPERATURE PROFILE EVOLUTION (TPE) i0 i0 i0 i0 I0 I0 i0 I0 18 18 18 18 18 18 18 18 26 26 27 27 27 27 27 27 34 34 34 35 35 35 35 35 42 42 42 42 42 43 43 43 50 50 50 50 50 50 50 51 58 58 58 58 58 58 58 58 65 66 66 66 66 66 66 66 73 73 73 74 74 74 74 74 81 81 81 81 82 82 82 82 89 89 89 89 90 90 90 90 97 97 97 97 98 98 98 98 105 105 105 106 106 106 106 106 113 113 114 114 114 114 114 114 122 122 122 122 122 122 122 122 130 130 130 130 130 130 130 130 138 138 138 138 138 138 138 138 146 146 146 146 146 146 146 146 153 153 153 153 153 153 153 153 161 161 161 161 161 161 161 161 168 168 168 168 168 168 168 168 175 175 175 175 175 175 175 175 182 182 182 182 182 182 182 182 188 188 188 188 188 188 188 188 195 195 195 195 195 195 195 195 201 201 201 201 201 201 201 201 207 207 207 208 208 208 208 208 214 214 214 214 214 214 214 214 220 220 220 220 221 221 221 221 227 227 227 227 227 227 227 227 233 233 233 234 234 234 234 234 240 240 240 240 240 240 241 241 247 247 247 247 247 247 247 247 253 254 254 254 254 254 254 254 260 260 260 261 261 261 261 261 267 267 267 267 267 267 268 268 274 274 274 274 274 274 274 274 281 281 281 281 281 281 281 281 287 287 287 287 287 287 287 287 294 294 294 294 294 294 294 294 300 300 300 300 300 300 300 300
i0 1B 27 35 43 51 58 66 74 82 90 98 106 114 122 130 138 146 153 161 168 175 182 189 195 202 208 214 221 227 234 241 247 254 261 268 274 281 287 294 300
i0 18 27 35 43 51 59 66 74 82 90 98 106 114 122 130 138 146 153 161 168 175 182 189 195 202 208 215 221 228 234 241 248 254 261 268 274 281 287 294 300
N
N
Enter n-value o f l e f t m o s t 4O
column
n: age:
40 37
41 38
42 39
44 41
il OI I) 21 31 4t 5f 61 71 81 9[ 10r III 121 13i 141 151 161 17t 181 19f 201 211 221 231 241 25[ 261 27[ 281 29] 301 311 32i 331 341 351 361 37[ 381 39i 40[
i0 18 26 34 42 50 57 65 73 81 89 97 105 113 121 129 138 146 153 161 168 175 182 188 195 201 207 213 220 226 233 239 246 253 260 267 274 280 287 294 300
I0 18 26 34 42 50 57 65 73 81 89 97 105 113 121 130 138 146 153 161 168 175 182 188 195 201 207 213 220 226 233 240 246 253 260 267 274 281 287 294 300
I0 18 26 34 42 50 58 65 73 81 89 97 105 113 121 130 138 146 153 161 168 175 182 188 195 201 207 214 220 226 233 240 247 253 260 267 274 281 287 294 300
N
N
equil?
43 40
N
N
188 194 200
188 194 200
206
207
213 219 225 232 239 246 253 260 267 273 280 287 294 300 N
689 188 5750 1 9 5 ' 6000 2 0 1 6250 207 6500 2131 6750 220~ 7000 226 7250 233 7500 240 7750 246 8000 253 8250 260 8500 267 8750 274 9000 281 9250 287 9500 294 9750 300 I0000
(il to exit):
45 42
N
46 43
N
47 44
N
N
N
N
kyrs
I z i01 0 191 250 271 500 351 750 43J 1000 511 1250 591 1500 671 1750 741 2000 82( 2250 901 2500 98# 2750 106 3000 114 3250 122 3500 130 3750 138 4000 146 4250 153 4500 161 4750 168 5000 175 5250 182 5500 189 575O 195 6000 202 6250 208 6500 215 6750 221 7000 228 7250 234 7500 241 7750 248 8000 254 8250 261 8500 268 8750 274 9000 281 9250 288 9500 294 9750 300 i0000 N