A finite-element algorithm for modelling the subsidence and thermal history of extensional basins

.I. Geo&mnmicsVol. 26, No. 1, pp. I-25, 1998 % 1998 Elsevier Science Ltd All rights reserved. Printed in Great Britain PII: SO264-3707(97)00021-5 02643707/98 $19.00+0.00

A FINITE-ELEMENT ALGORITHM FOR MODELLING THE SUBSIDENCE AND THERMAL HISTORY OF EXTENSIONAL BASINS GUOJIANG

LIU* and GIORGIO RANALLI

Department of Earth Sciences and Ottawa-Carleton Ontario, Canada K 1S 5B6

Geoscience Centre, Carleton University, Ottawa,

(Received 15 July 1996; accepted 5 February 1997)

Abstract-We present a finite-element model for the prediction of the subsidence history and thermal evolution of extensional basins. The model takes into account the inhomogeneity of the lithosphere, and allows the space formed by lithospheric extension to be filled with sediments or water or both. Local isostatic compensation is assumed. The stability and accuracy of the algorithm are investigated. Time increments of 0.054.1 Ma are most appropriate for triangular elements of size 13 km2. For an initially old continental lithosphere (thickness: 130 km), thermal equilibrium is reached in approximately 100 Ma for a basin of initial width 100 km with a maximum /?-factor of 2.0. The duration of extension for a given /Hactor affects the heat flow and subsidence. For extension periods of 10 and 20 Ma, faster extension results in less subsidence and more prolonged high heat flow after the cessation of extension. If the basin is filled with sediments (as opposed to water), thermal blanketing effects are significant. The algorithm can be used to complement our understanding of the effects of thermal blanketing, lateral heat conduction, and sediment loading on the subsidence and thermal evolution of extensional basins. Xj 1998 Elsevier Science Ltd. All rights reserved

INTRODUCTION Extensional basin modelling has been rapidly advancing since the publication of the analysis of McKenzie (1978) who postulated instantaneous lithospheric thinning, resulting in initial subsidence followed by thermal relaxation. Subsidence curves constructed from wells or pseudo-wells from seismic sections normally confirm the pattern of rapid initial subsidence, followed by slow thermal subsidence (Royden and Keen, 1980; Sclater et al., 1980; Barton and Wood, 1984; Thorne and Watts, 1989; Kooi and Cloetingh, 1992; Lippard and Liu, 1992; Skogseid et al., 1992; Van der Beek and Cloetingh, 1992; Cloetingh et al., 1993; Cloetingh, 1995).

*Present address:

SpaceBridge

Networks

Corporation,

115 Champlain

Street, Hull, Quebec,

Canada,

JSX 3RI.

Since McKenzie‘s simple quantitative model. several algorithms have been proposed for the simulation of the kinematics. dynamics and thermal evolution of extensional basins. For kinematic models with prescribed velocity fields. the most commonly applied numerical methods are tinite-difference schemes (Stephenson C/ r/l., 19X7; Hermanrud. 1992). Finite difference schemes, however. encounter difficulties wshen the boundaries of computational domains are irregular, or the material properties within the domain change drastically. Finite-element methods. however. can deal with these problems. Most papers describing the dynamics oflithospheric extension adopt finite-clement methods (Braun and Beaumont, 1987; Lynch and Morgan, 1987: Ricke and Mechie. 19X9: ChCry cl/ (11.. 1990: Keen and Boutilier. 1990; Covers and Wortel. 1993; Boutilier and Keen. 1994; Giilkes PI t/l.. 1994). These papers. by and large, deal with the stress and strain fields in extensional basins. The subsidence and thermal histories of ;I sedimentary basin are of crucial importance to the prediction of hydrocarbon generation. migration, entrapment. and maturation. De Bremaecker (1983) was probably the first to describe :I finite-element model for this simulation of subsidence history. thermal evolution and thermal maturation in extensional basins. De Bremaecker’s algorithm assumed uniform. pure shear extension. without incorporation of the radiogenic heat production. Liu (199.3) developed a more flexible finiteclement algorithm for one-dimensional modelling. allowing for both uniform and nonuniform modes of extension. The objective ofthe present paper is to prehent a two-dimensional tinite-element algorithm specifically designed to model the subsidence history and thermal evolution of extensional basins. including an analysis of the numerical stability of the algorithm in the context of spatial and temporal discretization. The algorithm is simple and useful for kinematic models (prescribed velocity field) with local Airy isostasy. The focus is on methodology rather than modelling of actual structures. However. to give some examples, the algorithm is applied to examine the etfect of the magnitude of stretching factor/j and of the duration of extension on subsidence and temperature histories. Subsidence and thermal evolution in the cases of water. sediment. and ‘non-thermal’ sediment loading are compared in order to assess thermal blanketing effects. TI lb M0I)E.l

