Some thermal and mechanical consequences of rapid uplift: an example from the Southern Alps, New Zealand

Earth and Planetary Science Letters, 86 (1987) 307-319 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

307

[5]

Some thermal and mechanical consequences of rapid uplift: an example from the Southern Alps, New Zealand P.O. Koons

Geology Department, University of Otago, Box 56, Dunedin (New Zealand) Received December 17, 1986; revised version received August 28, 1987 The thermal evolution of continental crust during active collision is modeled through numerical solutions of the two-dimensional heat conduction equation for a rapidly moving medium. The boundary conditions used in the modeling are derived from geological and geophysical observations from the active collision zone in the South Island of New Zealand. The problem domain over which the solutions are obtained consists of a 40 km horizontal by 25 km vertical spatial plane with a vertical discontinuity at 10 km from the western boundary. To the east of this discontinuity, vertical uplift rates of up to 10 m m / a occur over a timespan of up to 4 Ma. Temperature distributions are calculated at 10 ka intervals over the 4 Ma duration. A two-dimensional high-temperature region is established upon initiation of uplift of the eastern block due to the advective component carrying heat upwards more rapidly than it can be dissipated laterally from the problem domain. Temperatures within the upper 5 km are greater than 400 o C after 2.25 Ma with geothermal gradients of up to 200 o C / k m attained within the upper 3 km. At times greater than 2.5 Ma, the vertical temperature distribution changes little while the anomalously high temperatures spread laterally into the stationary crust. Using rheological equations to describe the brittle behaviour of a water-saturated upper crust and the ductile behaviour of a quartz-dominated lower crust, together with the thermal distribution of the conduction models, the mechanical evolution of a collision zone is investigated. In addition to high crustal temperatures and associated high heat flow, rapid uplift produces a weakening of the crust by raising of the depth of transition from brittle to ductile behaviour. Within the zone of most rapid uplift, the brittle-ductile transition rises from - 13 km to less than 5 km after 1.5 Ma of uplift. Further uplift reduces the brittle layer t o - 3 km thickness and causes lateral spreading of the low-strength zone. The reductions in crustal strength caused by the thermal weakening produce a high-strain zone within the region of maximum uplift which is incapable of sustaining large differential stresses. This causes horizontal and vertical stress transfer and results in shallow seismicity increases in the adjacent crust as well as in intermediate depth seismicity within the high-strength upper mantle. Because the thermal and mechanical anomalies discussed are a function of rapid uplift, all regions of active continental collision may be expected to exhibit similar behaviour. Some mechanical and thermal characteristics of the Himalayan collision zone are briefly examined in light of the numerical modeling.

1. Introduction Thermal evolution of continental lithosphere d u r i n g a n d a f t e r c o n t i n e n t a l c o l l i s i o n has b e e n d i s c u s s e d in r e c e n t y e a r s b y n u m e r o u s a u t h o r s concerned with diverse aspects of continental coll i s i o n [1]. T h e p h e n o m e n a i n v e s t i g a t e d i n c l u d e the g e n e r a t i o n o f m e t a m o r p h i c facies d u r i n g c o l l i s i o n (e.g. [2-4]), l i t h o s p h e r i c v i s c o s i t y effects [5] a n d p a r t i a l m e l t i n g o f t h e c r u s t d u r i n g a n d a f t e r collis i o n (e.g. [6,7]). T h e m o d e l s p r o d u c e d in these studies are primarily concerned with the characterisation of thermal profiles during and after c o n t i n e n t a l c o l l i s i o n b y c a l c u l a t i o n o f t h e r m a l re0012-821X/87/$03.50

© 1987 Elsevier Science Publishers B.V.

l a x a t i o n f o l l o w i n g a discreet, i n s t a n t a n e o u s episode of thrusting. Although these models require r a p i d t e c t o n i c b e h a v i o u r to p r o d u c e the initial l i t h o s p h e r e p e r t u r b a t i o n , m o s t o f these m o d e l s use t h e a s s u m p t i o n t h a t u p l i f t rates are e x p l i c i t l y o r i m p l i c i t l y c o n t r o l l e d b y e r o s i o n rates t h r o u g h a n isostatic response mechanism. Because the authors r e s t r i c t the u p l i f t r a t e s to less t h a n 1 m m / a , the t h e r m a l i m p a c t o f u p l i f t in t h e s e slow u p l i f t m o d els is p r i m a r i l y to r e m o v e the h e a t p r o d u c i n g e l e m e n t s t h r o u g h e r o s i o n a n d t h e r e f o r e serves as a n e t h e a t sink. If the v i s c o e l a s t i c r e p o n s e o f the a e s t h e n o s p h e r e a n d l i t h o s p h e r e to v e r t i c a l l o a d i n g is c o n s i d e r e d to c o n t r o l u p l i f t r a t e s f o l l o w i n g c o n -

308 tinental collision, then conductive time constants in thermal models of continental collision are on the order of 250 Ma [8]. England and Richardson [3] have investigated the thermal effects of erosion and uplift on continental heat flow over time scales of 50-200 Ma. Subsequently, England and Thompson [4] and Thompson and England [9] used these time constants to determine the pressure-temperature-time paths of continental material undergoing metamorphism during burial and subsequent exhumation. As shown b y England and Thompson [4, Appendix C] uplift rates of ~< 1 m m / a result in low Peclet numbers and domination of the thermal budget by the two thermal components of (1) heat production due to radioactive decay, and (2) the depletion of the heat-producing elements due to erosion. Therefore, the geotherms in the burial-exhumation thermal models remain, in the general case, below the equilibrium geotherm (cf. [4, fig. 2]). The surface heat flow generated in the slow uplift models is consistently low and consequently regions of active collision are often assumed by authors interested in other phenomena to have depressed thermal gradients at depth [10,11]. An increasing amount of evidence from different collision zones and compiled from interpretations of fission track data (Himalaya [12]), uplifted terraces (New Zealand [13]; Alaska [14]) and isotope systematics (Scottish Dalradians [15]) demonstrates that uplift rates as high as 10 m m / a are common during orogeny and that these uplift rates may be sustained for millions of years [16]. An examination of regions of continental collision reveals the common association of thermal springs with regions of active collision, providing evidence of higher than average crustal heat flow [1]. The lateral dimensions of high uplift zones are apparently on the order of 100-300 km by 10-50 km and as such do not represent the entire deformation zone of a collision belt. In each of these examples the sustained high uplift rates are due to tectonic uplift and not primarily due to an isostatic, viscoelastic response function (cf. [8]). The mechanisms for the partitioning of strain into overthrnst and underplating components are discussed elsewhere. The effect on the thermal budget of the time scale shift from one of hundreds of millions of years to one of less than five million years is to

