• Keine Ergebnisse gefunden

Numerical modelling of melting processes and induced diapirism in the lower crust

N/A
N/A
Protected

Academic year: 2022

Aktie "Numerical modelling of melting processes and induced diapirism in the lower crust"

Copied!
12
0
0

Wird geladen.... (Jetzt Volltext ansehen)

Volltext

(1)

Geopkys. J . Int. (1995) 123,59-70

Numerical modelling of melting processes and induced diapirism in the lower crust

D. Bittner' and H. Schmeling2

Jean-Paul-Str. 21, D-95444 Bayreutk, Germany

'lnstitut fur Meteorologie und Geophysik, Feldbergstr. 47, D-60323 Frankfurt am Main, Germany

Accepted 1995 March 27. Received 1995 January 31; in original form 1994 June 27

S U M M A R Y

The processes of partial melting and magmatic diapirism within the lower crust are evaluated using a numerical underplating model. Fully molten basalt ( T = 1200°C) is emplaced at the Moho beneath a solid granite ( T = 750OC) in order that a melt front grows into the granite. If diapirism does not occur, this melt front in the granite reaches a minimal depth in the crust before (like in the molten basalt) crystallization takes place. The density contrast between the partially molten granite layer and the overlying solid granite can lead to a Rayleigh-Taylor instability (RTI) which results in diapiric rise of the partially molten granite. Assuming a binary eutectic system for both the granite and the underplating basalt and a temperature- and stress-dependent rheology for the granite, we numerically solve the governing equations and find (a) that diapirism occurs only within a certain but possibly realistic range of parameters, and (b) that if diapirs occur, they do not rise to levels shallower than 15 or perhaps 12km. The growth rate depends on the degree of melting and the thickness of the partially molten layer, as well as the viscosity of the solid and the partially molten granite. From a comparison of the growth rate with the velocity of a Stefan front it is possible to predict whether a melt front will become unstable and result in diapiric ascent or whether a partially molten layer is created, which remains at depth. We carry out such a comparison using our thermodynamically and thermomechanically consistent model of melting and diapirism.

Key words: diapirism, granites, melting processes, phase-driven convection.

INTRODUCTION

While the occurrence of granitic intrusions at both shallow (Pitcher & Cobbing 1985) and deep (Ocala & Myer 1972;

Schwartz et af. 1984) crustal levels has been described for many regions, the physical processes of their formation, rise and emplacement are not well understood. In principle, large amounts of felsic melts in the lower crust can be produced by:

(1) a thickening of the crust in thrust regions (Buck &

Toksoz 1983) associated with heating due to thermal relaxation;

(2) a regional increase of heat flow from the mantle (De Yoreo, Lux & Guidotti 19891, due to upwelling of asthenosphere (Bird 1978);

(3) basaltic underplating (Huppert & Sparks 1988);

(4) internal heating due to anomalous distribution of radio- active elements at mid-crustal depths (Lathop, Blum &

Chamberlain 1994).

The amount of granitic melt due to basaltic underplating depends mainly on the volume of the underplating body. From

a simple 1-D analysis one can derive that basaltic magma intruding at Moho depths has the potential to partially melt the overlying felsic crust up to a thickness comparable to that of the basaltic layer (see also Bergantz 1989). In the case of lithospheric root delamination (Bird 1979) the asthenosphere might rise up to levels as shallow as the Moho, inducing extensive melting of the lower crust (Lorenz & Nicholls 1984;

Matte & Nicolls 1986).

Different mechanisms of melt removal from the source region and further ascent of felsic magma have been proposed.

Viscosities of felsic magmas are usually too high to allow for significant distances of rise by segregation (McKenzie 1985;

Wickham 1986). However, segregation through vein-like net- works might help to accumulate the melt within the upper part of an initially partially molten granitic body, as observed in the, Ivrea zone (Oncken 1993, private communication).

Models of diapiric ascent of fluid spheres (Stokes' flow) into a crustal layer with temperature-dependent viscosity (Mahon

& Harrison 1988) or a non-Newtonian rheology (Weinberg &

Podladchikov 1994) show that the shallowest possible emplace- ment depth is about 15 km. Somewhat shallower emplacement

0 1995 RAS 59

(2)

60 D. Bittner and H. Schmeling

depths were found by Reuther & Neugebauer (1987) and Kukowski & Neugebauer (1990) who used a rather weak temperature-dependent viscosity of the crust. Stoping may be a viable mechanism, but it can explain only small distances of magmatic rise (Oxburgh & McRae 1984). Dyke transport has been shown to be a very efficient ascent mechanism by Clemens & Mawer (1992), while Rubin (1993) argues that dyke propagation might be inhibited due to high viscosity of felsic magma.

Most of the above-mentioned approaches do not consider melting/solidification processes during the rise. They also assume simplified rheologies of the ductile granitic crust. In the present study we focus on a thermodynamically and thermomechanically consistent model of granitic plutonism due to basaltic underplating. Of particular importance is the rheology of the solid and partially molten granite. Non- Newtonian rheology laws based on laboratory experiments (Kirby & Kronenberg 1987; Shaocheng & Pinglao 1993) are used.

If a granitic unmolten crust is underplated by a molten basaltic layer, a melting front migrates into the crust. In the absence of diapirism this front reaches a minimum crustal depth before solidification of both the granitic and basaltic layers occurs. During this process the density of the partially molten granitic layer is less than that of the overlying crust, giving rise to the possibility of a Rayleigh-Taylor instability (RTI) with subsequent diapiric rise of the granitic magma body. The growth rate of the RTI will depend on the degree of melt, the thickness of the partially molten layer and the viscosity of both the partially molten layer and the solid overburden. These parameters may be estimated from 1-D Stefan-problem-like analyses of the underplating problem (Bergantz 1989) and be used to estimate the growth rate. This growth rate can be compared with the velocity of the penetrat- ing or retreating Stefan front. From such comparisons it should be possible to predict whether a melting front penetrating the base of the crust will become unstable and rise as diapirs, or whether it will stay near the basaltic layer. We carry out such a comparison using the full solutions of the thermodynamical- thermomechanical problem of melting and diapirism.

THEORETICAL F O R M U L A T I O N

In the formulation of our binary, two-phase, multirock model we distinguish between N different rock types, each of which may have different physical properties and may undergo melting according to different binary phase diagrams. Two different rock types will be used in our model. The basic underplating rock will be referred to as basalt while the acidic one above may be called granite. The governing equations describe the conservation of mass, chemical components, momentum and energy:

ap/at

+ v

(pu) =

o

( 1 )

x = 1

( S i + M ’ ) = l ,

i = 1 ... N

o = - V P + - q

axj a ri -+- axj axi

& j ) + p g , (3) aT/at