The two-dimension~il physical model and the corresponding finite-element discrctiation arc shown in Fig. I. In the examples that the lithosphcrc at least sections direction into in the vertical based and pctrophysical as thermal specific radiogenic heat and density. the concerned domain blocks. of these thermal properties (Note. the algorithm deal with domains with propertics The left-most sections in order the boundary of lateral heat flow boundaries. is discretized into elements. of the computational is set at the of the lithosphere the maximum takes The boundary and the other are to during evolution. zone extension pure is at the centre The temperatures and upper arc constant. at the and right is /era.

A

finite-element algorithm and thermal history of extensional basins

3

(4 0

20

40 g m’ E g

6o 80

100

120 -00

-60

-40

-20

0

X-AXIS (km)

40

60

60 (b)

Fig. I. (a) Physical model of the lithosphere; (b) spatial discretization for the finite-element algorithm. rectangle in (b) is divided into two triangular elements. For clarity, these oblique lines are omitted.

Each

We assume local isostasy with the isostatic compensation level at the base of the lithosphere following the conventional definition (McKenzie, 1978). There is no doubt that the lithosphere has finite vertical strength (Cloetingh et al., 1994; Burov and Diament, 1995; Ranalli, 1995; Cloetingh and Burov, 1996). However, often local isostasy is a reasonable approximation to reality (Cabal and Fernandez, 1995; Ayala et al., 1996) especially when the level of necking is at upper-middle crustal depths (Ranalli, 1997).

4

G. Liu and G. Ranalli

The extensional process is divided into an active phase and a subsequent thermal relaxation phase. The active extensional phase is further divided into short time intervals At. Within each At, the extension is instantaneous and is followed by a thermal relaxation period. Assuming that the amount of extension is the same in each time interval (constant strain rate). extension at position .Yfor each time interval is A/$.X) = [I]”

““.

(1)

where /3(r) is the total b-factor at .I-. and PIEis the number of time intervals in the active extension phase. In the case of water loading, the amount of thinning during each time interval is referred to the original lithospheric thickness. In the case of sediment loading. the sediments accumulated during active extension (syn-rift sediments) become part of the thickness of the stretching lithosphere in the next time interval. In the application presented here, the p-factor does not vary in the vertical direction. Therefore, the model applies to pure shear extension. The general energy equation can be expressed as (White, 1991)

(2) where D/Dt is the material derivative, /r the thermal enthalpy defined in eqn. (3). k thermal conductivity,p pressure, p density, S radiogenic heat production per unit volume, U,velocity, s:, the stress tensor defined in eqn (4). and .Y,and t spatial coordinates and time, respectively. Enthalpy and stress for both Newtonian and non-Newtonian fluids are given by

(4) where e is the internal energy per unit mass, and /Land i material properties. For Newtonian fluids. p is the dynamic viscosity and n = I. For non-Newtonian fluids, n c 1 in crystalline silicates under most conditions (Ranalli, 1995). The deformation of Earth materials can be regarded as incompressible, i.e.

iu,

DP

where L’,,is specific heat at constant

pressure.

~ = 0. i\-,

nt = 0,

LIP= c,, d T.

Substituting

(5)

eqns (3))(5) into eqn (2). we have

(6) where

is the viscous dissipation (per unit volume). This term is neglected in most studies. This neglect could cause errors in temperature computations, especially near the centre of

A finite-element algorithm and thermal history of extensional basins

5

extension. The dissipation @ is of the order of ti,.$i where e,, is the strain rate tensor. In pure shear, this can be expressed as 6 = (AB- l)/At.

(8)

Taking At = 1 Ma and Ab = 1.1 (corresponding to an overall fl = 2.59 if m = lo), the strain rate according to eqn (8) is 3.17 x lo-l5 SC’. With 100 MPa as a typical value for t,, (Brace and Kohlstedt, 1980; Burov and Diament, 1995) @ sz 3.17 x lo-’ Wmp3. This is approximately of the same order as the heat diffusion term kAT/L2, where AT and L are the temperature difference between the upper and lower boundaries of the lithosphere and its thickness, respectively. However, @ decays rapidly as stress and strain rate decrease away from the centre of extension. Another simplification is the neglect of the term related to advective heat transfer in eqn (6). The active extension phase is divided into short intervals, within which the extension is instantaneous. Therefore, instead of direct evaluation of the advective terms a(pc0 ‘ax+”

a(PcT)

ay

9

where u and v are the velocities in the horizontal and vertical directions, we account for advection by fixing the temperature at the top and bottom of each element before and after deformation. This idea was first introduced by Odling (1992) and further refined by Liu ( 1993, 1994) for one-dimensional modelling. Dropping the terms of energy advection and viscous dissipation, eqn (6) becomes (9) where

(10) Equation (9) is discretized into finite-element equations. The detailed procedures are provided in Appendix A. After the assembly of the element matrices and vectors, the general form of the finite-element equation is d(T)

‘w~ dt

+[fl{T}-{F}

= 0,

(11)

where [MJ, [a, and {F} are global mass matrix, global stiffness matrix, and global load vector, respectively. The time derivative of the temperature is evaluated by a fully implicit scheme. Therefore, eqn (11) becomes

(12) (13)

6. Liu and G. Radii

6

Computation of initial subsidence is based on local isostasy. with compensation level at the base of the undeformed lithosphere. Thus. the initial subsidence can be expressed as .\,(.Y) = (/?I,, uz, ~~ru,),:[fj ,,,,- p,),(I ~-~XT*)].

(14)

where I)z,, and 01, are the masses per unit area of the lithospheric column before and after thinning, tq is the product of the density of the asthenospheric mantle and the difference in lithospheric thickness before and after thinning. /Jr,, and p,,, are the reference density of respectively, CAis the thermal expansion infilled sediments and of the asthenospherc. coefficient. and T* is the temperature at the base of the lithosphere. These terms can be written as

111,) =

I ,,J I - z 7‘(.\-..I.. I = 0)] d.1,.