strongly emphasise the component of thermal advection due to uplift while suppressing the thermal effects of heat production through radioactive decay. This in turn largely eliminates the cooling impact of erosion of heat producing elements and results in a marked increase in the Peclet number. The thermal effects of rapid uplift were recognised by Allis et al. [17] who modeled rapid uplift along the Alpine Fault in one-dimension and demonstrated the high thermal anomaly associated with such uplift. This paper, using uplift rates and petrological relationships observed in the Southern Alps of New Zealand to provide boundary conditions for the numerical modeling, examines the thermal state of a region which has undergone sustained rapid uplift through finite difference solutions to the two-dimensional equation for heat conduction in a moving medium:

3T/Ot = K(a2T/O2x + ~2T/~Zz) + v(OT/Dz) in which T = temperature (given here as o C), K = thermal diffusivity, v = uplift velocity and DT/~z = the geothermal gradient. Due to the relatively short time scale of the processes discussed in this paper (~< 4.0 Ma), the thermal component due to radioactive decay is ignored [4]. The implications of the thermal state for the mechanical behaviour of the continental crust are then developed. Due to the uncertainties in the~ material properties of the rocks within any collision zone and local strain variation, this model is not intended to reproduce precisely the thermal and mechanical conditions of the Southern Alps collision zone. Instead I shall use the relatively accessible and tectonically simple Southern Alps as an example of general thermal and mechanical behaviour in active collision zones undergoing rapid uplift. Consequences of rapid uplift on fluid flow and hydrothermal mineralisation, petrologic and isotopic systems and geomorphologic development are the subject of ongoing research and will be discussed elsewhere. 2. The physical and numerical model; description and rationalisation The tectonic setting along the Australian-Pacific plate boundary has been extensively described (e.g. [18,19]). Although some disagreement exists among

309

AUSTRALIPLAT~ AN

,21 ,i,2 PACIF TE i 0

,

I

lOOkm