+

( U V ) T = “ ~ ~ - ( p ~ ~ , ~ / p C , ) ( a ~ / a t

+

( ~ v ) M ) . (4) The definitions of variables are given in Table 1

Table 1. Variables.

a := thermal expansivity T :=temperature p :=density v :=velocity q :=viscosity

Tsol :=solidus temperature Tliq :=liquidus temperature g := gravitational acceleration k :=thermal diffusivity

E :=energy P :=pressure X :=mass in percent M :=liquid mass in percent S :=solid mass in percent

I :=1 ... N number of rock types

L :=latent heat

ZII := 2nd invariant of deviatoric stress tensor

PO :=min(psolid.,Pliqud

In the momentum equation inertial terms are neglected (high Prandtl number approximation). The Boussinesq approximation is used, i.e. density is kept constant except for the buoyancy term.

The equation of state is given as

( 5 )

=

C

(pi.solidSi + pi,liquidh,ii) - p o c t ~ ,

i = I ... N

The viscosity for solid and partially molten rock with up to 10 per cent melt is assumed to be stress- and temperature- dependent:

q = Az:,-”exp(E/RT). (6)

At higher melt fractions ( M > 10 per cent) an empirical rheological law depending on the melt fraction is used (Pinkerton & Stevenson 1992):

q = ~ ~ e x p ( 2 . 5

+

[(I - M ) / M ] O . ~ ~ ) ( 1 - M ) ) (7) The amount of melt M = M ( T ,

Kol,

&,,L) is a function of temperature, solidus and liquidus temperature and the latent heat. It is calculated for a fully binary or a eutectic system using analytical equations for the solidus and liquidus curves (Paufler 1981) and the lever principle.

N U M E R I C A L F O R M U L A T I O N

While the heat capacity and the thermal conductivity k are different for granite and basalt, they are assumed to be temperature-independent. In the model we assume no multiphase flow or segregation.

A finite-difference (FD) code with uniform grid size is used.

For solid rock an explicit FD formulation of the heat eq. (4) is used (forward in time, central in space). If melt is present, an explicit formulation of the time derivatives is no longer possible because of the implicit form of M ( T ) . An iteration algorithm has to be used: first the heat equation is solved without melting terms, to get a lower or upper bound for the temperature; then eq. (4) and M = M ( T ,

KO,, xi¶,

L) from the phase relation are solved iteratively to obtain the temperature at the new time step.

The conservation of momentum (eq. 3) is solved by a stream- function approach which fulfils also the conservation of mass (eq. 1). The mass transport of different rock types is simulated by markers. The movement of these markers is calculated by a Runge-Kutta algorithm of fourth order.

A problem in the numerical calculation of this underplating model arises from the large viscosity contrast between liquid basalt (q < lo5 Pas) and solid or partially molten granite ( q > 1014 Pas). Accordingly, the Rayleigh numbers in the basaltic and granitic regions might differ significantly, and 0 1995 RAS, GJI 123, 59-70

typo: left paranthesis missing

(3)

Modelling of melting processes and diapirism 61

viscosity of granite

23

