Mathematical simulation of geotemperature and heat flow patterns

JOURNAL OF GEODYNAMICS4, 45--61 (1985)

MATHEMATICAL SIMULATION HEAT FLOW PATTERNS

45

OF G E O T E M P E R A T U R E

AND

X I O N G L1ANG-PING, Z H A N G J U - M I N G A N D S U N H U I - W E N

Institute ~?/' Geology. Academ/a Sinica, Be~jing, China (Accepted July 4, 1985 )

ABSTRACT Xiong kiung-ping, Zhang Ju-ming and Sun Hui-wcn, 1985. Mathematical simulation of geotcmperaturc and heat flow patterns. In: k. Rybach editor}, Heat Flow and Geothermal Processes. Journal o~ Geo~lvnamics, 4:45 61. Geotemperature and heat flow patterns m a large-scale Meso-Cenozoic basin such as the North China Basin are strongly affected by the relief of the basement, and controlled by the contrast of thermal conductivity between basement rock and sedimentary cover. Usually, heat flow observed at the surface of a basement uplift is greater than that of a basement depression. Calculation revealed, that the ratio of the former and the latter is determined by the uplifted height (H) of the bed-rock roof of the basement and the thickness (h) of the sedimentary cover. The relief of the basement also disturbs the geotemperature and, hence, the heat flow patterns at shallow depth. Consequently, the more or less "'uniform" one dimensional heat flow from the deep interior of the Earth becomes two dimensional at shallow depth with great lateral and vertical variations. The extent of the disturbed zone is also controlled by the contrast of the thermal conductivity between basement rock and sedimentary cover as well as the uplifted heigh (Hi of the bed-rock roof of the basement. Numerical computation demonstrated that the disturbed depth (Ze) is usually about 3 6 times of the uplifted height (tt) of a basement uplift.

INTRODUCTION

Just as Jaeger pointed out in early 1960's, that the heat flow measured on land in boreholes may be influenced by transient disturbances, mass movement and material inhomogeneity (Jaeger, 1965). In North China Basin, this problem is especially prominent because the relief of basement beneath thick Cenozoic sediments is unflat and rugged, as well as faults and fractures are extensive, causing groun water circulation and inhomogeneity of physical properties of the medium. Observations indicate that geotemperature and heat flow patterns in North China Basin are strongly affected by the relief of basement rock and controlled by the contrast of the thermal conductivity between basement rock and sedimentary cover (Wang Ji-yang et al., 1981; Geothermal Research Division, 1980). 0264-3707/'85/$3.00

~" 1985 Geophysical Press Ltd.

46

XIONG, ZHANG, AND SUN

In this paper, mathematical simulation of geotemperature and heat flow patterns by finite element technique has been conducted and some results are discussed. METHOD USED IN MATHEMATICAL SIMULATION

Heat transfer in the uppermost part of the Earth's crust is almost entirely realized by two forms: conduction and convection. In cases where the heat flowing towards the Earth's surface has been transferred by conducdtion, the general equation governing heat conducdtion in two dimensional form is:

0 /KOT\ Where

0 /OT\

~T

7 " temperature t--time K-- thermal conductivity Q--heat generation ~--heat capacity

Since early 1970's, finite element technique has been applied in geothermal studies by several authors (Geertsma, 1971; Lee and Henyey, 1974; Lee, 1975; Jones and Oxburgh, 1979). Recently, we have employed this technique to investigate various problems of geothermics and some results were obtained (Zhang Ju-ming et al., 1982, Xiong Liang-ping et al., 1983, 1984). As the finite element technique is well known and widely accepted by geothermal community nowadays, we only deal with some points in our computation.

1. Determination of anomalous areas Since the physical properties of the medium considered may vary, we regarded the whole region as homogeneous and gave it a set of physical parameters which we call the regional physical parameters. The areas where the physical properties are different from the regional ones were delimited by some closed curved lines. These areas are called anomalous areas, each of which has its corresponding anomalous physical parameters. For each quadrilateral, if its centroid falls into an anomalous area, the two triangular element of this quadrilateral will get the physical parameters of corresponding anomalous area. If the centroid is out of all these areas, the two elements will get the regional physical parameters. We can use the method described below to judge the position of the cen-

MATHEMATICAL S I M U L A T I O N O F G E O T E M P E R A T U R E

I

Fig. 1.

i

47

[

Judgemenl of anomalous areas. Centroid 1 and 3 in the area, 2 out of the area.

troid to see if it is in or out of an anomalous area. It is imaginable that there would be a vertical straight line through this centroid and the line would intersect the closed curve surrounding the anomalous area. Then we count the intersections over the centroid. If this is an odd number, it shows that the centroid falls into the anomalous area, otherwise the centroid is out of this area. In order to simplify computation, these closed curves will actually be replaced by broken lines (polygons) approximating these curved lines (Fig. 1).

2. Constraining the boundary and initial conditions The N nodal points on the upper row are the upper boundary nodal points which are approaching the Earth's surface. We give each of them a given temperature or vertical heat flux for each point.

3. Plotting the .final results There are often great amount of data set in computation such as the coordinate values of each nodal point, the temperature value on each point, geothermal gradient value and heat flux value within each element. These data may all be plotted with X-Y plotter. MATHEMATICAL M O D E L AND RESULTS O F C O M P U T A T I O N

To understand and evaluate the effect of relief of basement rock on geotemperature and heat flow patterns in North China Basin, a mathematical model has been established. This model is composed of two layers (basement rock and sedimentary cover) with different conductivity,

48

XIONG, ZHANG, AND SUN

To h

K1

K2

q* Fig. 2.

Mathematical model simulating uplift and depression of basement rock.

thickness and configuration, simulating uplift (right in Fig. 2) and depression (left in Fig. 2) of the basement rock. The upper boundary has been regarded as an isothermal surface with To = 15°C in all computations, whereas the heat flux at lower boundary (q*) is allowed to vary from 33.5 to 83.7 m W / m 2 in different computations. The basement relief is described by the thickness (h) of sedimentary cover, the uplifted heigh (H) of basement rock and dip angle (0) of boundary between these two strata. The results obtained are as follows:

1. Effect of the thickness of sedimentary cover Field observations in North China Basin exhibit that the geothermal gradient of Cenozoic sedimentary cover strata is closely correlated to the thickness of these strata (Fig. 3). Based on numerous measurements of rock samples, the mean thermal conductivity of Cenozoic sedimentary cover (K~) in North China Basin is typically between 1.25 and 2.09 W/re.K, and of Pre-Cenozoic basement rock (K2), typically between 2.50 and 4.20 W/m.K. In our computation, a ratio Kz/K~ = 2.86 has been chosen simulating the actual situation of Rongcheng-Niotuozhen uplit where KI = 1.46 and K2--4.18 W/m.K were determined. The dip angle (O) of boundary between sedimentary cover and basement rock has been defined as tan 0 = 1.5 and the uplifted height (H) of basement rock was given as H = constant. Assuming the depth for model computation H * = 6 H , h/H=O.025, 0.125, 0.25, 0.35, 0.55 and 0.65 were selected in various computations respectively. In addition, the model has been considered to be an axis symmetrical one.

MATHEMATICAL SIMULATION OF G E O T E M P E R A T U R E

Geothermal 1

2

3

gradient 4

5

49

of C z s e d i m e n t s 6

7

8

9

°C/

100m

0

7 3 m

•

e-

2000 ~Z ©

¢~ -~

3000

¢-

z

E

~: ~

4000 o

5000

m

: ,

Fig. 3. Geothermal gradient of Cenozoic sedimentary cover versus depth of bedrock roof of lhc basement in North China Basin. I. A small subuplift in Beijing Depression;2. Northern part of Niotuozhcn Uplift:3. Southern part of Niotuozhen Uplift;4. Ningjin-Xinhc Uplift;5. Dachcng Uplift

Results thus obtained revealed that the calculated average surface heat flow on top of a basement uplift (q,) is always greater than that on top of a depression (q,). Furthermore, the values of q,, q,/q., and q,,/q* are decreasing with the increase of the relative thickness of sedimentary cover (h/H). whereas q., increases with incresing of h/H. It has been noticed that q, and q, tend to reach the value of 62.8 mW/m 2, which equals the heat flux at lower boundary (q*) in this computation while h/H= 0.83. This implies that

50

XIONG, ZHANG, AND SUN

mW/m2

1

90

if_ 70

~

"~ "" "

qu s

50

qu l q s qu/q*

40

K2/K1

30 0

'

O.l

q* - - - 6 2 . 8

--- 2.86 '

0.2

'

0.3

'

0.4

'

0.5

'

0.6

mW / m2 '

0.7

'

0.8

~

0.9

h/H l

Fig. 4. Averagedsurface heat flow on top of basement uplift and depression versus relative thickness of sedimentary cover.

the effect of basement relief upon surface heat flow diminishes to the minimum, i.e. the effect of basement relief already could not been examined at the surface, while h = 0.83H. It must be stressed that the temperature observed at the same level in basement uplift zone is always lower than that in depression zone, but the heat flow in basement uplift zone seems to be still higher than that in depression zone, while there is no sediments on top of a basement uplift ( h = 0 ) . In this case, a reverse picture of isotherm appeared (For further details see Xiong Liang-ping and Gao Wei-an, 1982).

2. Effect of the ratio of thermal conductivity of basement rock and sedimentarv cover

As mentioned above, the ratio of thermal conductivity of basement rock and sedimentary cover (K2/K 1) is an important parameter affecting geotemperature and heat flow patterns in North China Basin. Many investigators (Lachenbruch and Marshall, 1966; Hyndmann and Sass, 1966) have repor-

MATHEMATICAL SIMULATION OF GEOTEMPERATURE

51

qu/qs %

2.0

/

/

/

/

1.0

/

/

/

Fig. 5.

/

/

/

/

/

/

//

|

0

/

i

1.0

-

2.0

/

3.0

K2

!

K1

Surface heat flow in uplift and depression zones versus ratio of thermal conductivity for the

model with h/H = 1/3 and q* = 62.8 m W / m 2.

ted that the ratio of heat flow at both sides of the boundary between two rock strata with different thermal conductivities should equal the ratio of thermal conductivity of these rock strata. In our model, it means that the ratio of heat flow observed in basement uplift and depression zones should equal the ratio of thermal conductivity of rock strata at both sides of the boundary between basement and sediments. According to the restllts of our model calculation, this is only the particular situation for the case while there is no sedimentary cover on top of basement, i.e. basement rock exposed at the Earth's surface ( h = 0 ) . The general case is that q,,/q, increases with increasing of K2/K~, and even the ratio of the maximum heat flow appeared in basement uplift and depression zones is always less than the ratio of K2/KI, while exists sedimentary cover on top of a basement uplift (h -¢0) (Fig. 5 and Table 1).

3. Effect of the dip angle q[ boundary The effect of the dip angle (cf. Fig. 2) is also related to the thickness of sedimentary cover and the ratio of thermal conductivity. In our calculation, TABLE 1

Heat flow variation with the ratio of thermal conductivity K~.K I

1.14

1.71

2.86

3.43

4.00

5.00

q,, max/'q.....

1.04

1.18

1.33

1.38

1.66

1.82

52

XIONG, ZHANG, AND SUN

h/H= 1/3, K2/K 1=2.86 and q * = 62.8 mW/m 2 were chosen. Under these conditions, temperature in the sedimentary cover is increasing with diminishing of dip angle0, and the maximum temperature increase (about 15°C) occurs at the boundary (Fig. 6). The effect of bsement relief on geotemperature and heat flow patterns is used to be explained by heat flow "refraction" and "boundary effect" Temperature

"C

profile

I

I

90, ~G

- - - - - - ' -. .~. . .

. ._. .

3

"'-'----" 1 80

Temperature

profile

70

~ ,i.

II

II

5 _ -2-'_ . . . . . . .

3_

,7

5o---~4~7.'-.

.

2

4

6

.

.

8

t ~ b l

I

[

II

.

.

.

10

b2,,

.

.

.

12

t

.

14

16

1

•

~

x ~ K1

K2

Fig. 6. Temperature profile with different dip angle of boundary at different depth. (Figures 1, 2, 3, 4, 5 represent models with different dip angles of boundary; b--top width of the uplifted part of a basement uplift. In our calculation, assume bl = b2 = b3 = b4 = b s }

MATHEMATICAL SIMULATION OF GEOTEMPERATURE

53

(Lachenbruch and Marshall, 1966; Hyndmann and Sass, 1966). In our opinion, the problem is far more complicated than simple "refraction" or "boundary effect", because the heat always flows towards the place with lowest total thermal resistance, and hence, the geotemperature and heat flow patterns are constrained by the whole thermal resistance field rather than the partial resistance field near the boundary. From the conceptual model in Fig. 7 it is clear that in the basement, the heat flow from below is directed towards the uplift zone, where the total thermal resistance is much less than in the depression zone with thick sediments of low thermal conductivity. The isotherm configuration in basement (lower part in Fig. 7), therefore, is concave, i.e. at a certain level the temperature in depression (T~) is higher than that in uplift (T,,), and vice versa for the heat flow, The difference between Tu and T, as well as qu and q,. reaches the maximum at the interface between basement and sediments where qu>>q.,.. In the sediments above a basement uplift (upper part in Fig. 7) and at a certain

JTI

Ts
he

-~

h

]_ Be

~

qs <

qu

K2 Ts > T u Ze

q* Fig. 7. zones.

Conceptual modelshowingthe redistributionof heat flow in basementuplift and depression

54

XIONG, ZHANG, AND SUN

level, 7",> T~., because there the q, still greater than q,~. This lateral difference of temperature in sediments above the basement uplift forced heat flow from uplift back to depression zone again, until equilibrium was established. We have defined the whole process, that is, the heat flows towards uplift zone in basement rock and vice versa in sedimentary cover, as "redistribution of heat flow". Because of redistribution of heat flow, the more or less uniform one dimensional heat flow from the deep interior of the Earth became two dimensional flow at shallow depth. The two dimensional heat flow zone resulted from the effect of basement relief in the uppermost part of the Earth's crust has been named as "disturbed heat flow zone" (Xiong Liang-ping and Zhang Ju-ming, 1984). To specify the so called "disturbed heat flow zone" is of great importance in estimating the temperature in deep Earth's crust, usually derived from extrapolation of near-surface heat flow measurements, which are always biased by the relief of basement rock and "redistributed" at shallow depth. Within disturbed heat flow zone, any horizontal plane would not be treated an isotherm or uniform heat flow boundary. Otherwise, some difficulties and even mistakes may arise. The extent of disturbed heat flow zone could be described by the disturbed depth (Z,,), width (L,~) and height (h,,). Another words, the zone within depth (h,, + H + Z,,) and width (L,, + d) has been regarded as disturbed heat flow zone. In the area of (h-h,,), while h > h,,, the temperature and heat

K2/K

1

6.0 ~.0'

4.0 3.0 ¸ ~°0" i . O.

0

110

210

310

,

4.0

,

5.0

,

6.0

/

Ze! H

Fig. 8. Relationship between Z,., H and Kz/K~ in a basement uplift. (Explanation see text)

MATHEMATICAL SIMULATION OF GEOTEMPERATURE

55

flow are already beyond the effect of basement relief (Fig. 7). Meanwhile, the heat flow conducted from below the depth of ( h e + H + Z , , ) might be thought to be unaffected by the relief of basement, and hence, any horizontal plane below that depth could then be treated as uniform heat flow boundary or isothermal surface in deep temperature calculation. Our further results show that the disturbed depth (Z,,) is mainly controlled by the uplifted height (H) of basement and the ratio of thermal conductivity of basement rock and sedimentary cover (K2/K1). In most cases, Z,, is of about 3 to 6 times of H, for instance, Z,, equals 4.6 H in case that K2/KI = 2.86 (Fig. 8). The disturbed height (he) is also correlated with the ratio of K2/K~. It is clear from Fig. 4 that h,, equals about 0.83 H, while K2/K~ -- 2.86. Within depth (h - h,,), while h > h,,, the temperature and heat flow at a certain level above basement uplift and/or depression zones are nearly the same, the difference of which is less than 5% according to our calculation. In practice, this difference is negligible, and therefore, the area of (h-h,,) near the Earth's surface could be thought of an undisturbed heat flow zone too. The disturbed width (L,,) is of 1.5 to 2 times of H resulted from our previous work (Xiong Liang-ping and Gao Wei-an, 1982).

4. Effect q[' topography This problem may also be treated as redistribution of heat flow at shallow depth and at the Earth's surface, because the total thermal resistance at positive topograghy is always larger than that at negative one, and hence the heat flow would be concentrated in zones of negative topograghy. Our calculation for terrain correction of heat flow measured in Yang Zhuoyong Lake, Tibet suggested that the effect depth of topography on geotemperature and heat flow patterns could reach a depth of 5 Km, that is of about 6 times of the amplitude of topographic variation. The amount of terrain correction for heat flow measured in this lake appears to be within 10% (Fig. 9).

5, Effect of heat generation Our calculation reveals that the radiogenic heat, produced in the uppermost part of the Earth's crust, has also been flowing with the heat from the depth towards uplift zone. Essentially, this phenomenon may be thought of as the integrated part of the whole process of heat transfer in the Earth's crust and upper mantle. In Liaohe Basin of North China, the radiogenic heat generated within a depth of 10 Km and concentrated in the uplift zone is of about 20.93 mW/m 2, whereas in the depression zone, 15.07 mW/m 2.

56

XIONG, ZHANG, AND SUN /

0.998 H *

'('/k !

•

mW

//

m2

'

m

-60

Ii I

P ii

c

4o

'07 :C" ~-.-q"'U'~:"~22~.£"T--4x 0,992 H *

!•

f s°

. l l' j.

I

~-4o

'1

^

_

80"v'-' " ' N . / ' ~ ' - , ' ~ / "r'\ fJ "'f °0t . . . . - . . . . . . . . i ~ - ~ - ........

{;

~

k

120 ! F 4{I

0. 988 H *

8o~.j..f.,,~-Ji._, r.~.fG I-f2° 4o4- . . . . . . . . ' . . . . . . /""'"-" . . . . . . . ~ . . , qv | , , x . . . r---- ~ .X v---~"- ' ~'" 0. 905 H * 8{}~'--'[ ' " - ~ ' - -

40 T- . . . . .

"-/ ' " "~" " ~ ' . . . .

i , ~. . . . . . . __ • . -.,-. . . . . . . . .. .. . .

G

. . .

~30

O. 567 H * ............

G

........................... ,if}

I

,

1

,

E 3

,

,

,

5

,

qv ,

7

,

9

Temerature

1"0"~

~2{]

qq× v

'

,

11

'

~

13

''

km

profile

~

E2.0 ~

50 60

3.0

H*

70 80 90

4.0 100 °C

Fig. 9.

Effect of t o p o g r a p h y on geotemperature and heat flow patterns in Yang Z h u o y o n g Lake, Tibet.

MATHEMATICAL S I M U L A T I O N O F G E O T E M P E R A T U R E

57

CASE STUDIES

For better understanding the effect of basement relief upon geotemperature and heat flow patterns, a case study for Ji Yang Depression in North China has been carried out. The Ji Yang Depression is believed to be the deepest (12 Km) in North China Basin, and consists of Mesozoic and Cenozoic sediments and Pre-Mz basement rocks with different thermal conductivities (Table 2). TABLE 2 Thermal conductivity of rock strata in Ji Yang Depression, North China Strata

Conductivity ( W/m. K )

Q+ N

1.55

E

1.67

E,, + E k

1.88

Mz

2.09

P22

Pzl

AnZ

(Upper Paleozoic)

(Lower Paleozoic )

(PreSinian )

2.09

2.51

2.93

A simplified model for Ji Yang Depression was established in terms of geological and geophysical evidences so far obtained in the region studied. Depth of the model has been taken as 13 Km, and width of the model, 75 Km respectively. Having considered to be the upper boundary, the Earth's surface is regarded as isothermal, and the mean annual temperature (14°C) has been used in computation. The lower boundary is characterized by a variable heat flux at 13 Km, deduced from calculation for a simplified model with a depth of 30 Km, where the heat from the deep interior of the Earth must be considered to be uniform and constant (62.8 mW/m2). All the data computed are automatically plotted on a X - Y plotter, and a geothermal profile with heat flow, isotherm, geothermal gradient at different depth (100, 500, 1000, 2000, 4500 and 13,000m respectively) is presented in Fig. 10. Results obtained are as follows: 1. The effect of basement relief upon geotemperature and heat flow patterns is prominent within a depth of 13 Km. It can be seen in Fig. 10 that isotherm configuration is much more influenced at the lower part of the Depression, caused by strong undulation of Pre-Sinian (AnZ) basement rock with high thermal conductivity (Fig. 10 and Table 2). 2. The maximum depth of buried Pre-Sinian basement rock in Ji Yang Depression reaches 12 Km, and hence the effect depth of basement relief seems to be more than 40 Km. Therefore, heat flow at lower boundary (13 Km in deep) of this model is not uniform and constant, it varies from 50.24 to 66.99 mW/m 2 according to our calculation.

E

1

i

l(I ]"

{a

'('

-

1

n

'

lIj

12q~l

!

aUt . . . . .

[61 -

11 '

10t

807 -101

m ~,~/m

it)

'

/

15

•

vO

~

J

25

30

,

,

35 I0 Temperature

,

r

r

15 50 profile

1301)0 m

2000 m

100(I m

5()[)m

] OU m

:~:)

'

'

6'0

ti5

,

--

- - - -

80

70

7'5

,

r

IO

200

60

20

40

~C

k m

G f 20 30 q)

t; I~0 qy

qGy

4U I2 0

O ~50 q y f 30

G 15030 qy

c; [ 5(I q x [30

C'/km

z

7

z

N 'w

5Z

X

MATHEMATICAL SIMULATION O F G E O T E M P E R A T U R E

59

. As shown in Table 3 and Fig. 10 that the lateral variation of heat flow across the profile is much more clear at rather shallow depth than at surface or at depth, the maximum variation (qmax--qm~n= 36.67 mW/m 2) has been observed at 1000 m depth, where the qma× is of 1.67 qmi." TABLE 3 Lateral variation of heat flow with depth in Ji Yang Depression, North China Depth (m)

100

200

1000

1500

2000

2750

4500

8000

13000

q ..... (mW/m ~)

75.45

88.04

91.36

80.00

76.95

78.38

69.46

67.70

67.00

qmin (mW,,'m e

54.85

54.81

54.68

54.51

53.59

51.76

46.43

50.83

53.46

qinax qmm (mW,,'m 2

20.60

33.23

36.68

25.49

23.36

26.62

23.03

16.87

13.54

. Heat flow varies with depth, and the amplitude of variation in basement uplift zone appears to be greater than that in depression zone (Table 4). TABLE 4 Heat flow variation with depth in basement uplift zone (representative at point of 8.5 Km from left side in Fig. 10) and depression zone (37.5 Km from left side in Fig. 10) (q~--- heat flow in uplift zone; q~-- heat flow in depression zone) Depth (m) qu

100

500

1000

1500

2000

2750

4500

8000

13000

71.05

62.84

91.36

73.35

73.10

71.22

68.41

66.74

65.86

55.39

55.39

55.39

55.10

54.74

54.22

53.34

51.83

52.46

mW,,m 2) q, (mW,,'m 2 )

Fig. 10. Heat flow, geothermal gradient and isotherm calculated for the simplified model of Ji Yang Depression in North China Basin. (For comparison, the observed values of geothermal gradient are plotted on top)

XIONG, ZHANG, AND SUN

60

5. With only a few exceptions, the calculated values are in good accordance with the observations in the field within 1 Km depth (Fig. 10). CONCLUSION

I. Geotemperature and heat flow patterns in sedimentary basins are strongly affected by the basement relief and by the contrast of thermal conductivity between basement rock and sedimentary cover; 2. Sedimentary cover with lower thermal conducitvity on top of a basement uplift plays an important role in forming the total thermal resistance of rock strata, and hence the resulted geotemperature and heat flow patterns; 3. Redistribution of heat flow in the uppermost part of the Earth's crust is a far more complicated phenomenon than simple refraction of heat flow and/or boundary effect, and is controlled by the total thermal resistance of rock strata; 4. To specify the so called "disturbed heat flow zone" is of great importance in estimating the temperature at depth; 5. Effect of topography upon geotemperature and heat flow patterns could also be treated as problem of redistribution of heat flow at shallow depth and at the Earth's surface; 6. As an integrated part of the whole process of heat transfer in the Earth's crust and upper mantle, radiogenic heat produced in the uppermost part of the Earth's crust has also been flowing with the heat arising from the deep interior of the Earth towards uplift zones; 7. As convincingly stated in this paper, finite element technique may be thought of to be a powerful approach in regional heat flow studies.

REFERENCES Chen Mo-xiang, Huang Ge-shan and Zhang Wen-ren, 1982. Distributive characteristics of the geothermal resources in the northern area of the North China Plain, in: Research on Geology ( 1 ), Institute of Geology, Academia Sinica, Cultural Relics Publishing House, pp. 72 79. (in Chinese) Geertsma, J., 1971. Finite element analysis of shallow temperature anomalies, Geophys. Prospect. 19: 662-681. Geothermal Research Division, Institute of Geology, Academia Sinica, 1980. Geothermal characteristics of North China Plain and its surroundings, in Formation and Development of North Fault Block Region, Science Press, Beijing, pp. 368-378. (in Chinese) Hyndmann, R. D. and J. H. Sass, 1966. Geothermal measurement at Mount Isa~ Queensland, J. Geophys. Res., 71(2): 587 601.

MATHEMATICAL SIMULATION OF GEOTEMPERATURE

61

Jaeger, J. C., 1965. Application of the theory of heat conduction to geothermal measurements, in: Terrestrial Heat Flow, Am. Geophys. Union, Monograph 8, Washington, D.C., pp. 7' 23. Jones, F. W. and E. R. Oxburgh, 1979. Two-dimensional thermal conductivity anomalie~ and vertical heat flow variations, in: Terrestrial Heat Flow in Europe, Inter-Union Commission of Geodynamics Scientific Report No. 58, Springer-Verlag, Berlin-New York, pp. 98-106. Lachenbruch, A. H. and B. V. Marshall, 1966. Heat flow through the Arctic Ocean floor: Canada BasinAlpha Rise boundary, J. Geophys. Res., 71(4): 1223 1248. Lee, T. C. and T. L. Henyey, 1974. Heat flow refraction across dissimilar media, Geophys. J. R. Astron. Soc., 39: 319-333. LEE, Y. C., 1975. Focusing and defocusing of heat flow by a buried sphere, Geophys. J. R. Astron. Sot., 43: 635-641. Lee, T. C., Rudman, A. J. and A. Sjoreen, 1980. Application of finite element analysis to terrestrial heat flow, Occasional Paper 29, Department of Natural Resources, Geological Survey of State Indiana, U.S.A. Wang Ji-yang, Chen Mo-xiang, Wang Ji-an, Deng Xiao, Wang Jun, Shen Hsien-chieh, Xiong Liangping, Yang Shu-zhen, Fan Zhi-cheng, Liu Xiu-wen, Huang Ge-shan, Zhang Wen-ten, Shao Hai-hui and Zhang Rong-yan, 1981. Geothermal studies in China, J. of Volcanology and Geothermal Research, 9: 57-76. Xiong Liang-ping and Gao Wei-an, 1982. Characteristics of geotherm in uplift and depression, Acta Geophysica Sinica, 25(5): 448~456. {in Chinese) Xiong Liang-ping, Zhang Ju-ming and Sun Hui-wen, 1983. Study of the effect of regional slructure form on the thermal field by finite element analysis, Acta Geophysica Sinica, supplement 26:751 763 (in Chinese ). Xiong Liang-ping and Zhang Ju-ming, 1984. Mathematical simulation of heat flow refraction and redistribution, Scientia Geologica Sinica, 4:445~454 (in Chinese). Zhang Ju-ming, Sun Hui-wen and Xiong Liang-ping, 1982. Mathematical simulation of regional geothermal field and case analysis, Scientia Geologica Sinica, 3:315 321. (in Chinese)