Fig. 1. Location diagram illustrating observed geological and geophysical characteristics of the central South Island which provide the constraints on the boundary conditions for the numerical models [18,19]. The horizontal dimension of the two-dimensional models is shown by the heavy bar across the Alpine Fault and spans the region from 10 km to the west of the Alpine Fault to 30 km east of the fault. Schist and gneiss belonging to the Haast Schist are exposed within the region of rapid uplift with the highest metamorphic grade exposed adjacent to the fault. Greywackes of the Torlesse Supergroup are exposed on the eastern edge of the problem domain (e.g. [20-23]. The Ostler Fault lies within the region covered by the Pukaki Microseismic Network [22,46].

those authors, there is widespead consensus that transpressive continental collision is currently occurring along the Alpine Fault. Plate tectonic estimates of - 46 m m / a along the plate boundary may be resolved i n t o - 22 m m / a normal to the fault a n d - 40 m m / a parallel to the fault [19] (Fig. 1). Movement normal to the fault has been occurring for at least the last 2.5 Ma and perhaps for the last 5 Ma. This component of compression has resulted in relatively high uplift rates with determinations ranging from - 20 m m / a [20] to 5 m m / a [13]. Bull and Cooper [13] present a thorough review of the uplift data and on the basis of their own study of elevated marine terrace sequences suggest current uplift rates of 7.8 m m / a near the region of greatest elevation and slightly lower rates to the north and south. Uplift rates are less well constrained at the

eastern edge of the problem domain, however, fault-monitoring geodetic arrays at the active Ostler Fault, 60 km east of the Alpine Fault yield recent uplift rates of ~< 1 m m / a [21]. On a longer time scale, Wellman [22] has estimated uplift rates o f - 1 m m / a from displaced lake terraces in the same region. Cooper [23], in his description of a kyanite-talc assemblage from the Haast schist, estimated that the total Kaikouran uplift has exposed a crustal section of greater than 25 km thickness from a prehnite-pumpelleyite upper level to a depth capable of developing the pressures responsible for the formation of the kyanite-talc assemblage. The initiation of uplift along the Alpine Fault is dated by the first appearance of greywacke cobbles derived from the Torlesse Supergroup on the Pacific Plate side of the Alpine fault in West Coast basin conglomerates at 2.5-3.0 Ma [24]. Rattenbury [25] has suggested that uplift began in the late Pleistocene rather than the early Pleistocene on the basis of Haast Schist overthrust onto Pleistocene gravels. The problem domain over which the numerical solution is applied is a region which extends over 40 km horizontally and 25 km vertically (see Fig. 2). A vertical discontinuity exists at 10 km from the western edge along which uplift of the eastern block occurs. The reference frame is of a Lagrangian type and therefore not intended to trace the horizontal displacement of the tectonic and thermal disturbance during the solutions. Because the vertical velocity component alone is examined, the position and orientation of the vertical discontinuity is free to move as a function of the horizontal velocity. This frame of reference causes the magnitude of the anomalies to be insensitive to fault dip, although the relative positon of the thermal anomalies will shift as a function of fault dip. Using the above observations on the timing and extent of continental collision, the following uplift scheme is used to produce the thermal perturbation in the numerical models. The uplifted region consists of a 5 km-wide region immediately east of the fault at x = 10 with a maximum uplift rate of 10 m m / a . Further to the east, the uplift rate decays linearly to 0 m m / a at 25 km from the fault (x = 35). The remaining regions undergo no uplift. A very large number of possible geometries and velocities may be used in the uplift models and I

310 have not attempted to reproduce them all. Instead, I shall discuss the thermal and mechanical results of this one uplift scheme together with a brief evaluation of the effects of varying uplift parameters. A geothermal gradient of 25 ° C / k m on both sides of the discontinuity has been assumed to provide the initial thermal structure of the problem domain. No-flow conditions are constructed on the western and eastern edges of the domain (x = 0 and x = 40). The upper (z = 1) and lower (z = 25) boundaries are maintained at the constant temperatures of 2 5 ° C and 625°C respectively. Throughout the remainder of the area, heat is transferred by conduction alone. The physical rationalisation of these conditions is discussed below. The material parameters of the thermal model are less well constrained. Uncertainties in the thermal diffusivities are discussed by England and Thompson [4]. I have used their intermediate value of 28 km2/Ma for the calculations and recognise that possible vertical gradients in conductivity, heat capacity and density may substantially alter the absolute temperature values in the calculations. However, the relative thermal trends remain valid for a wide choice of material parameters. Maintenance of the upper and lower boundaries at constant temperatures implies the assumption of definite cooling and material transfer mechanisms. In the upper layer, maintenance of the thermal boundary (z = 1) at 25°C implies that removal of heat at that boundary is accomplished by convection through access of meteoric waters to the hot rock (see also [17]). The presence of hot springs within the Southern Alps is evidence for both free and forced convective flow and has been modeled as such [26]. Structural evidence coupled with fluid inclusion studies of the quartz- and gold-bearing veins in the schist adjacent to the Alpine Fault indicates that the convective systems associated with uplift are restricted to the upper 5 km of the crust [27]. The assumption of constant temperature at the lower boundary is more controversial as it implies that as the rocks are uplifted from a defined depth they are replaced horizontally from the east by crustal material of the same thermal structure. This tectonic configuration, which is validated by plate tectonic evidence (e.g. [18,19]), requires that

the lower boundary condition be one of constant temperature. Constant flux, a more commonly used boundary condition in thermal modeling, is considered inappropriate for modeling the observed tectonic regime. Although I have attempted to avoid discussion of the dynamics of collision in this study, the assumption of a constant temperature lower boundary requires the existence of a d6collement zone and therefore a delamination within the crust of the upward bound material and the remainder of the lithosphere. The extensive occurrence of mylonites adjacent to the Alpine Fault is evidence for concentration of high-strain deformation at depths below the brittle-ductile transition [28,25]. However, the spatial relationships of the mylonites to the rest of the crust are not clear. On the basis of petrologic arguments [23] that approximately 25 km of crust is currently exposed along the Alpine Fault, I have chosen 25 km as the depth of the d6collement. This depth is clearly open to argument. The presence of crustal delamination was suggested previously by Wellman [22]. Two reasonable alternatives to this d6collement are: (1) the crust may be interlayered with crustal packets from the Australian plate [29], or (2) the uplifted material may be replaced vertically by deeper material. The first alternative would lead to cooling through interlayering of Australian and Pacific crust and because of the distinctly different rock types of the two plates in this region should be obvious in the local geology. One of the marked characteristics of the Haast Schist near the Alpine Fault is the monotonous lithological association of quartzofeldspathic schists and greenschists with occasional ultramafic pods and the absence of structural units from other than the Pacific plate (e.g. [30]). The geological evidence therefore precludes interlayering within at least t h e - 25 km of crust currently exposed. The second alternative, that of vertical replacement of the uplifted crust by deeper crust and mantle material, would not yet be evident in the surface geology, but the possibility of near-surface mantle material is excluded by the presence in this region of a large negative bouguer and isostatic anomaly [31]. Shear heating is not considered as a heat source for the numerical modeling. Field and experimental investigations [32,33] demonstrate that, even

311

with the uncertainties regarding the efficiency of shear heating, the strain rates and stresses discussed below would result in significant transient

heating adjacent to the Alpine Fault, but the duration and magnitude of a transient, shear heating pulse is very poorly constrained. The omission

A 5

10

15

~,,,,,j,,,~,

D 20

25

, , , , ~ , , , ~ , , ,

50

5`5

40

5

, , ,,,,,,,,,,~

10

1S

20

25

~

50

5S I

' ''

40 200

'~'''

:OoO

10

~

10

.

400

20

5

_

..--4uu

0

0

~

`

5

0

0

~

~

_=0

~

~

-2C

B , , ,5

-

10

1,5

~

E ,

20

25

.50

, ,5,`5

,,

4.0

S

10

15

20

25

SO

55

40200

-10

0

-so

,~

-as~

'

~

~

' 'S"

as' ~ -

ss~%7-=~

-~o

~,

_

-~s

C 5

.

.

10

.

15

.

.

4

F 20

.

~

25

.

.

30

.

35

.

.

S

40

10

15

2°°

20 - _

25

~

.

-

30 35 , 40 ' '~-~--~-q ~oo

~

~

~0 4

10

-2`5

....

~''

'10''

lb

''~0

25

.Su

6b

40

Fig. 2. T w o - d i m e n s i o n a l c o n t o u r d i a g r a m s of t e m p e r a t u r e in the p r o b l e m d o m a i n f o r v a r i o u s times f r o m t = 0 to t = 4.0 M a . C o n t o u r i n t e r v a l is 20 o C a n d the s p a t i a l d i m e n s i o n s a r e in k i l o m e t r e s . T h e g e o m e t r y in all o f the f i g u r e s is t h e s a m e w i t h a vertical d i s c o n t i n u i t y a t x = 10 kin. T h e b l o c k to the w e s t o f this d i s c o n t i n u i t y r e m a i n s s t a t i o n a r y w i t h n o uplift. S u p e r i m p o s e d o n the t w o - d i m e n s i o n a l t h e r m a l p a t t e r n s a r e the 2 0 0 ° C, 300 o C, 4 0 0 o C a n d 500 o C i s o t h e r m s c a l c u l a t e d f r o m a o n e - d i m e n s i o n a l s o l u t i o n a t a c o n s t a n t u p l i f t r a t e o f 10 m m / a . T h e m a x i m u m u p l i f t a t e a c h o f the t i m e s is: (A) t = 0, u p l i f t = 0 k m ; (B) t = 1.0 M a , uplift = 10 k m ; (C) t = 2.0 M a , u p l i f t = 2 0 k m ; ( D ) t = 2.5 M a , uplift = 25 kin; (E) t = 3.0 M a , uplift = 30 k m ; (F) t = 4.0 M a , uplift = 4 0 k m .

312

of shear heating as a heat source causes the anomalous temperatures calculated below to be due solely to uplift and therefore probably minim u m temperatures.

3. Results of the thermal modeling The results of the two-dimensional thermal modeling are presented in Fig. 2 for times varying from t = 0 (initial thermal state) to t = 4.0 Ma. Contours are at 20 o C intervals over the 40 km by 25 km x/z plane. Superimposed on the two-dimensional thermal profiles of Fig. 2 are selected isotherms from a one-dimensional solution with an uplift rate of 10 m m / a . At t = 0 the undisturbed thermal profile has a constant geothermal gradient of 25 o C/kin. After 1.0 Ma and 10 km of uplift, the pattern of a high geothermal gradient in the upper 6 km and spreading of the isotherms at depth is well established. The temperature at 5 km depth within the region of m a x i m u m uplift has been elevated from an initial value of 1 2 5 ° C to >/270 o C. The 500 ° C isotherm has risen from 20 km at t = 0 to < 14 km. The thermal highs at this stage correspond clearly to the regions of high uplift with little thermal diffusion into the western block. Continued uplift until 2.0 Ma accentuates the earlier established pattern with the maximum temperature at 5 km increasing to - 340 ° C and the 400 ° C isotherm rising above 7 kin. In addition, two-dimensional broadening of the thermal anomaly in the upper 2 - 3 km has occurred, producing flattened isotherms within this upper region. At 2.5 Ma, maximum total uplift has reached 25 kin, the temperature at 5 km is in excess of 350 ° C and the 400 ° C isotherm has risen - 6 kin. The geothermal gradient within the thermally compressed upper 3 km follows an error function shape which, over the upper 3 km, averages 100 o C / k m . Thermal diffusion has resulted in further widening of the thermal anomaly into the western segment of the problem domain as well as broadening of the high-temperature region to those areas of lower total uplift. At t > 2.5 Ma, the vertical thermal pattern in the region of most rapid uplift does not alter a great deal from that of t = 2.5 Ma. Lateral thermal changes do occur with continued widening of the thermal anomaly. Fig. 2F shows the tempera-

tures at 4.0 Ma in which the thermal anomaly has spread laterally into the regions which have not undergone uplift. Variations in the thermal diffusivity used in the numerical solutions from 38 k m 2 / M a to 19 k m 2 / M a alters the upper thermal structure by approximately one km on either side of the values presented above for K = 28 km2/Ma. This corresponds to increasing the effects of the advective thermal component if the lower diffusivities are used or increasing the lateral conduction of heat if higher diffusivities are used. The impact of changes in the geometry on the resultant thermal structure was investigated in order to evaluate the importance of the lateral diffusion term. A decrease in the width of the rapidly uplifted region from 5 km to 1 km lowers the isotherms by somewhat more than 1 km due to the increased relative importance of lateral heat conduction. Alternatively, any increase in the lateral dimensions of the high uplift zone from the initial 5 km width diminishes the relative importance of lateral heat loss. A high uplift zone of width greater than 15 km causes the thermal anomaly to approach that of a one-dimensional model [17] for which analytical solutions exist (e.g. [34]). The 200 o C, 300 o C, 400 o C, 500 o C isotherms for a one-dimensional solution are plotted together with the two-dimensional isotherms in Fig. 2 for each of the time steps. The one-dimensional thermal pattern exhibits little variation for t > 1.8 Ma. This compares with the two-dimensional vertical pattern which varies little after 2.5 Ma. The greatest overall discrepancies between the one-dimensional and two-dimensional patterns lie within this time period of - 1.5 Ma ~< t ~< 2.5 Ma. At all times, the 2 0 0 ° C isotherms for both solutions lie within 1 km of each other, while the greatest variation occurs in the positions of the 400 ° C and 500 ° C isotherms at 2.0 Ma ~< t ~< 3.5 Ma. In this time span, the one-dimensional isotherms are displaced upwards from the positions of the two-dimensional maxima by 2 and 3 km, respectively. The comparison of the one-dimensional and two-dimensional isotherm patterns in Fig. 2 indicates that although the one-dimensional solution provides a good approximation for the near-surface isotherms of~< 250°C, the one-dimensional thermal pattern is not reliable for the positioning of isotherms of >f 300 ° C. The one-di-

313

D

A 5

10

15

~0

25

30

35

5 10 15 [ ' ' ' Lg._rl ' ' '2 . . . . . . . . . . . . .

'40

-5 -10

-10

-15

-15

10

1S

25 30 35 P~Nt . . . . . . . . . . . .

4-0

I

-5

10

15

_20 L

20

-2sS

25

B 5

20

£ 20

25

~0

55

40

5

10

15

20

25

30

3S

40

-S

-5

-10

-10

-15

-15

20

20

-25

-25

C 5

10

15

20

25

50

55

S

40

10

15

~'0 . . . .

~

20

25

50

.55

40

-5

i

15 -20

10

-10

15

-IS

_~0

_20 I

_~5

-2s:

i

,,

,~ ....

'

s

'

~

....

~'s' ' ~ 4

Fig. 3. Two-dimensional contour diagrams of maximum differential stress within the crust derived from the rheological equations of Brace and Kohlstedt [36] for the thermal distributions shown in Fig. 2. The contour interval is 20 MPa. The transition from brittle to ductile behaviour is marked by the steep negative stress gradient. The stepwise roughness of some contours near steep gradients is due to a contouring algorithm deficiency. The times and maximum uplift values for each of the figures correspond to those of Fig. 2.

mensional solution, while useful for describing the thermal behaviour of broad ( > 20 km) regions of uplift, does not yield any information on the thermal behaviour within the adjacent non-up-

lifted material. Knowledge of the thermal evolution of the region outside of the uplifted block is essential for the following discussion on the mechanical behaviour of a collision zone.

314

4. Effects of rapid uplift on crustal stress levels and lithospheric strength Perturbation of the initial steady state thermal profile due to sustained rapid uplift results in significant changes in the physical properties of the crustal material in and adjacent to the two-dimensional region of uplift. The strength of the crust is primarily a function of the thickness of the brittle upper layer [35] and therefore any processes which thin or thicken this layer will effect the overall crustal strength. Although the importance of crustal composition has been discussed (e.g. [10]) and may be of local significance, the thermal structure of the crust is the major factor determining the deformation mechanism of a homogeneous quartzofeldspathic crust. In a material in which the stength is dominantly a function of temperature, the above thermal perturbation causes redistribution of stress and strain throughout the lithosphere. Fig. 3 shows contour diagrams of the maximum differential stress sustainable within the crust at the same starting conditions and uplift rates as those used to derive the thermal profiles of Fig. 2. The contours of maximum differential stress were derived for a constant strain rate (5 × 10 -15 s -1 for comparison of results with work by Sonder and England [35]) using the rheological equations of Brace and Kohlstedt [36] for brittle failure of material with Pnuid = Phydrostatic and ductile failure of crustal material dominated by a quartz rheology. In the ductile regions of the crust, the strain rate is a function of stress and a constant value of strain rate is inappropriate for a rigorous description of stress-strain relationships with varying stress levels. However, the depth to the brittle-ductile transition is not particularly sensitive to strain rate variations of several orders of magnitude. The role of strain focussing in time and space due to non-homogeneous behaviour of the crust is probably very important yet outside the scope of this paper. The consequences of horizontally and vertically varying stress patterns are discussed qualitatively below. The literature abounds with differing values of maximum differential crustal stress levels ranging from experimentally derived values of several hundreds of MPa [36], through heatflow estimates of tens of MPa to values estimated through seismic stress drops of << 10 MPa (e.g. [37]). While recognising

that the absolute values of stress may be in error, the constitutive equations are used with the hope that crustal materials deform in a manner which is at least qualitatively described by these equations (see also [38]), The differential stress values thus derived are likely to be maximum values. The initial differential stress values for an unperturbed crust are shown in Fig. 3A with a contour interval of 20 MPa. As noted by several previous authors (e.g. [10,36,39]) in a homogeneous crust, maximum crustal stresses are developed immediately above the transition between brittle and ductile behaviour. The strength of the crust then decays rapidly with depth below this transition. With a geothermal gradient of 25 ° C, the brittle-ductile transition occurs a t - 13 km depth with a differential stress maximum of - 650 MPa. Because the depth to the brittle-ductile transition is primarily a function of temperature, the two-dimensional stress maxima patterns due to uplift assume the shape of the modified error function of the thermal curves and are thoroughly predictable from the thermal profiles. As in the thermal profiles, the most rapid decrease in crustal strength occurs in the early stages of uplift. 1 Ma of uplift results in a shallowing of the brittle-ductile transition from - 13 km to < 8 km and reduces the maximum differential stress in the weakest region b y - 30% from 650 MPa to 450 MPa. After 1.5 Ma of uplift the brittle layer is less than 6 km thick at its thinnest point. The brittle layer continues to thin to < 5 km at t ~< 2.5 Ma and is capable of maintaining a maximum differential stress level of < 275 MPa. Further uplift does not cause greater thinning of the brittle layer, but the low strength zone does increase in size by lateral spreading. The width of the low strength zone at 5 km depth increases from 4 km at t = 2.5 Ma to 9 km at t = 4.0 Ma. Although the majority of this lateral weakening occurs within the uplifted block, weakening at depth due to lateral conduction into the stationary block also occurs, reinforcing the two-dimensional character of the anomaly (Fig. 3F), Variations in the thermal constants and geometries used in the numerical modeling as discussed above, alter the temperature-dependent depth of the brittle-ductile transition by a corresponding, minor amount. Change in the uplift rate does alter the transition depth. Solutions were obtained for

315 K = 28 k m 2 / M a and the standard geometry with the m a x i m u m uplift rates of 5 m m / a and 2 m m / a . At an uplift rate of 5 m m / a the brittle-ductile transition rises substantially from about 13 k m to less than 8 km after 3.0 Ma and 15 km of uplift. Uplift rates of 2 m m / a have little impact on the overall thermal or mechanical structure of the upper crust at times up to 4.0 Ma which corresponds to 8 km of total uplift. The numerical solutions were not carried beyond 4.0 Ma. 5. Effects of thermal weakening in the uplifted zone on the overall stress and strain patterns of the collision zone The consequences of a time-variant crustal weakening due to uplift are important for the evolution of the consequent mountain belt and for the stress-strain relationships within the two- and three-dimensional collision zone. In the initial stages of convergence across a fault, a similar symmetrical stress field exists on either side of the fault. The lateral extent of this idealised stress field is determined by diffusion of stress away from the fault zone and by heterogeneities within the colliding blocks. Strain is accomodated by brittle failure within the elastic crustal layer on both sides of the fault and the zone of deformation will be diffuse. Following the initiation of uplift on one side of the fault, weakening of the crust by thinning of the elastic layer will cause concentration of strain within this weakened zone adjacent to the fault. Deformation on the side of the fault which has undergone no uplift may be expected to decay with time as strain is transferred laterally to the thermally weakened region. The stationary block remains a relatively highstrength region and as such assumes the characteristics of an idealised rigid indentor. As the maxim u m sustainable stress within the uplifted zone decreases with time due to the decreasing depth of the brittle-ductile transition, the overall stress levels within the adjacent crust will also decrease but to a lesser extent than the stress decrease within the weakened zone. This will lead to a horizontal zone of high-strain concentrated within the thermally weakened region and a horizontal zone of relatively high stress within the adjacent, unweakened crust. In addition to lateral stress transfer, the devel-

opment of a two-dimensional zone of low strength within the crust causes transferral of stress vertically to the high strength upper mantle. Increased stresses in the upper mantle may result in brittle failure within this high-strength region and thereby generate intermediate-depth seismicity. The lateral extent of this zone of intermediate-depth seismicity will be related to the mechanism by which convergent strain is accommodated within the lower crust and upper mantle. The mechanism for strain within this lower crust-upper mantle region is not considered in this paper, but underplating by crust-mantle packets of 12-20 km thickness similar to those described in the Himalaya [40] has been suggested for the New Zealand collision zone [41,421. Evolution of the stress regime as described above has implications for the development of seismicity patterns with time. Meissner and Strehlau [39] and Sibson [10] examined the focal depths for a number of intraplate earthquakes in regions of differing surface heat flow and related these focal depths to the depth of the brittle-ductile transition. Both studies have shown that the depth to the seismogenic zone [10] is a function of heat flow. Meissner and Strehlau [39] also suggested that, in a very general way, shallower earthquakes in a given region tended to have smaller magnitudes than the deeper earthquakes, reflecting the lower differential stress levels sustainable at higher layers within the seismogenic zone. F r o m these observations, it follows that a shallowing of the brittle-ductile transition due to rapid uplift will cause a subsequent shallowing of the seismogenic zone and a consequent reduction in the seismic moment of earthquakes generated within the uplifted zone. Scholz et al. [18] have interpreted the low seismicity and the absence of any earthquake greater than magnitude 7 along the Alpine Fault in recorded times to be indicative of the accumulation of large elastic strain. However, others [17,19] have demonstrated that the strain along the plate boundary is likely to be taken up by ductile rather than elastic deformation. The stress profiles developed above (Fig. 3) reinforce the ductile deformation hypothesis and indicate that, a f t e r - 1.5 Ma of uplift at 10 m m / a , the elastic thickness of the crust is less than 5 km. On the evidence of displaced alluvial surfaces along the

316 Alpine Fault, Adams [20] predicted that significant elastic strain is accumulating and that the lack of large earthquakes in historical times is a function of the brevity of the historical record and not of ductile deformation. Unfortunately, it is in the prediction of actual seismicity that our ignorance concerning the earthquake generation process and the absolute levels of crustal stress becomes most obvious. Recognising that stress drops of major earthquakes are of the order of less than 10 MPa [37] and that these stress levels are well below those predicted above even for elastic thicknesses of less than 3 km, it remains impossible to choose between the two cases presented above of stored elastic strain versus continual ductile strain release. However, given the reduction in the vertical dimension of the seismogenic zone in the thermally weakened zone it is unlikely that earthquakes of large seismic moment will occur within this region. Chen and Molnar [43] have demonstrated through compilation of focal depths for intra-continental earthquakes that sub-crustal, intermediate depth seismicity is widespread in active convergence zones. They assigned the origin of the seismicity to failure within the high-strength upper mantle. As discussed above, thermal erosion of the high-strength brittle layer of the crust results in transferral of stress vertically as well as horizontally. Therefore, convergence zones in which the crust has been thermally weakened by rapid uplift may be expected to exhibit intermediate-depth seismicity due to failure within the upper mantle. Sub-crustal seismicity has been described within the New Zealand convergence zone [44-46] and may be due to this mechanism of vertical stress transfer. In a study of the microseisms of the central South Island, Reyners [46] has described both a westward dipping zone of intermediatedepth earthquakes deepening towards the Alpine Fault and a zone of crustal earthquakes which shallows towards the Alpine Fault.

models have therefore produced a generally cooler lithosphere than existed prior to collision. The boundary condition in the present model of constant temperature and and hence d~collement within the lower crust causes some volume problems within the material not uplifted. Equivalent material in the Himalaya appears to be underplated in crust-mantle packets o f - 1 2 - 2 0 km thickness [40] and there is some seismic and gravity evidence of a similar mechanism occurring in New Zealand [41,42]. The thermal structure within a region of underthrust packets of lower crust-upper mantle in which uplift occurs as a response to isostatic rebound following the major tectonic event is the same as that calculated for most collision belts (e.g. [3]). This underthrusting and related cooling of the middle lithosphere which accompanies rapid uplift in collision zones has the effect of thickening the high-strength layer within the mantle. Sonder and England [35], using a thin plate approximation for the lithosphere and the constituitive equations of Brace and Kohlstedt [36], have produced estimates of vertically integrated lithospheric strength under varying pore pressures, Moho temperatures and stress exponents. In their calculations, the contribution of the crust to the total lithospheric strength is dominantly a function of the depth to the brittle-ductile transition. A reduction in the thickness of the brittle layer from 20 km to effectively zero thickness (X = 1) results, in their models, in a reduction of total lithospheric strength of approximately 1.5 orders of magnitude at a strain rate of 0.5 x 10 -15 [35, fig. 1]. It must be emphasised that the calculations of Sonder and England [35] are based upon the thin plate approximation and are not intended to account for the contribution to the total lithospheric strength due to thickening of the highstrength upper mantle member during collision. Rapid uplift will thin the upper crust high strength layer but will increase the thickness of the lower high-strength layer.

6. Thermal and mechanical behaviour of the remaining lithosphere

7. Other collision zones

Most other heat flow models of continental collision have been concerned with crustal and lithospheric thickening through displacement downwards of hotter material by cooler. These

The New Zealand example of convergence is particularly simple both in geometry and in timing relative to other collision zones. However, some thermal and mechanical aspects of other collision

317

zones may also be the result of sustained rapid uplift. Uplift rates have been estimated on the basis of geological observations [47] to be >/5 m m / a within the Himalaya, with the most rapid uplift rates near the highest elevations. Molnar [48] briefly reviewed uplift rates within the Himalayan collision zone and concluded that the high uplift rates ( > 5 m m / a ) observed in the Greater and Lesser Himalaya pertained to those regions in which they were measured and not to the entire collision zone. Molnar et al. [6] also suggest that the majority of uplift in the Greater Himalaya is less than 10 and perhaps less than 5 Ma old. With the exception of one recent study [49] which delineated a region of high heat flow south of the Yalu-Zangbo Suture, no rigorous heat flow measurements are available for the collision zone. However, thermal activity throughout the collision zone in the form of hot springs has been documented [5,50] indicating high crustal heat flow. Francheteau et al. [49] and Jaupart and Francheteau [51] have suggested that a thermal anomaly measured in the Tibetan lakes is the result of recent igneous intrusion of granitic material. Their hypothesis is entirely compatible with local geology and the uplift rates which they use in their erosion calculation. These rates were derived by averaging total regional uplift over 35 Ma yielding uplift rates of 0.2 m m / a [52]. The heat flow pattern ascribed [51] to igneous activity is similar to the pattern which would be observed in a region of generally low uplift rates which contained discrete zones of rapid uplift. LyonCaen and Molnar [53] have suggested that Himalayan uplift has occurred within the last 10 and perhaps the last 5 Ma. This short duration of uplift would produce uplift rates which are compatible with those observed in other mountainous regions ( - 10 m m / a ) which in turn would result in regions ot anomalously high heat flow. Molnar arid Tapponnier [38] discussed the distribution of active strain within the Himalayan zone and showed that the differences in the degree of tectonic activity were related to the relative ages of orogenic activity. They ascribed this to a thermal weakening of the crust due to an increased mantle heat flux in those areas of active tectonism. The thermal and differential stress profiles produced during sustained rapid uplift may provide

an alternative mechanism for lithospheric weakening and thereby concentration of strain within the weakened zones. 8. Conclusions In continental collision, rapid uplift of crustal blocks will produce anomalously hot regions within the crust and consequent surface heat flow anomalies. Thermal erosion of the brittle upper crustal layer to less than 5 km results in significant reduction of crustal strength. The zone of crustal weakening becomes a high-strain zone while stress is transferred horizontally to the stationary highstrength crust and vertically to the high-strength region within the upper mantle. Seismicity patterns within and adjacent to the rapid uplift zone are predicted to change with time as weakening of the upper crust progresses. A decrease in the depth and frequency of large earthquakes is expected within the weakened crust. Intermediate depth seismicity which has been observed in regions of continental collision [43-46] is predicted in this model due to the vertical transferral of stress from the crust to the upper mantle. Although the crust is weakened by uplift, the effect of rapid uplift on the total lithospheric strength is not obvious due to the displacement downwards of crustal material from below the crustal drcollement. The drcollement necessitated by the sustained rapid uplift in regions of horizontal convergence implies discontinuities of strain rates as well as strain rate gradients with depth. This is completely compatible with geological evidence gathered from nappe zones throughout the world in which large variations in the degree of vertical coupling may be observed. The thermal modeling presented in this paper is based upon observed uplift rates which are consistent w~th general rates of tectonic processes involved in collision. The thermal anomalies generated by rapid uplift are a consequence of tectonic processes occurring much more rapidly than does dissipation of heat by conduction. Because of the thermal inertia of lithospheric materials, the thermal and mechanical behaviour of the lithosphere discussed in this paper in the context of continental collision may be expected to be produced in any rapid tectonic event. These may include extensional as well as compressional processes.

318

The rapid uplift involved in the early mountain building phase of collision is a potent heat engine and as such is capable of exerting strong influence on thermal, mechanical and geochemical systems operating within the lithosphere. Evidence from active collision zones indicates that the thermal and mechanical effects of tectonic processes with timescales of less than 10 Ma may determine lithospheric behaviour and yet not be easily visible in the geological record comprising hundreds of millions of years.

Acknowledgements Constant discussions and wranglings with my colleagues and students at the University of Otago have caused me to attempt to quantify the concepts in this paper. I am particularly indebted to Richard Norris, Alan Cooper and Dave Craw for ideas as well as comments on the manuscript. The suggestions and criticisms of C. Beaumont, P. England and an anonymous reviewer are appreciated. This is not to imply that they share any of manuscript's failings. The University of Otago Research Committee provided all funding for the project. References 1 P. Morgan, The thermal structure and thermal evolution of the continental lithosphere, Phys. Chem. Earth 15, 107-196, 1984. 2 E.R. Oxburgh and D.L. Turcotte, Thermal gradients and regional metamorphism in overthrust terrains with special reference to the Eastern Alps, Schweiz. Mineral. Petrogr. Mitt. 54, 641-662, 1974. 3 P.C. England and S.W. Richardson, Erosion and the age dependence of continental heat flow, Geophys. J.R. Astron. Soc. 62, 421-437, 1980. 4 P.C. England and A.B. Thompson, Pressure-temperaturetime paths of regional metamorphism, I. Heat transfer during the evolution of regions of thickened continental crust, J. Petrol. 25, 894-928, 1984. 5 W.R. Buck and M.N. Toks~Sz, Thermal effects of continental collision. Thickening a variable viscosity lithosphere, Tectonophysics 100, 53-69, 1983. 6 P. Molnar, W.-P. Chen, and E. Padovani, Temperatures in overthrust terrains and the origin of Tertiary granites in the greater Himalaya, J. Geophys. Res. 88, 6415-6429, 1983. 7 M.N. Tokst~z and P. Bird, Modelling of temperatures in continental convergence zones, Tectonophysics 41,181-193, 1977. 8 R. Stephenson and K. Lambeck, Erosion-isostatic rebound models for uplift: an application to south-eastern Australia, Geophys. J.R. Astron. Soc. 82, 31-55, 1985.

9 A.B. Thompson and P.C. England, Pressure-temperaturetime paths of regional metamorphism II. Their inference and interpretation using mineral assemblages in metamorphic rocks, J. Petrol. 25,929-955, 1984. 10 R.H, Sibson, Roughness at the base of the Seismogenic Zone, contributing factors, J. Geophys. Res. 89, 5791-5799, 1982. 11 N.J. Kusznir and R.G. Park, Intraplate lithosphere deformation and the strength of the lithosphere, Geophys. J.R. Astron. Soc. 79, 513-538, 1984. 12 P.K. Zeitler, N.M. Johnson, C.W. Naeser and R.A.K. Tahirkheli, Fission-track evidence for Quaternary uplift of the Nanga Parbat region, Pakistan, Nature 298, 255-257, 1982. 13 W.B. Bull and A.F. Cooper, Marine terraces along the Alpine Fault, Science, in press, 1986. 14 G. Plafker and M. Rubin, Uplift history and ear',hquake recurrence as deduced from marine terraces on Middleton Island, Alaska: U.S. Geol. Surv. Open-File Rep. 78-943, 687-721, 1978. 15 T.J. Dempster, Uplift patterns and orogenic evolution in the Scottish Dalradian, J. Geol. Soc. London 142, 111-128, 1985. 16 E.A. Silver, M.J. Ellis and N.A. Breen, Comments on the growth of accretionary wedges, Geology 13, 6-9, 1985. 17 R.G. Allis, R.W. Henley and A.F. Carman, The thermal regime beneath the Southern Alps, R. Soc. N.Z. Bull. 18, 79-85, 1979. 18 C.H. Scholz, J.M.W. Rynn, R.W. Weed and C. Frohlich, Detailed seismicity of the Alpine Fault Zone and Fiordland region, New Zealand, Geophys. Res. 84, 6770-6982, 1973. 19 R.I. Walcott, Present tectonics and late Cenozoic evolution of New Zealand, Geophys. J.R. Astron. Soc. 52, 137-164, 1978. 20 J. Adams, in: Origin of the Southern Alps, R.I. Walcott and M.M. Cresswell, eds., Bull. R. Soc. N.Z. 18, 47-54, 1979. 21 G.H. Brick, S.A.L. Read and P.T. Hall, Ostler Fault Zone. Progress report on surveying results 1966-1985, levelling and tilt levelling, N.Z. Geol. Surv. Rep. EDS 103, 1985. 22 W.H. Wellman, An uplift map for the South Island of New Zealand, and a model for uplift of the Southern Alps, in: Origin of the Southern Alps, R.I. Walcott and M.M. Cresswell, eds. R. Soc. N.Z. Bull. 18, 13-20, 1979. 23 A.F. Cooper, Retrograde alteration of chromian kyanite in metachert and amphibolite whiteschist from the Southern Alps, New Zealand with implications for uplift on the Alpine Fault, Contrib. Mineral. Petrol. 75, 153-164, 1980. 24 S. Nathan, H.J. Anderson, R.A. Cook, R.H. Herzer, R.H. Hoskins, J.I. Raine and D. Smale, Cretaceous and Cenozoic sedimentary basins of the West Coast region, South Island, New Zealand, in: N.Z. Geological Survey Basin Studies, I, DSIR, Wellington, 1986. 25 M.S. Rattenbury, Late low-angle thrusting and the Alpine Fault, central Westland, New Zealand, N.Z.J. Geol. Geophys. 29, 437-446, 1986. 26 I.G. Donaldson, Temperature gradients in the upper layers of the Earth's crust due to convective water flow, J. Geophys. Res. 67, 3449-3459, 1962. 27 D. Craw, M.S. Rattenbury and R.D. Johnstone, Structural

319

28

29 30

31

32

33 34

35

36

37

38

39

geology and vein mineralisation in the Callery River Headwaters, Southern Alps, New Zealand, N.Z.J. Geol. Geophys. 30, in press, 1987. R.H. Sibson, S.H. White and B.K. Atkinson, Fault rock distribution and structure within the Alpine Fault Zone: a preliminary account, in: Origin of the Southern Alps, R.I. Walcott and M.M. Cresswell, eds., R. Soc. N.Z. Bull. 18, 55-65, 1979. R.G. Allis, Continental underthrusting beneath the Southern Alps of New Zealand, Geology 9, 303-307, 1981. R.H. Findlay, Summary of structural geology of Haast Schist terrane, central Southern Alps, New Zealand: implications of structures for uplift of and deformation within Southern Alps, in: Origin of the Southern Alps, R.I. Walcott and M.M. Cresswell, eds., R. Soc. N.Z. Bull. 18, 113-120, 1979. D.J. Woodward, The crustal structure of the Southern Alps, New Zealand, as determined by gravity, in: Origin of the Southern Alps, R.I. Walcott, and M.M. Cresswell, eds., R. Soc. N.Z. Bull. 18, 95-98, 1979. R.H. Sibson, Generation of pseudotachylyte by ancient seismic faulting, Geophys. J. R. Astron. Soc. 43, 775-795, 1975. D.A. Lockner and P.G. Okubo, Measurements of frictional heating in granite, J. Geophys. Res. 88, 4313-4320, 1983. L. Royden and K.V. Hodges, A technique for analyzing the thermal and uplift histories of eroding orogenic belts: a Scandinavian example, J. Geophys. Res. 89, 7091-7096, 1984. L.J. Sonder and P. England, Vertical averages of rheology of the continental lithosphere: relation to thin sheet parameters, Earth Planet. Sci. Lett. 77, 81-90, 1986. W.F. Brace and D.L. Kohlstedt, Limits on lithospheric stress imposed by laboratory experiments, J. Geophys. Res. 85, 6248-6252, 1980. T.C. Hanks, b Values and seismic source models: implications for tectonic stress variations along active crustal fault zones and the estimation of high-frequency strong ground motion, J. Geophys. Res. 84, 2235-2242, 1979. P. Molnar and P. Tapponnier, A possible dependence of tectonic strength on the age of the crust in Asia, Earth Planet. Sci. Lett. 52, 107-114, 1982. R. Meissner and J. Strehlau, Limits of stresses in continental crusts and their relation to the depth-frequency distribution of shallow earthquakes, Tectonics 1, 73-89, 1982.

40 A. Him, G. Poupinet, G. Wittlinger, X.Z. Zin, G.E. Yuan, T.J. Wen, X.S. Bai, M.R. Pandey and J.M. Tater, Crustal structure and variability of the Himalayan border of Tibet, Nature 307, 23-25, 1984. 41 P.O. Koons, Implications of continental collision in southern New Zealand, Proc. N.Z. Geol. Soc. Annu. Conf., 1984. 42 R.G. Allis, Mode of crustal shortening adjacent to the Alpine Fault, New Zealand. Tectonics 5, 15-32, 1986. 43 W.-P. Chen and P. Molnar, Focal depths of intracontinental and intraplate earthquakes and their implications for the thermal and mechanical properties of the crust, J. Geophys. Res. 88, 4183-4214, 1983. 44 I.M. Calhaem, A.J. Haines and M.A. Lowry, An intermediate depth earthquake in the central region of the South Island used to determine a local crustal thickness, N.Z.J. Geol. Geophys. 20, 353-361, 1977. 45 M. Reyners, C.E. Ingham and B.G. Ferris, Microseismicity of the upper Clutha valley, South Island, New Zealand, N.Z.J. Geol. Geophys. 26, 1-6, 1983. 46 M. Reyners, Subcrnstal earthquakes in the central South Island, New Zealand, and the root of the Southern Alps, Geology, submitted, 1987. 47 A. Gansser, The morphogenic phase of mountain building, in: Mountain Building Processes, K.J. Hsu, ed., pp. 221-222, Academic Press, London, 1982. 48 P. Molnar, Structure and tectonics of the Himalaya: constraints and implications of geophysical data, Annu. Rev. Earth Planet. Sci. 12, 489-518, 1983. 49 J. Francheteau, C. Jaupart, X-J. Shen, W.-H. Kang, D.-L. Lee, J.-C. Bai, H.-P. Wei and H.-Y. Deng, High heat flow in southern Tibet, Nature 307, 32-36, 1984. 50 D.R. Bhattarai, Some geothermal springs of Nepal, Tectonophysics 62, 7-11, 1980. 51 C. Jaupart and J. Francheteau, On the thermal structure of the southern Tibetan crust, Geophys. J.R. Astron. Soc. 81, 131-155, 1985. 52 W-L. Zhao and W.J. Morgan, The uplift history of the Tibetan Plateau from a "Fault Block" reconstruction, EOS, Trans. Am. Geophys. Union 64, 860, 1983. 53 H.J. Lyon-Caen and P. Molnar, Constraints on the structure of the Himalaya from an analysis of gravity anomalies and a flexural model of the lithosphere, J. Geophys. Res. 8171-8191, 1983.