n (II

CI 0

F

21

n.

m 19

>r 17

.- 3

._

+.’

15

> 13 11 9

- -

Haplogranite real power-law 0.1 MPa -power-law 1 MPa -power-law 10 MPa -melt in model

Q____*

520 640 760 880 temperature [“C]

Figure 1. Viscosity of granite as used in the models. In the first region where T < T,,,,, a power-law rheology is used (eq. 6 ) . Here the viscosity calculated for stresses of 0.1 MPa, 1 MPa and 10MPa are shown. The stresses created in the model are less than 10MPa. In the second region where T z T,,,,, the viscosity is calculated from eq. (7) for melt concentrations ( M ) greater than 10 per cent. Since within the present resolution, the Rayleigh numbers larger than 105-106 lead to numerical instabilities, we artificially increase the viscosity of the highly molten regions and do not use the data from Dingwell, Knoche & Webb (1992) shown as a dashed curve.

values of Ra > 10” are easily possible within a molten basaltic layer of several kilometres thickness. Since within the present resolution Ra numbers larger than 10s-106 lead to numerical instabilities, we artificially increase the viscosity of the highly molten regions.

The rheologies chosen for granite are shown in Fig. 1.

Below melt concentrations of 10 per cent the constants A = 6 x 10’ MPa” s, n = 3.2 and E = 144 kJ mol-’ are used for granite (based on Kirby & Kronenberg 1987). Effective vis- cosities are shown for flow stresses of 0.1-10MPa. As can be seen in the diagram, at stresses consistent with diapiric dynam-

as low as 1OI6 Pa s. At melt concentrations above 10 per cent, eq. (7) is used with qo = 5 x 1014Pas for granite (Fig. 1, curve

‘model’) and qo = loi3 P a s for basalt (not shown). The values of qo have been chosen arbitrarily in order to ensure numerical stability in the liquid regions. This high viscosity implies that the time-scale in the liquid rock is too long (i.e. the number of overturns per time will be smaller than that in natural systems), but we can simulate the mechanical behaviour of the ductile granite (up to 10 per cent melt) and in particular the rise of the diapirs with some confidence. For comparison, an exper- imentally derived curve for molten granite ( M > 30 per cent) ics (- 10 MPa)the viscosity in the granite may d r i p to ialues is shown (‘haplogranite’).

-

M -T relation Ra pa kivi-g ran ite

+P = 5 kbar --S-..P = 2 kbar

0 20 40 60 80 1 00

melt

proportion M

(wt. %)

Figure 2. M-T relation of Rapakivi granite ( P = 2 kbar and P = 5 kbar with 2.8 per cent H,O) after Nekvasil (1991). The effects of the different components of the system can be seen as curve segments with different slopes.

0 1995 RAS, G J I 123, 59-70

(4)

62 D. Bittner and H . Schmeling

density basalt solid latent heat granite

PHASE D I A G R A M S A N D D E N S I T Y VARIATIONS U S E D I N THE MODEL

Lower crustal acidic or basic rocks may exhibit large variations in melting behaviour, but can be generally characterized by eutectic systems. The actual solidus and liquidus temperatures and the eutectic point depend on pressure, water content, C 0 2 content, and composition. Since we do not intend to simulate particular petrological systems, the actual choice of melting parameters is not important, as long as the general distinction between acidic (granite) and basic (basalt) rocks is made.

Melting may control the dynamic and thermal processes in the lower crust by influencing density and by the release of latent heat. Thus it is essential to know the melt fraction-tem- perature ( M - T ) relation. This relation can be derived from experiments or phase diagrams of multicomponent systems.

Typical granite systems can be characterized by eutectic melting

2900 k@m3 density basalt liquid 2800 k@m3 300 kJkg latent heat basalt 380 kJ/kg

phase diagram

thermal conductivity granite

0" 1100

2 1000

0 Y

3

Q) c,

900 800 700

c, Q)

2.5 W/(m*K) thermal conductivity basalt 3,5 W/(m*K)

0 0,2 0,4 0,6 0,8 1 X (composition)

M-T relation

a 1 100

Y

2 1000

c,

a

900

Q)

800

*

Q)

700

0 095 1

melt concentration M

Figure 3. Three different phase diagrams (Type 1, solid line; Type 2, dashed line; Type 3, dotted line) are used for granite in the models.

The start composition of the granite at X = 0.8 is shown in Fig. 3(a).

A melt concentration-temperature (M-T) relationship (Fig. 3b) can be calculated from these phase diagrams (Fig. 3a). Going from Type 1 to Type 3, the starting composition becomes closer to the eutectic point. This results in a higher maximal melt concentration at the eutectic temperature.

(i.e. increasing the melt fraction at constant temperature) followed by an increase of melt with increasing temperature.

As an example, Fig. 2 shows the melting behaviour of Rapakivi granite (Nekvasil 1991). The effect of the different components of the system can be seen as curve segments with different slopes.

To model the essential features of experimental eutectic melting curves (i.e. horizontal branches of different lengths followed by an increasing M-T curve) it is sufficient to use a eutectic binary system. We assume that the kinks in the M-T curves as seen in real multicomponent systems (Fig. 2) are of minor importance to the dynamics of the model compared to the effect of the eutectic part of the M-T relation.

The three phase diagrams used in the model to represent the granite and the resulting M-T relationships are shown in Fig. 3. There are two different types of M-T dependences seen in the eutectic system. Near the eutectic temperature there is an almost horizontal segment which is assumed to increase by 1 "C along its entire length. This increase is sufficient to yield an unequivocal relation between M and T , required by the numerical formulation, Above the eutectic temperature the curve increases in a convex manner. Variations in the length of the horizontal part of the curve control the melting behav- iour and therefore the dynamics of the model. This variation in the M-T curve can be due either to different eutectic systems as shown in the diagram or to different starting compositions within the same eutectic system. Fixing the eutectic point and varying the starting composition would lead to the same M-T curves as those shown in Fig. 3(b). Only the eutectic point in the granite phase diagram is varied.

The basalt is approximated using a eutectic phase diagram of diopside-anorthite at 10 kbar water pressure. Using a start- ing composition of X = 0.6 anorthite, the liquidus and solidus temperatures are 1100 "C and 950 "C, respectively.

An additional parameter which has a large influence on the dynamics of the model is the density difference between the solid and molten granite. The typical density contrast lies between 50kgm-3 and 200kgm-3. The values used in the model are listed in Table 2.

B O U N D A R Y A N D I N I T I A L C O N D I T I O N S As initial conditions (Fig. 4), the lowermost 5 km of the model are occupied by a totally molten basalt with a temperature of 1200 "C (100 "C above liquidus temperature). In the overlying solid granite a linear temperature profile with a gradient of 28 "C km-' is assumed, representing a heat flow of 70 mW m ~ '.

At 27km depth the solid granite is 10°C below the eutectic temperature. The thermal boundary conditions are constant Table 2. Parameters used in the calculations (Johannes & Holtz 1990;

KTB Feldlabor 1994; Pribnow & Umsonst 1993).