1 I I

1~1~ =

p,,,(

.,

I

rT*)L(t)(

I - I ‘A/i).

(15)

(17)

where J’,+ I and J’, are the top and bottom of layer i. p, is the density of the ith lithospheric layer, and L(t) is the lithospheric thickness at time [. Thermal subsidence caused by thermal cooling between time t and I + AI is calculated by

A.\, =

The heat flow at position

.\ is calculated f/(f) =

by Fourier’s _.

I,, -A,! ,’ 8, I

where hi,, is the thermal conductivity of the topmost J‘,~+ I. J‘,, are temperatures and vertical coordinates layer.

NUMEKICAL

Law of heat conduction

I;,

,’ .

(19)

,/

layer of the lithosphere. and T,,+ 1. T,,. at the top and the bottom of the top

C‘ONSIDERATIONS

The convergence, stability. and accuracy of a numerical model are areas of concern. The fully implicit scheme for the discretization of the temperature derivative with respect to time is unconditionally stable (Press CJ~N/., 1992). However, the accuracy of numerical results can be affected by the time step At and the size of the elements. We have tested the

A finite-element algorithm and thermal history of extensional basins 1. Thermal

Layer no.

Depth (km)

k

c (10’)

p (102)

s

C!(10 -‘)

Sediments 1 2 3 4

&variable Variable-l 0 lo- 20 20-30 30-120

2.0 2.0 2.5 3.0 4.0

1.0 1.0 1.o 1.0

2.1 2.1 2.8 2.9 3.3

0.8 0.8 0.5 0.1 0

3.28 3.28 3.28 3.28 3.28

I .o

Within each layer, the properties are assumed to be constant. Symbols are as follows: k, thermal conductivity (Wrn~-’ K ‘), c: specific heat (J kg-’ Km’); p. density (kgm- ‘); S, radiogenic heat production [~Wrn- ‘I; and c(, thermal expansion coefficient (K-l).

accuracy of the procedure by comparing numerical results with analytical ones for the case of two-dimensional steady-state heat conduction. We also checked the effects of varying At, length Ax, and width Ay of the elements on the temperature profile at the centre of the basin. The finite-element discretization of the computational domain is shown in Fig. 1. The input parameters used in the tests are listed in Table 1. Initially, the lithosphere is divided vertically into four layers and horizontally into three sections. The thickness of the lithosphere and crust are 130 km and 30 km, respectively. The crust is divided into three layers (and sediments where present), overlying the subcrustal lithosphere. Within each of these layers, thermal and petrophysical properties are constant. Consequently, the isotherms before extension will be horizontal, and the temperature distribution can be obtained by analytical methods (Liu, 1993). The vertical temperature distribution calculated analytically and the discrepancy between numerical and analytical solution are illustrated in Fig. 2. After 400 iterations with over-relaxation factor w= 1.75 (see Appendix A), the accuracy can reach 3 x lop3 K, which is more than satisfactory. Fig. 3 shows the relative error caused by varying At and Ax for the temperature profile at the centre of the basin 10 Ma after the cessation of extension, assumed to last 2 Ma with maximum /I = 2.0 at the centre. The relative error is defined as 0 = (T,- r,)/(T,,,T,,,,), where T,, Th are the temperatures for test a and test b, and T,,,,,, T,,, are the maximum and minimum temperatures in the computational domain. The results indicate that for At < 0.05 Ma, with the spatial discretization of Fig. 1 (where the larger elements have an area of 13 km*), the relative errors in temperature are less than 5 x 10P4. As At increases, the error increases. However, even for At = 0.2Ma, the relative error is still less than 3 x 10 3. However, as Ax increases from 5 to 10 km, the error increases drastically. It should be noted that as the elapsed time increases, the errors can also increase (Johnson, 1987). In this study, we use the standard spatial discretization shown in Fig. 1 and 0.1 Ma as time increment. MODEL TESTS Several test runs have been carried out in order to clarify the thermal evolution subsidence history of a basin as a function of b-factor and duration of extension,

and the

G. Liu and G. Ranalli

TEMPERATURE

DIFFERENCE (10” “C)

(“C)

0 200 400 600 80010001200 O',

\

10:

I

'

'

I

-

T

.\

201

10 : '* \ \ \ . ,

30 { 40 :

20 1 30 : 40 y I . . .

50 { 60 : 70 :

50 { 60 1 70 ;

*

80 :

.

80 :

90 :

90 :

100 y

100 1

110 1

110 1

. .

120 {

120 : 7

130 -I Fig. 2. Temperature

,

at the centre of hthosphere domam (analyttcal and numertcal solutions.

0

IO

20

30

40

50

60

70

solution)

80

90

and difference

100

110

between

120

analytical

130

DEPTH (km) Fig. 3. Relative error in temperature distributton at the centre of the basin for varions Ar and As after 10 Ma 01 of Fig. 4a. The solid line represents thermal cooling. The extension lasts 2 Ma with the initial o-factor distribution the relative errors between At equal to 0.01 and 0.05 Ma. the dotted line between 0.01 and 0.1 Ma, and the dashed line between 0.01 and 0.2Ma. All above tests use the same spatial discretization shown in Fig. lb. The dashdotted line shows the relative error between temperature distributions for Ax equal to 5 and 10 km and with the constant A/=0. I Ma.

9

A finite-element algorithm and thermal history of extensional basins

differences in thermal and subsidence history for a basin with water and sediment loading, and the lateral heat transfer for symmetric and asymmetric basins. In the first test, we consider a lithospheric domain 180 km wide and 130 km thick with the thermal and petrophysical properties listed in Table 1. The space formed by extension is assumed to be filled to sea level by sediments. Extension lasts 20Ma, followed by 150 Ma of cooling. The prescribed p-factor distribution is shown in Fig. 4a. (Note that this distribution refers to the initial position of nodal points. As a consequence of extension, the width of the extending region increases with time.) The subsidence and heat flow histories are shown in Fig. 5. Both subsidence and heat flow reach equilibrium within about 100 Ma. The equilibrium heat flow at the centre of the basin is lower than that in the initial state, because of the reduction in radiogenic heat production caused by the thinning of the

-90 -60 -70 -60 -50 -40 -30 -20 -10

0

IO

20

30

40

50

60

70

80

X-AXIS (km)

2.5 -~"."'.""""."""""""""',.",",""",',"',""'."'.""".

1.**cm @)

2.0 -

:

0.5 -

0.0 -~p,,,,,,,,,,,,,,m_I -60 -70 -60 -50 -40 -30 -20 -10

0

IO

20

30

40

50

60

70

60

X-AXIS (km) Fig. 4. (a) b-factor distribution

for Figs 5-7; (b) b-factor

distribution

for Figs g-10

90

IO

(a)

crust. Note also that the thermal etf‘ects of thinning arc more important than the effects of sediments. as shown by the syn-extensional incrcasc in heat How. Fig. 6 shows the mesh at the end of the extension (the initial nodal positions were shown in Fig. I). The convex lines between the two thick lines were initially horizontal. The mantle part ofthe lithosphere is thinned and lifted up. The temperature fields at the end ofextension and at 50 and 100 Ma after cessation of extension are illustrated in Fig. 7. In this test, WC may have over-estimated the heat flow and under-estimated the subsidence, in particular during and immediately after extension, as the isotherms at the end of extension are not yet horizontal near the vertical boundaries of the model. Thermal equilibrium is reached within 100 Ma from the end of extension. In order to evaluate the effect of the /i-factor. in particular on lateral heat conduction. a second test was conducted with the prescribed /i-factor distribution shown in Fig. 4b. In this test. the width of the undeformed sections is set to 60 km. and the initial width of the

II

A finite-element algorithm and thermal history of extensional basins

-

1000

L

4000

3 G 2ooo n

3000

5000

6000 ~,,I,,,,,,.,,,.,,,,,,,.,.,,.,,,..,....,.,,.,.,.,,....,.,,.,,,,,,.,~ -120 -100 -80 -60 -40 -20 0 X-AXIS

20

40

60

80

100

120

(km)

(b)

60 E %

80

100

-120 -100

-80

-60

-40

0

20

X-AXIS

(km)

-20

40

60

80

100

Fig. 6. (a) Deflection of basement, and (b) deformed grid at the end of extension for the case illustrated Curves between thick lines were initially horizontal: basin is filled with sediments.

120

in Fig. 5.

stretched section to 40 km. The extension lasts 10 Ma, followed by 150 Ma of cooling. Fig. Sa, b show the subsidence history and heat flow evolution for this model. The temperature field at the end of extension is shown in Fig. 9. The lateral temperature gradient is approximately zero near the vertical boundaries. It is expected that the higher the magnitude of extension, the larger the area that will be affected by lateral heat conduction. Subsidence and heat flow at the centre of the basin with sediment loading for dif’ferent

0 10 20 30 40

g

g

5o

60 70

B 60 90

100 110 120 130

-120-100

-80 -60

-40

-20

0

X-AXIS

20

40

60

80

100 120

(km)

(b)

-120~100 -60

-60

-40 -20

0

X-AXIS

-120 -100

-80

-60

-40

-20

0

20

40

60

60

100 120

40

60

60

100 120

(km)

20

X-AXIS

Fig. 7. Temperature

lields (a) at the end of exlenson. (b) 50 Ma. and (c) 100 Ma after the end of extension. the case illustrated m Fig. 5

tbl

A finite-element

algorithm

and thermal

history of extensional

basins

13

Fig. 8. (a) Subsidence history, and (b) surface heat flow for input parameters as in Table I, sedimentary filling, active extension lasting 10 Ma, and p-factor distribution given in Fig. 4b.

G. Liu and G. Ranalli 0 10 20 30 40 50 60 70 80 90 100 110 120 130

....._____________........~~~~

-..-............_............

I’.~‘,““,““,“‘.I”“I”~‘I”“I”“I’~

-80

-60

-40

-20

0

20

40

60

80

extension periods (10 and 20 Ma). and /&factor distribution as given in Fig. 4a, are compared in Fig. IO. A longer extensional phase predicts a slightly larger subsidence than a shorter phase. This may be caused by the combination of the longer cooling period and thermal blanketing. Thermal blanketing effects are noticeable in Fig. IO where high heat flow can be maintained more than IO Ma after the cessation of the extension for the shorter extension phase. A comparison of subsidence history for water loading and sediment loading with a maximum /1=2.5 is shown in Fig. I I. Although the sediment loading case predicts more than two times more subsidence than the water loading case, the shoulder uplift is similar. In order to isolate the thermal effects of sediment loading. we also consider a case of ‘nonthermal sediments’. i.e. a material with the same density of sediments (see Table I). but with a tictitiously high thermal conductivity of SO0 W m ’ K ‘. The differences between this case and the sediment loading case are due to thermal cffccts only. Fig. I2 shows the predicted subsidence and heat flow at the basin centre for water, sediment, and non-thermal sediment loading. The thermal blanketing effects are clearly visible in the subsidence curves. The water loading model predicts a higher heat flow than the sediment loading model. The maximum difference can be as much as 20 mW m ’ for a maximum p-factor of 2.5. This result is mainly due to two factors: the thinning of the conduction length and the thermal blanketing effect. The lithosphere is thinner in the case of water loading because the space originated by extension at the top is not occupied by sediments. (The higher heat flow in the case of non-thermal sediments is an artifact due to the larger initial subsidence.)

A finite-element algorithm and thermal history of extensional basins

0

10

20

30

40

50

60

70

80

1.5

90 100 110 120 130 140 150 160 170

TIME (Ma) Fig. 10. Comparison

of (a) subsidence and (b) heat flow at the centre of the basin for /?-factor distribution in Fig. 4(b) and different duration of the active extension phase.

as given

All previous tests are concerned with a symmetrical distribution of the p-factor. We have also run a test with asymmetrical B-factor distribution as shown in Fig. 13. The extension lasts 10 Ma, followed by 50 Ma of thermal relaxation. The basin looks like a half graben. The maximum uplift at the left shoulder is about 400 m. If erosion is instantaneous (Liu et al., 1992) the maximum thickness of the material removed from the shoulder can be as much as 1500m. The heat flow evolution for this asymmetrical p-factor distribution is illustrated in Fig. 14 as a function of position and time. CONCLUSION The two-dimensional finite-element algorithm presented in this paper is simple, flexible, and takes into account the different thermophysical properties of different layers of the lithosphere. It can be used to predict subsidence and heat flow patterns in kinematic models of lithosphere extension with prescribed b-factor distribution and local Airy isotasy. Several

16

G. Liu and G. Ranalli

-125

-100

-75

-50

-25

0

X-AXIS

25

50

75

100

125

(km)

-500 0 g

500

tf:

1000

g n 5;

2000

$

2500

1500

3000

X-AXIS (km)

0

3

1000

G--2ooo

g

3000

g

4000

3

5000

,3

6000 7000 8000 -125

-100

-75

-50

-25

0

25

X-AXIS (km)

50

75

100

125

A finite-element algorithm and thermal history of extensional basins

(4

: 1

-mRMALsEDIMEN-f

_

“‘~~~~.........._._._._.,.,,_,,,,

17

WATER

-.

NONTHERMAL SEDIMENT

THERMAL. SEDIMENT WATER NONTHERMU SEDIMENT

0

10

20

30

40

50

60

TIME (Ma) Fig. 12. Comparison

of (a) subsidence and (b) heat flow at the basin centre for the case illustrated water. sediment, and non-thermal sediment loading.

in Fig.

11for

C;. Liu and G. Ranalli

IX 2.5

-80

-60

-40

-20

0

20

X-AXIS

40

60

80

100

(km)

-

2.5Ma 5.0 7.5 10.0

--

15.0

-

-

- 20.0

- 25.0

__ ” t:

--’ -.-

6000 7000

-80

-60

-40

-20

0

20

X-AXIS

40

60

80

30.0 35.0 40.0 45.0 50.0 55.0 - 60.0

100

(km)

SL

g ii

4000 1

$ E

6000 5000 -'

3ooo I

\

\:t

1

.---_

.--_

-

-

-.-

i

-.---__

_

20

I 7 '7 30

TIME

(Ma)

_

_,._. -

F---..-_-

,TT--, IO

_ ---

\

8000 7000 ! --,.,,. 0

-

_._.-

w-w 1 "

40

'.

1 50

”

“1

60

-

A’

-

B'

-.-

C’

-

D’

A finite-element algorithm and thermal history of extensional basins

Fig. 14. Heat flow evolution for the case illustrated

19

in Fig. 13

test runs have shown that time increments equal to, or less than, 0.1 Ma are appropriate for triangle element size of the order of 10 km’. Numerical results are sensitive to the size of the elements. Simple applications have been run to investigate the effects of the p-factor distribution (including asymmetric stretching) and the duration of the extension phase on the subsidence and thermal history. A comparison has also been made of the relative effects of water, sediments, and ‘non-thermal sediments’, which allows thermal blanketing effects to be isolated. Although the examples presented in this paper are very simplified representations of reality, the algorithm is applicable to more complex spatial distribution of thermophysical properties, and consequently can be used to mode1 geological situations in detail. A~knnlvlrdgements-The work reported in this paper has been supported by an NSERC (Natural Sciences and Engineering Research Council of Canada) grant to G.R. Wolf Jacoby. Mane1 Fernandez, and an anonymous reviewer are thanked for useful criticism.

APPENDIX

A: ALGORITHM

A weighted residual method of finite-element formulation was adopted. A domain is discretized into triangular elements. If the number of elements is E, the total residual can be expressed as

(Al) and the residual

for each element

is

G. Liu and G. Ranalli

20

(A21

The vector of shape functions

is N,(s. .v) N&.x, .I’)

A” and II are the area and total nodal 17= 3. After some algebraic manipulation, d(T; IV” = [M] dt The mass and stiffness matrices

number of element eqn (A2) becomes

(1. For triangular

elements,

+ [k]( T ;,

are (A5)

(Ah)

where (N,)

is the transpose

of 1N, ]. and

(A7)

(As) The load vectors

i,fj,Jand

if;; in eqn (A4) can be written

as (A9)

n

(.f
(AlO)

where

(Al 1)

A finite-element algorithm and thermal history of extensional basins

Fig. 15. Transformation

from real domain

(x,y) to reference

domain

21

(5, q)

In terms of reference coordinates [the mapping from real coordinate (x, y) to reference coordinate (4, q) is illustrated in Fig. 151, eqns (A3t_(All) can be rewritten as W<> =

6412)

[l-t-%Lfd

p-110,

(A13)

~=(-_1,01]

(A14)

is ’ ss I

[ml

=

I-:

(N,;> det [Jl dv dt

6415)

[41T[.Ar[~l[AP.ddet [Jl dvldt

6416)

{NJ

wp

0

0

1-r

[k]=

0 0

(A17)

(‘418) Yl -Y2

x2-x, det[J]

= 2A’

1

(A19) (A20)

G. Liu and G. Ranalli

(A21)

(A22) QIll

+

ii;: I

The local mass matrix is asymmetrical,

definite and positive.

h-,, = /<,(.I.~ .I';)( .I': -.L';) i K,(\i

kj, = I<,(,l’~~~‘:)(~l’; /\I 1 z x,,

are of the form

.\-J)(.\-? -.VJ)

~~.l.,)+/i,(.\-;~-_\-~)(.\-,-.\-;)

K,( -1‘:mm_l'j)(_l'~ -.I';)+ li,(.\3

-= li,(,l‘:m l',)(,l':

Its elements

.\J)(.\-2-_Yf

)

I',)+,<,(.\,

Y\-J)(.\-,-_Yq)

I,-)$-/<,( \s

.\-J.\-~-_\-,)

x,: = /<,(.I.:

I,,)( l’,

I<;: = K,(.l’,

.I,:)( .I,, -.l’,)+K,(.\,

(A23)

-_v,)(_Y~-_Y,).

It is assumed that the material properties within each element are constant. but may not can be different from the be isotropic. For instance. the vertical thermal conductivity horizontal one. After integration, the clement mass matrix [eqn (A IS)] becomes

(A231

and the element

load vector [eqn (A IO)] become5

(A25)

The only non-zero where the temperature Note that there arc given row as illustrated the iteration equation

nodes, terms of the load vector 1/,,I are those related to boundary is defined. Therefore, the load vector does not need to be evaluated. at most seven non-;rero entries in the matrices [M] and [Kj in any in Fig. 16. Therefore. for row j4 of the linear system of equations. can be written as (A26)

where F and 7‘are the temperatures

from the last and current

h,,l = c [(/II,,P A.oT,“-t F,] ’ 1, (I,

=

IH,~I..A~+ k,,/..

iterations.

respectively,

and (A27 (A28)

A finite-element

algorithm

and thermal

history of extensional

Fig. 16. Connection between nodes and elements. The value ofj, is determined

basins

by its six neighboring

23

nodes j, to,j,.

f--dFig. 17. A sketch showing conservation of mass during lithospheric extension with pure shear. The left part of the figure is a trapezium before extension and the right part after the extension. Each trapezium is divided into two triangular elementsq, = h,J/I,, q2= h,//L; b, and /I2 are extension factors at the left and right side of the trapezium.

where m,4r and k,?r are the entries of the global mass matrix and global stiffness matrix at row ,j4. The superscripts n and nf 1 indicate the times t and t+At. The over-relaxation factor w ranges in this case from 1 to 2 (Press et al., 1992). For conservation of mass, the width of each trapezium in the computational domain is calculated by:

(~29)

The symbols

in eqn (A29) are explained

in Fig. 17

24

G. Liu and G. Ranalli REFERENCES

Ayala, C. and Pous, J.. Tornk. M. ( 1996) 3D gravity model of lithospheric structure in Valencia Trough (abstract). A/in. Geoph~,s. 14, 103. Barton. P. and Wood, R. (1984) Tectonic evolution of the North Sea basin: Crustal stretching and subsidence. Geopl?,,.v. J. R. Astr. Sec. 79, 987-1022. Brace, W. F. and Kohlstedt, D. L. (1980) Limits oflithospheric stress imposed by laboratory experiments. J. Geoph?,.s. Res. 85, 6248-6252. Braun, J. and Beaumont, C. (1987) Styles of continental extension from dynamic models of lithosphere extension. Can. Sot. Petrol. Geol. Memoir 12, 41--258. Boutilier, R. R. and Keen. C. E. (I 994) Geodynamic models of fault-controlled extension. Tectonics 13, 439454. Burov, E.B. and Diament, M. (1995) The effective thickness (TC) of the continental lithosphere: What does it really mean? J. Geophj:r. Re.v. 100, 3905-3927. Cabal, J. and Fernandez, M. (1995) Heat flow and regional uplift at the north-eastern border of the Ebro Basin. NE Spain. Gcoph.~~.v.J. Int. 121, 393403. ChCry, J., Lucazeau. F., Daignikres. M. and Vilotte. J. P. (1990) The deformation of continental crust in extensional /ones: A numerical approach. In The Potentiul qj’ Deep Seismic Prcfiling ,fbr H~~drocarhon E_~plorution (Pint P. and Bois C. eds), pp. 157-166. Technip. Paris. Cloetingh, S. (1995) Ori~gin und Ewlution of Srdimentur~~ Basins, Proceedings of the VIII Summer School of Earth and Planetary Sciences, pp. 333 394. Siena. Cloetingh, S. and Burov, E. B. (1996) Thermomechanical structure of European continental lithosphere: Constraints from rheological profile and EET estimates. Geop$l:F. J. Znt. 124, 695-m723. Cloetingh, S., Eldholm, 0.. Larsen, B. T.. Gabrielsen. R. H. and Sassi. W. (1994) Dynamics of extensional basin formation and inversion: Introduction. Tectonophysics 240, 1-9. Cloetingh, S.. Sassi, W., Horvath. F. and Puigdefabregas, C. (1993) Basin analysis and dynamics of sedimentary basin evolution. Sediment. Geol. 86, 1~201. De Bremaecker, J.-C. (1983) Temperature and hydrocarbon maturation in extensional basins: A finite element model. Am. ASSOC. Petrol. Geol. Bull. 67, 1410-1414. Giilkes, M.. Cioetingh, S. and Fuchs. K. (1994) Finite-element modelling of pull-apart basin formation. Tec,tonoph!,,sic,.v 215, 167 185. Govers, R. and Wortel. M. J. R. (1993) Initiation of asymmetrical extension in continental lithosphere. Tectonoph!,.sic.v 223, 75 96. Hermanrud, C. (1992) Basin modelling techniyues an overview. In Basin Modelling Adt,ances and Applicutions (Do& A. G., Augustson J. H.. Hermanrud C., Steward D. J. and Sylta M. eds). Vol. 3, pp. I 34. Nor. Petrol. Sot. Spec. Publ. Johnson, C. ( 1987) Numericul Solution of’Purtiu1 D@rentiul Equations b?xthe Finite Element Method, p. 278. Cambridge University Press. Cambridge. Keen, C. E. and Boutilier, R. R. (1990) Geodynamic modeling of rifted basins: The synrift evolution of a simple half graben. In The Potentiul qf Deep Seismic Prqfiling,fiw Hvdrocurhon Exploration (Pint B. and Bois C. eds), pp. 23-34. Technip, Paris. Ko&. H. and Cloetingh, S. (1992) Lithospheric necking and regional isostasy at extensional basins 1: Subsidence and gravity modelling with an application to Gulf of Lions margin (SE France). J. GeophJ-s. Res. 97, 17553 17571. Lippard. S. J. and Liu. G. (1992) Tectonic modelling of northern North Sea using program RIFT. In Structural and Tectonic Modelling und Its Application to Petroleum Geolog>l (Larsen, R. M., Brekke, H.. Larsen, B. T. and Talleraas, E. eds), Vol. I, pp. 43~~54. Nor. Petrol. Sot. Spec. Pub. Liu. G. (1993) An algorithm and sensitivity study of one-dimensional tectonic modelling using a finite-element method. C’omput. Geosci. 19, 1101~~11I3.

A finite-element algorithm and thermal history of extensional basins

25

Liu, G. (1994) Tectonic modelling in the Bjarntirya West Basin of the west Barent Sea. Geophys. Prospect. 42,217-302. Liu, G., Lippard, S. J., Fanavoll, S., Sylta, M., Vassmyr, S. and DorC, A. G. (1992) Quantitative geodynamic modelling of Barents Sea uplift and erosion. Nor. J. Geol. 72, 313-316. Lynch, H. D. and Morgan, P. (1987) The tensile strength of the lithosphere and location of extension. In Contintental Extensional Tectonics (Coward M. P., Dewey J. F. and Hancock P. L. eds), Vol. 28, pp. 53-65. Geol. Sot. Lond. Spec. Pub. McKenzie, D. P. (1978) Some remarks on development of sedimentary basins. Earth Planet. Sci. Lett. 40,25-32. Odling, N. E. (1992) RIFT, a model of sedimentary basin evolution by finite rate nonuniform, pure shear extension of the lithosphere. In Structural and Tectonic Modelling and Its Application to Petroleum Geology (Larsen, R. M., Brekke, H., Larsen, B. T. and Talleraas, E. eds), Vol. 1, pp. 457468. Nor. Petrol. Sot. Spec. Pub. Press, W. H., Teukolsky, S. A., Vetterling, W. T. and Flannery, B. P. (1992) Numerical Recipes in FORTRAN (2nd edn). Cambridge University Press, Cambridge. Ranalli, G. (1995) The Rheology of the Earth (2nd edn). Chapman and Hall, London. Ranalli, G. (1997) Rheology and deep tectonics. Ann. Geojis., in press. Ricke, M. and Mechie, J. (1989) Finite element modelling of continental rift structure (Rhinegraben) with a large-deformation algorithm. Tectonophysics 165,8 l-9 1. Royden, L. and Keen, C. E. (1980) Extension processes and thermal evolution of the continental margin off eastern Canada determined from subsidence curves. Earth Planet. Sci. Lett. 51, 343-361. Sclater, G. T., Royden, L., Horvath, F., Burchfiel, B. C., Semken, S. and Stegena, L. (1980) The formation of intra-Carpathian basins as determined from subsidence data. Earth Planet. Sci. Lett. 51, 139-162. Skogseid, J., Petersen, T. and Larsen, V. B. (1992) Vsring Basin: Subsidence and tectonic evolution, In Structural and Tectonic Modelling and Its Application to Petroleum Geology (Larsen, R. M., Brekke, H., Larsen, B. T. and Talleraas, E. eds), Vol. 1, pp. 55-82, Nor. Petrol. Sot. Spec. Pub. Stephenson, R. A., Embry, A. F., Nakiboglu, S. M. and Hastaoglu, M. A. (1987) Riftinitiated Permian to Early Cretaceous subsidence of Sverdrup Basin. In Sedimentar), Basin Formation Mechanisms (Beamount C. and Tankard A. J. eds) 12,213-23 1. Thorne, J. A. and Watts, A. B. (1989) Quantitative analysis of North Sea subsidence. Am. Assoc. Petrol. Geol. Bull. 73, 88-l 16. Van der Beek, P. and Cloetingh, S. (1992) Lithospheric flexure and the tectonic evolution of the Beltic Cordilleras (SE Spain). Tectonophysics 203, 325-344. White, F. M. (1991) Viscous Fluid Flow (2nd edn). McGraw-Hill, New York.