density granite solid( 2450-2600 ks/rn31 density granite liquid( 2400 k@m3

1

manmum heat flow1 70 mW/rn2 I minimum heat flow1 30 mW/W

I

heat capacity solid( loo0 J/(kg'K)

I

heat capacity liquid/ 1500 Jl(k9.K)

I

0 1995 RAS, GJI 123, 59-70

(5)

Modelling of melting processes and diapirism 63

T=const=336"C Temperature at t=O

12 km 17km 22km 27km 32km

T

7OmW/W

u)mw/w

336°C 756°C 1200°C

q (heatflow)

Figure 4. Initial and thermal boundary conditions for the underplating model with solid granite and totally molten basalt.

temperature ( T = 336°C) at the top and a laterally varying heat flow between 70 mW m-' and 30 mW m-' at the bottom.

Free slip at the top and bottom is chosen as the mechanical boundary condition. Due to the high viscosity at the top, this condition is equivalent to the no-slip condition.

RESULTS

Several calculations have been made, varying the density of the solid granite (Table 2) and the eutectic composition of the granite in the phase diagram (Fig. 3). First the typical time evolution of the underplating model is described. For this a model was chosen in which diapirs grow over a period of c. 0.7Ma. The melt fraction M , temperature T , stream function @ and the chemical components X are shown in Fig. 5. An active marker-line (the basalt-granite boundary) and a passive marker-line within the granite have been used to represent X . The Type 1 phase diagram and a density variation of 100kgm-3 between the solid and molten granite is used in this model. During the first phase (here 0-0.2 Ma) the cooling of the basalt leads to strong convection while the growth of the melt front into the granite is purely conductive.

The first crystals begin to grow in the basalt. During the second time phase (here from 0.2 Ma to 0.7 Ma; Fig. 5a) the thickness of the partially molten granite is sufficiently high to allow for convection within the granite. This results in pertur- bations of the interface betweens solid and molten granite which can potentially form a RTI. The growth of diapirs is, however, hindered by the dominating melting process. The basalt begins to solidify during this phase. Diapirs begin to rise out of the partially molten granite in the third phase (here between 0.7Ma and 1.5Ma; Figs Sb, c and d). The diapirs reach their maximum height after 1.5Ma. The dominating process during phases 3 and 4 is the convection within the granite driven by the solid-liquid phase transition. As a characteristic feature of this phase-transition-driven convec- tion, there exists a spatially stable partially molten region (Figs 5c, d and e) that remains essentially undisturbed during a large number of overturns. As a result, granite particles

involved in this flow undergo multiple cycles of melting and crystallization. The crystallization of the basalt begins at the contact to the granite and a crystallization front grows from above into the basalt. Due to diapirism and the ensuing phase- transition-driven convection in the granite, the basalt solidifies faster than for the purely conductive case as the heat transport to the surface is more effective. Subsequently in the fourth phase (Figs 5e and f), the cooling process begins also in the granite. Similar time phases exist for other models, with the time-scales of heating, ascent of diapirs (when existent) and cooling being determined by convection of differing strengths.

In a first series of models the density contrast between solid and liquid granite is varied (Ap = 50-200 kg m-3). Figs 6(a)-(c) show three different models at the same instance of time, which corresponds to time phase 4 as discussed above.

It can be seen that rising diapirs do not develop for all cases. The occurrence of a RTI in our models depend on two time-scales, namely the characteristic growth time of an insta- bility and the time-scale associated with the existence of a granitic melt layer. At low Ap the first time-scale is too large compared with the second one, prohibiting the development of diapirs.

These two time-scales do not only decide whether or not diapirs ascend, they also control the time at which a diapir starts to rise. The different density variations lead to different growth rates of the RTI which results in the rise of diapirs at different times. The time at which a diapir rises given a high value of Ap is much earlier than for a small density difference.

For large values of Ap, therefore, the diapirs rise out of a thin partially molten layer of granite (Fig. 6c). This results in smaller volumes of melt, which rise by diapiric flow. Smaller melt volumes cool and solidify more effectively and no over- hangs are observed compared with cases where diapirs start rising at later times and with larger volumes (Fig. 6b). Further observations are pulsating flows for high Ap and near steady- state flows for a smaller density contrast. This behaviour can be compared to thermal convection, which is time-dependent at high Ra and steady state at low Ra.

Variations in the phase diagram as indicated in Fig. 3(a) 0 1995 RAS, G J I 123, 59-70

(6)

64 D. Bittner and H . Schmeling

amount of melt M temperature T stream function Q, chemical field X

12

(4

16

20 -20

24 - 24 - 24

Time 190 Ka

28 - 28 phase

2

32 -32

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15

12 -12 (bj

16 -16

20 -20 Time 730 Ka

phase 3

24 -24 -24

28 -28

32 - 32 -32 - 32

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15

12 - 12

(4

Time 1.2 Ma phase 3

-24

I8 - 16

20 -20

24 - 24

28 -28

32 -32 -32

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15

12 - 12

(4

16 - 16

20 -20 Time 1.4 Ma

28 -28 - 28 phase

3

24 -24 - 24

32 -32

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15

. .

12 - 12

(el

16 -16

20 -20

24 -24

Time 1.6 Ma

28 -28 phase

4

32 -32

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15

12 -12 - 12

r---

16 -16

20 -20

24 - 24

28 -28 -28

32 - 32 -32

(f) Time 2.1 Ma

phase 4

-24

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15

Figures. The development of a diapir in four phases. Parameters used are phase diagram Type 1, density contrast 100kgm-3, and the other parameters listed in Table 2. The width and depth of the models are given in kilometres.

0 1995 RAS, G J I 123, 59-70

(7)

Modelling of melting processes and diapirism 65

amount of melt M temperature T chemical field X stream function qj -

12

-

12

-20 -20 -20

-

16 -16

- 24 -24 -24

-28 -28

32 _I_c_____

32 -32 -32

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20

rnax. 0.0326730 r n n x . rnclt granitc 0 2605 time steps 4000.00

time in ycars 1.64631et-06 rnax I'crnp 995 565 r i i i n . viscosily -5.30103

amount of melt M temperature T -

12

-

12

-

16 -16

-20 - 20

- 24 -24

- 28 -28

- 32 - 32

chemical field X stream function @

1 - 1

time steps 4000.00 rnax. S t r e a m f . 1.15537 r n a x melt g r a n i t e 0.2254 time in y e a r s 1.64996e+06 m l B x ~ c m p . 807.731 rriitn. viscosity -5.30103

(c) Density Variation

P s l i d granite

-

Pliquid granite

= 200 kg/m3

amount of melt M temperature T chemical field X stream function

<p

_ - - ~ - -

-12

-

12 - 12

-

16

-

16

-20 -20 -20

-

24 -24

-

24

-28 -28 -28

-32 -32 -32 1 _ _ - _

0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20

tirnc s t c p s 4000.00 max Slrearnf. 3.1 1698 n m x . rrielt g r a n i t e 0.2313 t i m e i n y c a r s 1.64955e+06 r r , a x . T ~ : ~ ~ , . ii 968.502 min. viscosity - 5 , : j O I O : j Figure 6. Influence of density variation on diapirism. Phase diagram Type 1 and the parameters listed in Table 2 are used.

were used in further model calculations. The variations in phase diagram resulted in changes in the melt concentration- temperature ( M - T ) curves, which are shown in Fig. 3(b).

This causes a variation of the maximum melt concentrations at the eutectic temperature. The buoyancy increases with increasing melt concentration, and as a first result nearly the same behaviour is observed as in the model series with the

density variations (Fig. 6). Thus a general rule can be formu- lated: the closer the starting composition of the granite is to the eutectic point, the higher the potential diapiric rise.

From a comparison of the emplacement depths of the diapirs shown in Fig. 7 it is possible to see that there is a small increase in the ascent height for models with increasing density differences. The influence of changes in the phase diagram 0 1995 RAS, G J I 123, 59-70

(8)

66 D. Bittner and H . Schmeling

comparison between diapir-models

phase diagrams

Figure 7. Different ascent heights are reached by the diapirs based on the variation of the density contrast between the solid and liquid granite and the phase diagram. The 'wavenumber' (the number of diapirs across the width of the model, 40km) is indicated above each height bar. During the ascent of the diapirs, the phase diagrams result in average melt concentrations of M

-

20-30 per cent (Type l), M

-

35-45 per cent (Type 2), M

-

75-85 per cent (Type 3).

from a small melt fraction to a large melt fraction, however, results in a slight decrease in the ascent height. (Only models which resulted in diapirism are discussed here.) This is a surprising result as a larger volume fraction of melt leads to greater buoyancy. A closer analysis of the models, however, shows that with increasing eutectic behaviour the melt concen- tration strongly increases, while the thickness of the partially molten layer decreases. This leads to a decrease in the wave- length of the instability (see wavenumbers indicated in the diagram). The resulting smaller characteristic wavelength leads to smaller diapirs, which cool more efficiently and, if compared to large diapirs, build up stress fields with smaller amplitudes.

The resulting higher viscosity of the surrounding rock impedes high emplacement levels.

R I S I N G D I A P I R S OR H O R I Z O N T A L L A Y E R I N G

The density contrast between liquid and solid granite may cause a RTI. The typical result of such a RTI are rising diapirs (e.g. Ramberg 1981). Our results show that there are also solutions in which diapirs do not form. The question arises, which mechanism controls whether diapirism occurs or whether the partially molten rock remains horizontally layered.

A RTI can only develop in the presence of granitic melt with a density lower than that of solid granite. Thus the lifetime of granitic melt can be used to formulate a first condition that must be fulfilled before diapirs can rise. The amplitude A of a RTI of two constant viscosity layers at the early stage is

A = A, exp(yt) (8)

with A,=amplitude of the disturbance at time zero and y = growth rate

APg 1

2%k ('1126,

+

62)('112"1+

4'

mi = .-7

ai=

.

y = -

2khi cosh(2khi)

+

1

sinh(2khi) - 2khi' slnh(2khi) - 2khi'

where h i = thickness of layer i; qi=viscosity of layer i; q I 2 = q1/q2; k = wavenumber. Subscript i = 1 refers to solid granite;

i = 2 refers to partially molten granite (Woidt 1980).

The maximum time tmelmax available for the ascent of diapirs during the lifetime of the granite melt is dependent upon the heat input Q, the depth and the temperature gradient. At any time t < tmeltmax the remaining time tmelt = tmelmax - t corre- sponds to a critical growth rate of yc = l/tmeLt. Based on eq. (8) this growth rate would result in an amplification of the original amplitude by the factor e. Criterion 1 can be formulated as follows: for diapirism to occur during the remaining time for which the granite melt exists, the actual growth rate must be larger than the critical growth rate, i.e. y > yc or t , > tmcltmax - t.

It is important to note that the critical growth rate is not constant, but is dependent upon the remainder of time for which a granitic melt exists, i.e. the critical growth rate increases as a function of time. The values of the actual growth rate and the critical growth rate at any time, based on criterion 1, are shown in Fig. 8. The largest value for tmeltmax in our models is 4Mio years. This results in an initial critical growth rate of However, as the creation of a melt front and the growth of an RTI are competing with each other, a second criterion must be fulfilled before the diapir grows. This criterion is related to the speed of growth of both fronts.

At the early stage of a RTI the growth velocity of a perturbation of the interface partial melt-solid having an amplitude A is given as

(9)

y c = 7.9 x 10-15 ijs.

udiapir = dA/dt = A,d(exp(yt))/dt = Ay ,

where A = amplitude of the RTI and y = growth rate.

The waviness of the instability is, however, smoothed and the amplitude is outpaced by a growing melt front, if its velocity is higher than that of the instability. The velocity of the melting front is the time derivative of the height Hmelt of the granite melting front,

Umelt = dHmeIt/dt. (10)

0 1995 RAS, GJI 123, 59-70

(9)

Modelling of melting processes and diapirism 67

g r o w t h r a t e gamma

(a) - 1 1

. . , . . . , . , , . , , , , .

10

- 1 6

'I

0

0.0

0.2

0.4

0.6 0.8

1.0

relative time

g r o w t h rate gamma

(b) l o - l l

10-

s I -

0.0

0.2

0.4 0.6 0.8 1.0

relative time

Figure8. Growth rate as a function of non-dimensional time. For scaling, the time has to be multiplied by 2.16 x lo6 years. The actual growth rate y (dashed line) must be greater than the critical growth rate g, (solid line) in order that a diapir can rise during the existence of granitic melt. Curves in (a) and (b) are from models Runl and Run2 shown in Figs 6(a) and (b).

We can summarize the two criteria to predict diapiric ascent:

(1) Y > Y ~ ; ~,=l/'(tme~tmax--), (2) udiapir >

They are now tested with the help of the numerical model.

The parameters for the two velocities are obtained from the numerical results of each model. Owing to the dynamical behaviour of the models, average values are taken.

Owing to the variation in stress- and temperature-dependent viscosities in our model there are a number of coupled param- eters that can influence the growth rate curve. On one hand the viscosity is dependent upon the strain rate and therefore the strength of convection (i.e. diapirism), while on the other hand the growth rate of the RTI is controlled by the viscosity.

The resulting interrelation between the strength of convection and the growth rate of the RTI is responsible for the behaviour of the y curve for stress-dependent rheology.

In the numerical calculations of the velocity from eqs (9)

and (lo), the amplitude A remains a model-dependent free variable and can be adjusted to fulfil our two criteria. This amplitude was chosen to be 150m for a granite layer of 15 km, ensuring that the point of intersection of both velocity curves and the start of diapiric ascent coincides in our models.

Figures 8 and 9 show the results when criteria 1 and 2 are applied to the models. The first criterion y > y c is fulfilled in Runl (Fig. 8a) for only 1 Ma. Diapirism is possible for more than 2 Ma in the model Run2 shown in Fig. 8(b). In the model without diapirs (Runl, and Run4, Figs 9a and d ) there are no time ranges in which both criteria are fulfilled at the same time.

In the models Run2, Run3, Run5 and Run6 (Figs 9b, c, e and f ) , there are time regions in which both criteria are simultaneously fulfilled and diapiric ascent occurs within this time range.

There is a peak at early times for the RTI velocity to be seen in the velocity curves. This is caused by the speed at which convection starts in the basalt at the beginning of the calculation. Owing to the coupling between the strength of convection, the growth rate and the stress-dependent viscosity, the growth rate increases very quickly. This is followed by a decrease in both the growth rate and the strength of convection.

Shortly afterwards the velocity of the melt front reaches a constant value for a short while. The increased heat transport within the basalt during the early convection peak leads to an increase in the heat flow into the granite. As a result, the decrease in the growth velocity of the melt front is retarded for a short time.

With the two criteria discussed above, we now have the tools necessary to predict the occurrence of diapirism or the growth of horizontal layers of partially molten granite. The time-scale at which diapirs begin to ascend can also be calculated.

D I S C U S S I O N

As previously mentioned, there are a number of possible scenarios for the heating and melting of the lower crust. We wish to compare two examples and discuss their characteristic heat balances. The first example is the present one of underplat- ing by molten basalt. Based on the data used in the present models, an internal energy of E = 5.5 x 1OI8 kJ is stored in the basalt underlying the granite in the lower crust at a temperature of 760 "C. The model shows that this energy is transported to the granite within maximally 1 to 1.5Ma. The second case that we wish to discuss envisions the rise of hot unmolten asthenosphere of 1300°C up to Moho level. As mentioned above, this could possibly be the result of delamination of a lithospheric root during orogeny. For reasons of comparison we assume the asthenospheric layer of peridotite to have a thickness of 5 km, as in the case of basaltic underplating. The warm peridotite in this model has a heat content of E = 5.6 x 10l6 kJ compared with granite at 760 "C. As no convec- tion occurs in the solid peridotite, the time-scale on which the energy is transported to the granite is estimated by a character- istic diffusion time to be 2 to 2.5Ma, thus it is only slightly longer than in the case of basaltic underplating.

These calculations show that a 1300°C layer of crystalline peridotite will create a similar amount of granite melt as will underplating by basaltic magma. The realistic natural scenario is probably a mixture of these two extremes.

An important variable in the ascent of granitic diapirs is the 0 1995 RAS, G J I 123, 59-70

(10)

68

D.

Bittner and H. Schmeling

0.030. ' ' ' ' ' ' . ' ' ' ' ' . ' . ' . ' . 0.030r'

' ' " ' " ' ' 'I. ' ' ' ' ' ' . ' ' . 1. . .

-

0.020

i ;

\ E

Y

> 0.010:

emplacement depth, which may range between a few kilometres below the surface (Wilson 1989) to 10-15 km below the surface (Barton et al. 1988; Miller, Watson & Harrison 1988). The diapirs in our model reach a maximum height of 15 km (see Fig. 7) below the surface. This is the depth of the critical temperature of c. 500 "C below which the material can no longer flow. The deviatoric stress created by diapiric ascent at this depth is approximately 10MPa. This stress decreases as the surface is approached. Tectonic stresses up to 100MPa may be present in tectonic regimes like orogens or other plate boundary environments which are not included in our model.

These stresses will decrease the modelled rheology by a factor Aq = lo3, and based on an activation energy of granite of E = 150 kJ mol-' the critical temperature for the ductile-brittle transition will be decreased by c . 200 "C. Thus tectonic stresses decrease the effective viscosity, allowing granitic intrusives to be emplaced by diapiric ascent possibly to depths of up to

10-12 km.

In summary, our model shows that, based on realistic rheological and thermodynamic assumptions, diapirism may be an important mechanism within the ductile lower crust.

However, it is unlikely that any further ascent into shallower crustal levels may occur by diapirism. To explain shallower emplacement depths, other mechanisms like dyke propagation

0.030.' ' ' * ' . '

;

-

L 0.020 :

sh

\ E

> 0.010: ;

+.

I

(Clemens & Mawer 1992) or stoping (Marsh 1982) might take over.

Diapirism is possible in the rheological setting employed in the above calculations. This rheology is, however, the weakest that has been observed in experimental studies of crystalline granite (Kirby & Kronenberg 1987). At significantly higher subsolidus viscosities and using the same realistic values of density, etc., no diapiric ascent through solid host rock is expected. As diapiric deformation is observed in both granitic plutons and surrounding host rock (Brun & Pons 1981; Cruden 1988; Ramsay 1989) it may be concluded that the lower crust might have a weak rheology (Miller et al. 1988).

None of the models presented here showed significant entrainment of basaltic material into a rising granitic diapir.

Entrainment is controlled by both the density contrast and viscosity contrast at the basalt-granite interface. The density of both molten and unmolten basalt is more than 10 per cent higher than that of granite with the presently chosen param- eters. According to an empirical law based on isothermal models (Cruden, Koyi & Schmeling 1995), this density contrast would prohibit significant entrainment at any viscos- ity contrast. This prediction is confirmed now by our thermodynamically consistent models.

Our models show that diapirism is convection-driven by the

0.0307' ' " ' ' " ' ' . ' . ' ' ' ' ' '

-

0.020 ::

[

::

\ E

-

' 0.010:

.,

.

..

0 1995 RAS, G J I 123, 59-70

0.030.' ' " ' ' ' ' ' ' ' . . ' " ' . I ' . ~ ~ l ~ ' ~' . l - ~ ' '

-

0.020 ::

-

0.020 ::

i : i k

\ E

-

' 0.010:

.,

(11)

Modelling of melting processes and diapirism 69 solid-liquid phase transition. Granite particles may experience

multiple overturns and melting-solidification events, which should be detectable as zoned minerals. Another implication here is that heterogeneous lower crust could be homogenized by such a convective flow. Multiple overturns within diapirs and their deformation patterns have been discussed by Cruden (1988) and Schmeling, Cruden & Marquart ( 1988), being the result of viscous drag acting on a diapir rising over long distances. The present model suggests an alternative expla- nation leading to similar deformational structures but not requiring large distances of rise.

Our diapirism model is able to predict the ascent of large amounts of acidic melt (e.g. as has been inferred for the Variscan Orogeny), and with the introduction of tectonic stresses the ascent height of diapirs can be calculated in good agreement with geological observations.

It should be emphasized that in our model melting occurs as a consequence of heating from below. However, orogenic processes associated with crustal thickening might induce dehydration processes involving biotite or muscovite break- down at lower crustal depths. The release of water reduces solidus temperatures and melting might occur at isothermal conditions (Johannes & Holtz 1991). Increasing the lithostatic pressure within a thickening crust might be described by a pressure front penetrating the crust from below in the Lagrangian frame of reference. It has to be studied in future models if the effect of a moving pressure front is comparable to the effect of a temperature front as studied in this paper.

This numerical model is, however, only a predictive tool. In order to successfully model the natural processess occurring in the Earth it is necessary to better understand the relation- ships between experimental data, and modelled parameters and mechanisms.

CONCLUSIONS

We have modelled melting processes and diapiric flow within the lower crust in a fully self-consistent manner, i.e. by solving the energy and flow equations coupled with phase relations and rheological equations of state. The model leads us to the following conclusions.

( 1) Assuming realistic thermodynamic and rheological par- ameters, depths of emplacement of granitic diapirs of 15 km are possible, associated with diapiric stresses of up to 10 MPa.

In the presence of tectonic stresses of the order of 100MPa, weakening of the stress-dependent rheology may result in emplacement depths which are 3-5 km shallower.

(2) The model suggests a brittle-dutile transition tempera- ture of about 500°C. Below this temperature no viscous flow is observed.

(3) Convection driven by solid-liquid phase transformation may result in multiple overturns of granitic material, leading to multiple loops of the P T t path and multiple melting- solidification events within a rock unit.

(4) Viscous flow is not confined to the partial molten granite, but extends deeply into the unmolten surrounding rock.

( 5 ) Above the 300 "C isotherm, lateral temperature vari- ations are small (even in the models including the whole crust), making it difficult to identify active diapirism within the lower crust by heat-flow measurements at the surface.

A C K N O W L E D G M E N T S

This work profited significantly from stimulating discussions with Volker von Seckendorf, Don Dingwell and Sharon Webb.

Reviews by Alexander Cruden and Giancarlo Serri have been very helpful. We gratefully acknowledge financial support by Deutsche Forschungsgemeinschaft (Grant Schm 872-1).

R E F E R E N C E S

Barton, M.D., Battles, D.A., Bebout, G.E., Capo, R.C., Christensen, J.N., Davis, S.R., Hanson, R.B., Michelson, C.J. & Trim, H.E., 1988.

Mesozoic contact metamorphism in the western United States, in Metamorphism and Crustal Evolution of the Western United States, Vol. 7, pp. 110-178, eds Ernst, W.G. & Rubey, Prentice Hall, Englewood Cliffs, NJ.

Bergantz, G.W., 1989. Underplating and partial melting: Implications for melt generation and extraction, Science, 245, 1093-1095.

Bird, P., 1978. Initiation of intracontinental subduction in the Himalaya, J. geophys. Rex, 83(B10), 4975-4987.

Brid, P., 1979. Continental delamination and the Colorado plateau, J. geophys. Res., 84(B13), 7561-7571.

Brun, J.P. & Pons, J., 1981. Strain patterns of pluton emplacement in a crust undergoing noncoaxial deformation, Sierra Morene, Southern Spain, J . struct. Geol., 3(3), 219-229.

Buck, W.R. & Toksoz, M.N., 1983. Thermal effects of continental collisions: thickening a variable viscosity lithosphere, Tectonophysics, 100, 53-69.

Clemens, J.D. & Mawer, C.K., 1992. Granite magma transport by fracture propagation, Tectonophysics, 204, 339-360.

Cruden, A.R., 1988. Deformation around a rising diapir modeled by creeping flow past a sphere, Tectonics, 7, 1092-1101.

Cruden, A.R., Koyi, H. & Schmeling, H., 1995. Diapiric basal entrainment of mafic into felsic magma, Earth planet. Sci. Lett., in press.

De Yoreo, J.J., Lux, D.R. & Guidotti, C.V., 1989. The role of crustal anatexis and magma migration in the thermal evolution of regions of thickened continental crust, in Euolution of Metamorphic Belts, pp. 187-202, eds Daly, J.S., Cliff, R.A. & Yardley, B.W.D., Geological Society Special Publication No. 43.

Dingwell, D.B., Knoche, R. & Webb, S.L., 1992. The effect of P20, on the viscosity of haplogranitic liquid, J . Mineral., 5, 133-140.

Huppert, H.E. & Sparks, R.S., 1988. The generation of granitic magmas by intrusion of basalt into continental crust, J . Petrol., 29(3), Kirby, S.H. & Kronenberg, A.K., 1987. Rheology of the lithosphere:

selected topics, Reo. Geophys., 25, 1219-1244.

Kukowski, & Neugebauer, 1990. On the ascent and emplacement of granitoid magma bodies-dynamic-thermal numerical models, Geol. Rundsch., 79( 1). 227-239.

Lathop, A., Blum, J. & Chamberlain, C., 1994. Isotopic evidence for closed systems anatexis at midcrustal levels: An example from the Acadian Appalachians of New England, J. geophys. Res., 99, Lorenz, V. & Nicolls, LA., 1984. Plate and interplate processes of Hercynian Europe during the late Paleozoic, Tectonophysics, 107, 25-56.

Mahon, K.I. & Harrison, T.M., 1988. Ascent of a granitoid diapir in a temperature varying medium, Geophys. Res., 93( 32), 1174-1 188.

Marsh, B.D., 1982. On the mechanics of igneous diapirism, stopping and zone melting, Am. J . Sci., 282, 808-855.

Matte, P. 1986. Tectonics and platetectonics model for the variscan belt of Europe, Tectonophysics, 126, 329-374.

McKenzie, D., 1985. The extraction of magma from the crust and mantle, Earth planer. Sci. Lett., 75, 81-91.

Miller, C.F., Watson, M.E. & Harrison, T.M., 1988. Perspective on 599-624.

9453-9468.

Q 1995 RAS, G J I 123, 59-70

Referenzen

ÄHNLICHE DOKUMENTE

The petrophysical properties of the fractured granitic reservoir and the Buntsandstein and Muschelkalk units, sampled from exploration borehole EPS-1 at the Soultz-sous-Forˆ

Abstract At Montagna della Maiella and at Gola del Furlo (central Apennines) two discrete layers of bentonic clay are intercalated within the pelagic (Furlo) and tur-

We may thus conclude that both the viscosity and ro- tation suppress the instability of the superposed grav- itating streams when the streams rotate about an axis in the

• Implementation of direct environmental burdens of the energy system Solved (example: air emissions of power plants). • Implementation of indirect environmental burdens of the

Control electrical components Microcontroller Programming of microcontroller USB serial programmer Transmit ultra sound signal Ultra Sound Transducer Receive ultra sound signal

The time-change term in the continuity equation implies that this cannot be accomplished using measured surface profiles for direct input (as in the case of

– Total amount of all compensations paid to persons related to members of the BoD, executive board and advisory board; the names do not need to be disclosed. Compensations of

In order to assess the effects of preferential flow on soil water flow and streamflow in forested watershed, precipitation, amount of preferential flow, infiltrated flow and