- Email: [email protected]
Linear and weakly nonlinear problems of the theory of thermal convection in the Earth's mantle
92
Physics of the Earth and Planetary Interiors, 50 (1988) 92-98 Elsevier Science Publishers B.V., Amsterdam — Printed in The Netherlands
Linear and weakly nonlinear problems of the theory of thermal convection in the Earth’s mantle Boris I. Birger The Earth’s Physics Institute, U.S.S.R. Academy of Sciences, 10 Bol’shaya Gruzinskaya, Moscow 123810 (U.S.S.R.) (Received February 18, 1987) Birger, B.I., 1988. Linear and weakly nonlinear problems of the theory of thermal convection in the Earth’s mantle. Phys. Earth Planet. Inter., 50:92—98. The thermoconvective instability of the upper thermal boundary layer formed by the large-scale convection of the mantle is analysed. The creep in the mantle is assumed transient when strains are small. The linear hereditary medium is accepted as a rheological model which describes the mantle’s transient creep. At this rheology the instability is oscillatory and the supercriticality, as suggested by the postglacial rebound data, is small. The application of nonlinear stability theory permits derivation of an asymptotic solution in the form of a steady convective oscillation the amplitude of which is independent of the initial conditions. The supercriticality is considered as a small parameter and perturbation methods are used. The amplitude and the frequency of oscillation are found as functions of the supercriticality.
1. Introduction The simplest model of the mantle rheology is a Newtonian fluid in which the stress is proportional to the strain rate. This model is usually used when thermal convection in the Earth’s mantle is studied. The Newtonian model is, however, an oversimplification, the true rheology of the mantle being much more complex. Laboratory studies of rock deformation at temperatures and pressures approximating mantle conditions have shown that, firstly, steady creep is reached at a creep strain of several percent, secondly, strain rate depends nonlinearly on stress at steady creep whereas this dependence is linear at transient creep (Goetze and Brace, 1972; Murrell and Chakravarty, 1973). Thus, if creep strains do not exceed several percent during a dynamical process, there is transient creep. Only such processes will be considered in the present work. The first of them is postglacial rebound. Weertman (1978) pointed out that strains in the mantle related to ice loading are estimated to be considerably smaller than 1O~and, therefore, transient creep is involved in the postglacial rebound phenomenon.
The second process considered here is a smallamplitude convective oscillation. The large-scale convective circulation which causes plate creation and destruction is associated with large strains in the mantle and steady creep. It was investigated by Parmentier et al. (1976) and Christensen (1984) who used a non-Newtonian fluid model in which a power law relates the stress to the strain rate. Large-scale convection forms the thermal boundary layer near the Earth’s surface which moves with constant velocity. Within this layer, there is a steep temperature gradient but no strain. Parsons and McKenzie (1978) and Yuen et al. (1981) investigated the instability of the boundary layer within the framework of the Newtonian fluid rheological model. In the present work the boundary layer instability will be analysed within the framework of a rheological model describing transient creep of the mantle.
2. The governing equations The upper boundary layer of a large-scale convection cell is modelled as a homogeneous infinite
93
horizontal layer. Its physical properties are given by the following parameters: density p 3.5 g cm3 thermal diffusivity ic 10—2 cm2 s’; thermal expansion a 3~iO~ K~ gravitational 2 shear modulus p~ acceleration g thickness iO~cm sof the layer d iO~cm; 1012 dyn cm2 temperature gradient A 10 ~ K cm~. The origin of the coordinate system is situated on the upper boundary plane of the layer. The z-~.xisis directed vertically upwards. The y-axis is parallel to the direction of the large-scale flow and only perturbations independent of y are considered. In this case the constant velocity of the boundary layer is eliminated from the equations. The unearized equations governing the stability of the basic state of the layer are ap 8p —
a~, ao
3T~~
aT
3T,, xz+_~~~+Ra0~=0
aT
+
8v
—
where ejj is the deviatoric strain tensor, K( 1) the creep integral kernel. The first term on the righthand side of eq. 5 gives the elastic strain and the second one gives the creep strain. In laboratory studies of rock deformation the creep function is usually sought. It is defined as the creep strain normalized against the elastic strain at the constant stress applied from t = 0. For the mantle rocks at temperatures close to the melting point the creep function is given by the known Andrade law
(1)
f(t)=t’yf3
(2)
where y and /3 are dimensionless rheological parameters. If we take = const in eq. 5, we
(3)
ía2 + az2) a2~
creep. Then the constitutive equation which must be added to the set of eqs. 1—4 may be written in the form + J1K(l t 1)~~(t1)dti} (5) e.~= 2
(4)
where 0, v,, p, i~ are the nondimensionalized perturbations of temperature, velocity, pressure, and of the components of the deviatoric stress tensor. The scale of space coordinates is d = iO~ cm, time scale d2/ic = 1016 s, velocity scale ic/d 109cm s’, pressure and stress scale ~s= 1012 dyn cm2, temperature scale Ad= iO~ K. Ra is a dimensionless parameter
(6)
readily obtain df K(t)—— (7) dt Hence, the creep kernel associated with the Andrade law is
(
K t) = yt5 ~//3 (8) The rheological model given by eqs. 5 and 8 we shall call the Andrade medium. The upper surface of the layer is taken to be isothermal and traction-free, that is, the perturbation of temperature and the force per unit area disappear on it. These boundary conditions reduce
Ra=pgaAd2/ 1.t _ —
which is estimated as Ra i0~. If we 2/Kp, introduce then the reference kinematic viscosity v = ,.5d the definition of Rayleigh number Ra is identical to that used for a Newtonian fluid. In writing eqs. 1—4 we neglect the inertia terms in the momentum eqs. 1 and 2, since the Prandtl number is very large, and use the Boussinesq approximation (the validity of this approximation for the rheological model used in the present work has been studied by Birger, 1977). The linear hereditary rheological model is accepted as a description of the mantle transient
to z=0;
0=u~,
T~2=0,
—p+T~,+~u~=O
(9) where u~is a small vertical displacement of the surface, ~ = pgd/~tsthe dimensionless parameter equal to 3.5 10~ for the layer considered here. The lower boundary of the layer is supposed to be ‘free’, i.e. the vertical velocity, the tangential stress, and the temperature perturbation vanish on it .
Z =
—
1;
~
=
0,
T~
2 = 0, 0= 0 (10) Free boundary conditions for the lower surface simplify the analysis permitting us to obtain an
94
analytical solution without significantly changing the results. In the next sections only perturbations periodic in x will be considered
where u~,(0, 0) is the initial displacement of the upper surface (z = 0) at t = 0. Since the displacement u~(0,t) is related to the velocity by v~= du~/dt,the relation of their Laplace transforms is
b1 =b1(z,
v(0, s) =su’(O, s) —u~(0,0)
t)
cos kx,
b2=b2(z,
t)
sin kx (11)
where k is the dimensionless wavenumber, b1 is one of variables v~,0, p, ~ ;~ b2 is one of v~, ~
3. The estimation of rheological parameters by postglacial rebound data Representing the basic eqs. 1—5 in terms of their Laplace transforms and eliminating all the variables except for the vertical velocity, we find [(i
—
~~)sH*(s)L4
—
Rak2I ~*(z, s)
2O(z, 0) (12) Rak where the asterisk denotes the Laplace transform, s is the Laplace transform variable, 0( z, 0) the initial temperature perturbation at t = 0, L the differential operator L2 = D2 k2, D = dz * (s) relates to the creep kernel Laplace trans=
—
form by 11*(s)
=
1/s[1
(
K*(s)]
+
(13)
—
(~~)
where 11(t) is the memory function. Thus, 1~I*(s) in eq. 12 is the Laplace transform of the memory function. The boundary conditions given by eqs. 9 and 10 become 2 + k2) ~ = 0 z = 0; (D [sH*(D3 3k2D) ~pk2] v* = ~pk2u~(0, 0) —
—
—
Rak2]v~*= Rak2 u
[sH*L
2(0, 0)
(15)
Such a form of the initial conditions corresponds to the known problem of harmonic ice loading. Solving the problem given by eqs. 12, 15, 16 and 19 we find v~*(z,s) and then, using eq. 17, we obtain uo2kB(k)q.,_hH*(s) (20) U’ (0, s) = (1 + 2kB k) ~ ‘s11 * (s))
(
—
where B(k) = (2k + sinh 2k)/(cosh 2k 1) (21) Since the parameter Ra/i~= aAd 3• 10-2 is small, the temperature gradient present in the layer can be neglected in the postglacial rebound —
v~=D
Theis final step of inverting transform extremely difficult. Wethe shallLaplace seek only an asymptotic solution at large t using the known theorem of asymptotic expansion of the pre-image (see, e.g., von Doetsch, 1967). The Laplace transform given by eqs. 20 and 18 has the singularity s = 0. We expand the right-hand side of eq. 20 in the vicinity of s = 0 in a power series of s. According to the theorem, to find an asymptotic solution it is enough to invert the first term of the expansion. So we obtaindisplacement the asymptotic dependence of the upper surface on time u~(0,t)
=
u
02kB(k)~jf~’$t~/(yF(y)I’(1—y)) (22) For a Newtonian fluid the Laplace transform of the dimensionless memory function reduces(23) to 2
2v~*=D4v z= —1;
For the Andrade medium, as follows from eqs. 8 and 13 K*(s)=YF(y)/$sY, (18) fl*(~)=$s~’/(yF(y) +/3s~) where F(y) is the gamma-function of y. Now consider initial perturbations of the form 0(z,0)~0, u~(0,0)=uo (19)
problem.
Note that the constitutive eq. 5 can be replaced by = 2f11 t t 1) de~~ (14)
4
(17)
2~~=0
(16)
H*(s)=vpic/,5d
95
where p is the kinematic viscosity. Substituting eq. 23 into eq. 20 and inverting the Laplace transform, we obtain the known solution of the post-
gives ~~(z)
glacial rebound problem for the Newtonian rheological model of the mantle u~(0,t)=u0exp(—t/T) V
T =
K
2kB(k)~ ~
(24) (25)
where T is the relaxation time nondimensionalized 2/K. Unlike eq. 22 obtained for large t, eq. 24 is valid for arbitrary t. As follows on the time scale d from eq. 24, the relaxation time is the time at
=
2e~~(z)f tH(
—
‘~)exp[iw(t1
—
t)]
=
2ê1j(z)f
0
2)e”’12 dt2
(28)
where d,1(z) is the complex amplitude of the deviatoric strain rate tensor. As follows from eq. 28, the solution of the basic equations may have the form (27) only at large t when we can take this case = ~ as the upper limit of the latter integral. In T,~(z) =
2F(w)è 11(z)
which the upper surface displacement decreases by a factor e. Although in the framework of the Andrade model the dependence of displacement on time is not exponential, we define the relaxation time in this case by the equation u0/u~(0, T) = e too. Then, as follows from eq. 22 2kB(k)/1’$e/(’fI’(’f)F(l —v)) (26) The laboratory studies of rock deformation show that y= 1/3 y.is Equation a typical value of the us rheological parameter 26 permits to estimate the rheological parameter /3 at the fixed value of -y with the help of the postglacial rebound data. When the dimensions of ice load are nearly 100 km (k = ir) the relaxation time is 6 X 10~s (T 6~106). Substituting ‘f = 1/3, k = ~r, T 6~106 into eq. 26 we find that /3 5. i0’. For the same values of k and T eq. 25 gives v 1020 cm2 s~. This known estimate may be considered as the layer effective viscosity for postglacial rebound processes but cannot be used for investigation of the other geodynamical processes if the true rheology of the mantle is not Newtothan. =
4. Thermoconvective instability
(29)
where
F(c~)=
j
l1(t)e”
dt = 11*(iw)
As follows from eq. 29, the complex function F( ca) is an analogue of the viscosity. Hence, F( ~,) is naturally called the ‘complex viscosity’. For the Andrade medium the complex viscosity is
1/(’fr(y)
+ /3(iw~) (30) F(w) = Substituting eqs. 11 and 27 into the basic equations and boundary conditions given by eqs. 9 and 10 and eliminating the complex amplitudes of all physical variables except for the vertical velocity amplitude, we obtain a boundary-value problem given by the ordinary differential equation
$(1wy
1—
=0
1~—)iwF(ca)L4_Rak2]vz(z)
(31)
The conditions on the upper boundary are ~ = 0; (D2 + k2) v~= 0
[icaF(w)(D~ 3k2D) —
—
i~tik2]Vr
=
0
(32)
[IWF(w)L4_Rak2]vz=O
On the lower boundary z = —1; v~= D2v~= D4v 2= 0
In this section we shall seek a solution for the basic eqs. 1—5 in the form of an oscillation t (27) b1(z, t) = bi(z)e~~~t,b2(z, t) = b2(z)e” where the meaning of b 1, b2 is the same as in eq. 11. The dimensionless frequency a is a complex number. Substitution of eqs. 11 and 27 into eq. 14
dt1
tH(t
(33)
where L and D are the differential operators introduced in the previous section. The problem given by eqs. 31—33 is considered as an eigenvalue problem. The eigenvalue is the critical Rayleigh number which gives the stability criterion. When Ra = Racr the solution of the problem has the form (27) with a real ca. Suppose
96
that the frequency which is sought satisfies the restriction << (34~
Replacing the complex viscosity (30) by (35) means that we neglect the elasticity of the medium for the frequencies considered. Under free boundary conditions the eigenfunctions of the problem are
where v~(z) is given by eq. 36 in which n = 1. In the present section the asymptotic solution of the basic equations in the form of a thermoconvective oscillation has been obtained without use of the Laplace transform. Besides the singularity s = 0 considered in the previous section, the boundary-value problem given by eqs. 12, 15 and 16 has the singularity s = ia associated with the thermoconvective oscillation. The method based on the use of Laplace transforms allows us to find the amplitude of asymptotic oscillation as a function of the initial perturbations 0(z, 0), u.(0, 0). But this advantage of the Laplace transform method is not very important because it can be shown that at any initial perturbations the corre-
v~(z)= C sin nirz
sponding amplitudes differ by a constant factor.
‘
/
Then the upper traction-free boundary behaves as a free one and the eqs. 33 are valid at z = 0 too. Besides, under the restriction (34), since yF(y) 1 for y = 1/3 and ~p<< 1, we can write F(c~)as ,~
\
\5—1 / /3~:w) /yyl’yy))
~35j\
~.
=
~ = 1, 2, 3,
..
(36)
where C is an arbitrary factor. The dispersion relation is given by 2
+ ~2~2) +
F(
)( k2 +
2~ 2)3
—
Racr k2
=
0
(37)
As follows from eq. 37
~)
Ra cr = (k2 +k2 2~2) IReF( F( w) = Im F(“ ~ (k2 + n2ir2) ReF(~~)
2
—
(38) (39)
—
—
-
3 =
F( “m)
‘P
.
hereminimal the is nearlycritical two orders Rayleigh of magnitude number. larger Therefore, than is fornecessary investigation to solve of convective numerically flows the strongly in the layer nonit .
where Re denotes ‘the real part of’, Im is ‘the imaginary part of’. For F(w) given by eq. 35, we find that the critical Rayleigh number reaches its minimum Ram at n = 1, k = km = ~‘r/(1+ y)1/2 For y = 1/3 we obtain km 2.7; cam = ca(km) 30; Ram 144 /3. If Ra = Ram, there is a regime of boundary instability for which only the convec~tive oscillation with wavenumber km = 2.7 (the wavelength (2iT/km)d is 200 km) and frequency = 30 (the period 2TT/Wm d2/K is nearly 108 years) is not attenuated. It is easy to verify that the frequency equal to 30 satisfies the restriction (34) which was a priori assumed. Then the amplitude of vertical displacement of the traction-free upper surface can be found readily from the second equation of set (32) (0)
.
to be 1020 cm2 s’, as shown above. Then the actual Rayleigh number for the layer considered
3
~.
.
5. Finite-amplitude thermoconvective oscillation If a Newtonian fluid is accepted as the mantle’s rheological model, the kinematic viscosity proves
~ 3D ~ k~ —
Z
(40)
linear equations. But for the Andrade rheological model the situation is greatly changed: the study of convection in the layer leads to linear or weakly nonlinear problems. According to the results obtained in the previous sections, y = 1/3, = 5. iO~ may be accepted as the Andrade rheological parameters for the upper mantle boundary layer. At these values of y and /3 the minimal critical Rayleigh number Ram = 144/3 7 iO~. The actual Rayleigh number for the layer was estimated to be Ra iO~. At Ra
.
97
linear stability theory predicts the unlimited increase of oscillation amplitude in time. But for sufficiently large amplitudes the use of linearized equations is not valid. The application of nonlinear theory permits us to obtain an asymptotic solution in the form of a steady thermoconvective of initial conditions. To describe perturbaoscillation the amplitude of which finite is independent tions in the frame of the Boussinesq approximation we must replace the partial derivative of temperature with respect to time in eq. 4 by the total derivative. Thus, an additional nonlinear term ~~(ao/ax) + ~~(80/8z) appears in the heat equation only. Introduce the supercriticality as 1 Ra
Ram Ra~
The next problem is posed by the coefficient of The solvability condition of the problem gives ca~= 0. Then solving the problem we obtain v(2) = ~(2) = 0 2
B 2F(ca k~)2 0) 2 sin 2irz 2 ~ Ra m k m B2F(ca 0) e2~~~~tl + c.c. (45) + 41~2 + 2ica0
~(2)
=
[I
(~2 +
~
I
—
j
Consideration of the problem posed by the coefficient of ~ leads to the solvability condition F( ca0) d F I (icao + ~2 + k~) + I
[
J
—
dca I~=~
—
+ica2F(w0)+— ~,2 1B12[
(41)
F(ca0) + 2ica0
4~2
The supercriticality is assumed to be small. We use it as a small parameter to expand the perturbations. The vertical velocity is sought in the form 2v~2~(z, x, t) = v~’~(z, x, t) + x, t) + (42) where v~(z, x, t) are periodic in x and t. These functions are expressed as Fourier series. The period in x is fixed. It equals 2 iT/km where km ~S found in the linear stability theory. Similar expansions are assumed for the temperature and the ...
other variables. The frequency is expanded in e too = ca 2ca 0 + ca1 + e 2 + (43) ...
When we substitute these expansions into the equations describing finite perturbations and equate the coefficients of “ to zero, we obtain a hierarchy of equations which must be solved successively, The coefficient of e yields a problem which is identical with that of the linear stability theory. Hence, ca0 cam. But now we write the solution of the problem as v(1)
=
sin
rrz
cos
kmX +
c.c.
(44)
where c.c. denotes the complex conjugate, B is an amplitude factor which must be found, ca is given by the expansion (43).
+
ReF(
ca0)J
2~2
—
—0
(46)
Substituting complex viscosity given by eq. 35, y= 1/3, km the = 2.7, and ca 0 = 30 into eq. 46, we find that ca2 —40, B I 36. Thus, the amplitude of vertical velocity oscillation, which is arbitrary in the linear theory, proves to be 36. I
—
Solution (45) obtained from the equations of the second-order approximation in c permits us to find the relation of Nusselt number to supercriticality
(47)
2
Nu=1+a For the Andrade medium at y= 1/3 the coefficient a is given by a 2.4. The solution for the linear stability problem in the form of a thermoconvective wave has been found by Birger (1978, 1980). In the present work the solution in the form of a thermoconvective oscillation, i.e. standing thermoconvective wave is obtained. Unlike the linear theory, solutions haying the form of travelling and standing waves do not relate to each other trivially in the nonlinear theory. If we write the solution of the first-order approximation in wave v~ = Z
B sin 2
srz
exp
in the form of a travelling
i(cat + kmX) +
c.c.
(48)
98
then, using the method developed in the present section, we find that ca2 = 0, I B I 24 (Birger, 1978). —
6. Conclusions Geology testifies that the vertical movements of the Earth’s crust have occurred everywhere during the whole geological history of the Earth. The major characteristic of the vertical movements is a periodic change in sign of movement which caused them to be called oscillatory (Beloussov et al., 1974). On the Earth’s surface there always exists simultaneously areas of movement in different (upward and downward) directions. The lowestfrequency vertical movements of the Earth’s crust in the continental platform regions are supposed to have periods of 108_109 years. The regions which rise (or subside) have dimensions of several hundreds of kilometres. The amplitudes of crustal oscillations in the platform regions are about 1 km. The results of the present work permit us to assume that the cause of the lowest-frequency oscillatory crustal movements is thermoconvective oscillation of the upper thermal boundary layer in the mantle, —
References Beloussov, V.V., Reisner, G.I., Rudich, E.M. and Sholpo, V.N., 1974. Vertical movements of the Earth’s crust on the continents. Geophys. Surv., 1: 245—273.
Birger, B.I., 1977. Complex viscosity and Boussinesq approximation in analysis of convective instability of upper mantle. DokI. Akad. Nauk SSSR, 235: 794—797 (in Rus. sian). Birger, B.I., 1978. The finite amplitude thermoconvective waves in the Earth’s mantle. DokI. Akad. Nauk SSSR, 242: 559—562 (in Russian). Birger, B.I., 1980. Thermoconvective waves in the Earth’s mantle. Phys. Earth Planet. Inter., 22: 238—243. Christensen, U., 1984. Convection with pressure- and temperature-dependent non-Newtonian rheology. Geophys. J.R. Astron. Soc., 77: 343—384. Goetze, C. and Brace, W.F., 1972. Laboratory observations of high-temperature rheology of rocks. Tectonophysics, 13: 583—600. Murrell, S.A.F. and Chakravarty, S., 1973. Some new rheological experiments on igneous rocks at temperatures up to 1120°C. Geophys. J.R. Astron. Soc., 34: 211—250. Parmentier, E.M., Turcotte, D.L. and Torrance, K.E., 1976. Studies of finite amplitude non-Newtonian thermal convection with application to convection in the Earth’s mantle. J. Geophys. Res., 81: 1839—1846. Parsons, B. and McKenzie, D.P., 1978. Mantle convection and the thermal structure of the plates. J. Geophys. Res., 83: 4485-4496. Von Doetsch, G., 1967. Anleitung zum praktischen Gebrauch der Laplace-Transformation und der Z-Transformation. Springer, München, 365 pp. Weertman, J., 1978. Creep laws for the mantle of the Earth. Philos. Trans. R. Soc. London A, 288: 9—26. Yuen, D.A., Peltier, W.R. and Schubert, G., 1981. On the existence of a second scale of convection in the upper mantle. Geophys. J.R. Astron. Soc., 65: 171—190.
Physics of the Earth and Planetary Interiors, 50 (1988) 92-98 Elsevier Science Publishers B.V., Amsterdam — Printed in The Netherlands
Linear and weakly nonlinear problems of the theory of thermal convection in the Earth’s mantle Boris I. Birger The Earth’s Physics Institute, U.S.S.R. Academy of Sciences, 10 Bol’shaya Gruzinskaya, Moscow 123810 (U.S.S.R.) (Received February 18, 1987) Birger, B.I., 1988. Linear and weakly nonlinear problems of the theory of thermal convection in the Earth’s mantle. Phys. Earth Planet. Inter., 50:92—98. The thermoconvective instability of the upper thermal boundary layer formed by the large-scale convection of the mantle is analysed. The creep in the mantle is assumed transient when strains are small. The linear hereditary medium is accepted as a rheological model which describes the mantle’s transient creep. At this rheology the instability is oscillatory and the supercriticality, as suggested by the postglacial rebound data, is small. The application of nonlinear stability theory permits derivation of an asymptotic solution in the form of a steady convective oscillation the amplitude of which is independent of the initial conditions. The supercriticality is considered as a small parameter and perturbation methods are used. The amplitude and the frequency of oscillation are found as functions of the supercriticality.
1. Introduction The simplest model of the mantle rheology is a Newtonian fluid in which the stress is proportional to the strain rate. This model is usually used when thermal convection in the Earth’s mantle is studied. The Newtonian model is, however, an oversimplification, the true rheology of the mantle being much more complex. Laboratory studies of rock deformation at temperatures and pressures approximating mantle conditions have shown that, firstly, steady creep is reached at a creep strain of several percent, secondly, strain rate depends nonlinearly on stress at steady creep whereas this dependence is linear at transient creep (Goetze and Brace, 1972; Murrell and Chakravarty, 1973). Thus, if creep strains do not exceed several percent during a dynamical process, there is transient creep. Only such processes will be considered in the present work. The first of them is postglacial rebound. Weertman (1978) pointed out that strains in the mantle related to ice loading are estimated to be considerably smaller than 1O~and, therefore, transient creep is involved in the postglacial rebound phenomenon.
The second process considered here is a smallamplitude convective oscillation. The large-scale convective circulation which causes plate creation and destruction is associated with large strains in the mantle and steady creep. It was investigated by Parmentier et al. (1976) and Christensen (1984) who used a non-Newtonian fluid model in which a power law relates the stress to the strain rate. Large-scale convection forms the thermal boundary layer near the Earth’s surface which moves with constant velocity. Within this layer, there is a steep temperature gradient but no strain. Parsons and McKenzie (1978) and Yuen et al. (1981) investigated the instability of the boundary layer within the framework of the Newtonian fluid rheological model. In the present work the boundary layer instability will be analysed within the framework of a rheological model describing transient creep of the mantle.
2. The governing equations The upper boundary layer of a large-scale convection cell is modelled as a homogeneous infinite
93
horizontal layer. Its physical properties are given by the following parameters: density p 3.5 g cm3 thermal diffusivity ic 10—2 cm2 s’; thermal expansion a 3~iO~ K~ gravitational 2 shear modulus p~ acceleration g thickness iO~cm sof the layer d iO~cm; 1012 dyn cm2 temperature gradient A 10 ~ K cm~. The origin of the coordinate system is situated on the upper boundary plane of the layer. The z-~.xisis directed vertically upwards. The y-axis is parallel to the direction of the large-scale flow and only perturbations independent of y are considered. In this case the constant velocity of the boundary layer is eliminated from the equations. The unearized equations governing the stability of the basic state of the layer are ap 8p —
a~, ao
3T~~
aT
3T,, xz+_~~~+Ra0~=0
aT
+
8v
—
where ejj is the deviatoric strain tensor, K( 1) the creep integral kernel. The first term on the righthand side of eq. 5 gives the elastic strain and the second one gives the creep strain. In laboratory studies of rock deformation the creep function is usually sought. It is defined as the creep strain normalized against the elastic strain at the constant stress applied from t = 0. For the mantle rocks at temperatures close to the melting point the creep function is given by the known Andrade law
(1)
f(t)=t’yf3
(2)
where y and /3 are dimensionless rheological parameters. If we take = const in eq. 5, we
(3)
ía2 + az2) a2~
creep. Then the constitutive equation which must be added to the set of eqs. 1—4 may be written in the form + J1K(l t 1)~~(t1)dti} (5) e.~= 2
(4)
where 0, v,, p, i~ are the nondimensionalized perturbations of temperature, velocity, pressure, and of the components of the deviatoric stress tensor. The scale of space coordinates is d = iO~ cm, time scale d2/ic = 1016 s, velocity scale ic/d 109cm s’, pressure and stress scale ~s= 1012 dyn cm2, temperature scale Ad= iO~ K. Ra is a dimensionless parameter
(6)
readily obtain df K(t)—— (7) dt Hence, the creep kernel associated with the Andrade law is
(
K t) = yt5 ~//3 (8) The rheological model given by eqs. 5 and 8 we shall call the Andrade medium. The upper surface of the layer is taken to be isothermal and traction-free, that is, the perturbation of temperature and the force per unit area disappear on it. These boundary conditions reduce
Ra=pgaAd2/ 1.t _ —
which is estimated as Ra i0~. If we 2/Kp, introduce then the reference kinematic viscosity v = ,.5d the definition of Rayleigh number Ra is identical to that used for a Newtonian fluid. In writing eqs. 1—4 we neglect the inertia terms in the momentum eqs. 1 and 2, since the Prandtl number is very large, and use the Boussinesq approximation (the validity of this approximation for the rheological model used in the present work has been studied by Birger, 1977). The linear hereditary rheological model is accepted as a description of the mantle transient
to z=0;
0=u~,
T~2=0,
—p+T~,+~u~=O
(9) where u~is a small vertical displacement of the surface, ~ = pgd/~tsthe dimensionless parameter equal to 3.5 10~ for the layer considered here. The lower boundary of the layer is supposed to be ‘free’, i.e. the vertical velocity, the tangential stress, and the temperature perturbation vanish on it .
Z =
—
1;
~
=
0,
T~
2 = 0, 0= 0 (10) Free boundary conditions for the lower surface simplify the analysis permitting us to obtain an
94
analytical solution without significantly changing the results. In the next sections only perturbations periodic in x will be considered
where u~,(0, 0) is the initial displacement of the upper surface (z = 0) at t = 0. Since the displacement u~(0,t) is related to the velocity by v~= du~/dt,the relation of their Laplace transforms is
b1 =b1(z,
v(0, s) =su’(O, s) —u~(0,0)
t)
cos kx,
b2=b2(z,
t)
sin kx (11)
where k is the dimensionless wavenumber, b1 is one of variables v~,0, p, ~ ;~ b2 is one of v~, ~
3. The estimation of rheological parameters by postglacial rebound data Representing the basic eqs. 1—5 in terms of their Laplace transforms and eliminating all the variables except for the vertical velocity, we find [(i
—
~~)sH*(s)L4
—
Rak2I ~*(z, s)
2O(z, 0) (12) Rak where the asterisk denotes the Laplace transform, s is the Laplace transform variable, 0( z, 0) the initial temperature perturbation at t = 0, L the differential operator L2 = D2 k2, D = dz * (s) relates to the creep kernel Laplace trans=
—
form by 11*(s)
=
1/s[1
(
K*(s)]
+
(13)
—
(~~)
where 11(t) is the memory function. Thus, 1~I*(s) in eq. 12 is the Laplace transform of the memory function. The boundary conditions given by eqs. 9 and 10 become 2 + k2) ~ = 0 z = 0; (D [sH*(D3 3k2D) ~pk2] v* = ~pk2u~(0, 0) —
—
—
Rak2]v~*= Rak2 u
[sH*L
2(0, 0)
(15)
Such a form of the initial conditions corresponds to the known problem of harmonic ice loading. Solving the problem given by eqs. 12, 15, 16 and 19 we find v~*(z,s) and then, using eq. 17, we obtain uo2kB(k)q.,_hH*(s) (20) U’ (0, s) = (1 + 2kB k) ~ ‘s11 * (s))
(
—
where B(k) = (2k + sinh 2k)/(cosh 2k 1) (21) Since the parameter Ra/i~= aAd 3• 10-2 is small, the temperature gradient present in the layer can be neglected in the postglacial rebound —
v~=D
Theis final step of inverting transform extremely difficult. Wethe shallLaplace seek only an asymptotic solution at large t using the known theorem of asymptotic expansion of the pre-image (see, e.g., von Doetsch, 1967). The Laplace transform given by eqs. 20 and 18 has the singularity s = 0. We expand the right-hand side of eq. 20 in the vicinity of s = 0 in a power series of s. According to the theorem, to find an asymptotic solution it is enough to invert the first term of the expansion. So we obtaindisplacement the asymptotic dependence of the upper surface on time u~(0,t)
=
u
02kB(k)~jf~’$t~/(yF(y)I’(1—y)) (22) For a Newtonian fluid the Laplace transform of the dimensionless memory function reduces(23) to 2
2v~*=D4v z= —1;
For the Andrade medium, as follows from eqs. 8 and 13 K*(s)=YF(y)/$sY, (18) fl*(~)=$s~’/(yF(y) +/3s~) where F(y) is the gamma-function of y. Now consider initial perturbations of the form 0(z,0)~0, u~(0,0)=uo (19)
problem.
Note that the constitutive eq. 5 can be replaced by = 2f11 t t 1) de~~ (14)
4
(17)
2~~=0
(16)
H*(s)=vpic/,5d
95
where p is the kinematic viscosity. Substituting eq. 23 into eq. 20 and inverting the Laplace transform, we obtain the known solution of the post-
gives ~~(z)
glacial rebound problem for the Newtonian rheological model of the mantle u~(0,t)=u0exp(—t/T) V
T =
K
2kB(k)~ ~
(24) (25)
where T is the relaxation time nondimensionalized 2/K. Unlike eq. 22 obtained for large t, eq. 24 is valid for arbitrary t. As follows on the time scale d from eq. 24, the relaxation time is the time at
=
2e~~(z)f tH(
—
‘~)exp[iw(t1
—
t)]
=
2ê1j(z)f
0
2)e”’12 dt2
(28)
where d,1(z) is the complex amplitude of the deviatoric strain rate tensor. As follows from eq. 28, the solution of the basic equations may have the form (27) only at large t when we can take this case = ~ as the upper limit of the latter integral. In T,~(z) =
2F(w)è 11(z)
which the upper surface displacement decreases by a factor e. Although in the framework of the Andrade model the dependence of displacement on time is not exponential, we define the relaxation time in this case by the equation u0/u~(0, T) = e too. Then, as follows from eq. 22 2kB(k)/1’$e/(’fI’(’f)F(l —v)) (26) The laboratory studies of rock deformation show that y= 1/3 y.is Equation a typical value of the us rheological parameter 26 permits to estimate the rheological parameter /3 at the fixed value of -y with the help of the postglacial rebound data. When the dimensions of ice load are nearly 100 km (k = ir) the relaxation time is 6 X 10~s (T 6~106). Substituting ‘f = 1/3, k = ~r, T 6~106 into eq. 26 we find that /3 5. i0’. For the same values of k and T eq. 25 gives v 1020 cm2 s~. This known estimate may be considered as the layer effective viscosity for postglacial rebound processes but cannot be used for investigation of the other geodynamical processes if the true rheology of the mantle is not Newtothan. =
4. Thermoconvective instability
(29)
where
F(c~)=
j
l1(t)e”
dt = 11*(iw)
As follows from eq. 29, the complex function F( ca) is an analogue of the viscosity. Hence, F( ~,) is naturally called the ‘complex viscosity’. For the Andrade medium the complex viscosity is
1/(’fr(y)
+ /3(iw~) (30) F(w) = Substituting eqs. 11 and 27 into the basic equations and boundary conditions given by eqs. 9 and 10 and eliminating the complex amplitudes of all physical variables except for the vertical velocity amplitude, we obtain a boundary-value problem given by the ordinary differential equation
$(1wy
1—
=0
1~—)iwF(ca)L4_Rak2]vz(z)
(31)
The conditions on the upper boundary are ~ = 0; (D2 + k2) v~= 0
[icaF(w)(D~ 3k2D) —
—
i~tik2]Vr
=
0
(32)
[IWF(w)L4_Rak2]vz=O
On the lower boundary z = —1; v~= D2v~= D4v 2= 0
In this section we shall seek a solution for the basic eqs. 1—5 in the form of an oscillation t (27) b1(z, t) = bi(z)e~~~t,b2(z, t) = b2(z)e” where the meaning of b 1, b2 is the same as in eq. 11. The dimensionless frequency a is a complex number. Substitution of eqs. 11 and 27 into eq. 14
dt1
tH(t
(33)
where L and D are the differential operators introduced in the previous section. The problem given by eqs. 31—33 is considered as an eigenvalue problem. The eigenvalue is the critical Rayleigh number which gives the stability criterion. When Ra = Racr the solution of the problem has the form (27) with a real ca. Suppose
96
that the frequency which is sought satisfies the restriction << (34~
Replacing the complex viscosity (30) by (35) means that we neglect the elasticity of the medium for the frequencies considered. Under free boundary conditions the eigenfunctions of the problem are
where v~(z) is given by eq. 36 in which n = 1. In the present section the asymptotic solution of the basic equations in the form of a thermoconvective oscillation has been obtained without use of the Laplace transform. Besides the singularity s = 0 considered in the previous section, the boundary-value problem given by eqs. 12, 15 and 16 has the singularity s = ia associated with the thermoconvective oscillation. The method based on the use of Laplace transforms allows us to find the amplitude of asymptotic oscillation as a function of the initial perturbations 0(z, 0), u.(0, 0). But this advantage of the Laplace transform method is not very important because it can be shown that at any initial perturbations the corre-
v~(z)= C sin nirz
sponding amplitudes differ by a constant factor.
‘
/
Then the upper traction-free boundary behaves as a free one and the eqs. 33 are valid at z = 0 too. Besides, under the restriction (34), since yF(y) 1 for y = 1/3 and ~p<< 1, we can write F(c~)as ,~
\
\5—1 / /3~:w) /yyl’yy))
~35j\
~.
=
~ = 1, 2, 3,
..
(36)
where C is an arbitrary factor. The dispersion relation is given by 2
+ ~2~2) +
F(
)( k2 +
2~ 2)3
—
Racr k2
=
0
(37)
As follows from eq. 37
~)
Ra cr = (k2 +k2 2~2) IReF( F( w) = Im F(“ ~ (k2 + n2ir2) ReF(~~)
2
—
(38) (39)
—
—
-
3 =
F( “m)
‘P
.
hereminimal the is nearlycritical two orders Rayleigh of magnitude number. larger Therefore, than is fornecessary investigation to solve of convective numerically flows the strongly in the layer nonit .
where Re denotes ‘the real part of’, Im is ‘the imaginary part of’. For F(w) given by eq. 35, we find that the critical Rayleigh number reaches its minimum Ram at n = 1, k = km = ~‘r/(1+ y)1/2 For y = 1/3 we obtain km 2.7; cam = ca(km) 30; Ram 144 /3. If Ra = Ram, there is a regime of boundary instability for which only the convec~tive oscillation with wavenumber km = 2.7 (the wavelength (2iT/km)d is 200 km) and frequency = 30 (the period 2TT/Wm d2/K is nearly 108 years) is not attenuated. It is easy to verify that the frequency equal to 30 satisfies the restriction (34) which was a priori assumed. Then the amplitude of vertical displacement of the traction-free upper surface can be found readily from the second equation of set (32) (0)
.
to be 1020 cm2 s’, as shown above. Then the actual Rayleigh number for the layer considered
3
~.
.
5. Finite-amplitude thermoconvective oscillation If a Newtonian fluid is accepted as the mantle’s rheological model, the kinematic viscosity proves
~ 3D ~ k~ —
Z
(40)
linear equations. But for the Andrade rheological model the situation is greatly changed: the study of convection in the layer leads to linear or weakly nonlinear problems. According to the results obtained in the previous sections, y = 1/3, = 5. iO~ may be accepted as the Andrade rheological parameters for the upper mantle boundary layer. At these values of y and /3 the minimal critical Rayleigh number Ram = 144/3 7 iO~. The actual Rayleigh number for the layer was estimated to be Ra iO~. At Ra
.
97
linear stability theory predicts the unlimited increase of oscillation amplitude in time. But for sufficiently large amplitudes the use of linearized equations is not valid. The application of nonlinear theory permits us to obtain an asymptotic solution in the form of a steady thermoconvective of initial conditions. To describe perturbaoscillation the amplitude of which finite is independent tions in the frame of the Boussinesq approximation we must replace the partial derivative of temperature with respect to time in eq. 4 by the total derivative. Thus, an additional nonlinear term ~~(ao/ax) + ~~(80/8z) appears in the heat equation only. Introduce the supercriticality as 1 Ra
Ram Ra~
The next problem is posed by the coefficient of The solvability condition of the problem gives ca~= 0. Then solving the problem we obtain v(2) = ~(2) = 0 2
B 2F(ca k~)2 0) 2 sin 2irz 2 ~ Ra m k m B2F(ca 0) e2~~~~tl + c.c. (45) + 41~2 + 2ica0
~(2)
=
[I
(~2 +
~
I
—
j
Consideration of the problem posed by the coefficient of ~ leads to the solvability condition F( ca0) d F I (icao + ~2 + k~) + I
[
J
—
dca I~=~
—
+ica2F(w0)+— ~,2 1B12[
(41)
F(ca0) + 2ica0
4~2
The supercriticality is assumed to be small. We use it as a small parameter to expand the perturbations. The vertical velocity is sought in the form 2v~2~(z, x, t) = v~’~(z, x, t) + x, t) + (42) where v~(z, x, t) are periodic in x and t. These functions are expressed as Fourier series. The period in x is fixed. It equals 2 iT/km where km ~S found in the linear stability theory. Similar expansions are assumed for the temperature and the ...
other variables. The frequency is expanded in e too = ca 2ca 0 + ca1 + e 2 + (43) ...
When we substitute these expansions into the equations describing finite perturbations and equate the coefficients of “ to zero, we obtain a hierarchy of equations which must be solved successively, The coefficient of e yields a problem which is identical with that of the linear stability theory. Hence, ca0 cam. But now we write the solution of the problem as v(1)
=
sin
rrz
cos
kmX +
c.c.
(44)
where c.c. denotes the complex conjugate, B is an amplitude factor which must be found, ca is given by the expansion (43).
+
ReF(
ca0)J
2~2
—
—0
(46)
Substituting complex viscosity given by eq. 35, y= 1/3, km the = 2.7, and ca 0 = 30 into eq. 46, we find that ca2 —40, B I 36. Thus, the amplitude of vertical velocity oscillation, which is arbitrary in the linear theory, proves to be 36. I
—
Solution (45) obtained from the equations of the second-order approximation in c permits us to find the relation of Nusselt number to supercriticality
(47)
2
Nu=1+a For the Andrade medium at y= 1/3 the coefficient a is given by a 2.4. The solution for the linear stability problem in the form of a thermoconvective wave has been found by Birger (1978, 1980). In the present work the solution in the form of a thermoconvective oscillation, i.e. standing thermoconvective wave is obtained. Unlike the linear theory, solutions haying the form of travelling and standing waves do not relate to each other trivially in the nonlinear theory. If we write the solution of the first-order approximation in wave v~ = Z
B sin 2
srz
exp
in the form of a travelling
i(cat + kmX) +
c.c.
(48)
98
then, using the method developed in the present section, we find that ca2 = 0, I B I 24 (Birger, 1978). —
6. Conclusions Geology testifies that the vertical movements of the Earth’s crust have occurred everywhere during the whole geological history of the Earth. The major characteristic of the vertical movements is a periodic change in sign of movement which caused them to be called oscillatory (Beloussov et al., 1974). On the Earth’s surface there always exists simultaneously areas of movement in different (upward and downward) directions. The lowestfrequency vertical movements of the Earth’s crust in the continental platform regions are supposed to have periods of 108_109 years. The regions which rise (or subside) have dimensions of several hundreds of kilometres. The amplitudes of crustal oscillations in the platform regions are about 1 km. The results of the present work permit us to assume that the cause of the lowest-frequency oscillatory crustal movements is thermoconvective oscillation of the upper thermal boundary layer in the mantle, —
References Beloussov, V.V., Reisner, G.I., Rudich, E.M. and Sholpo, V.N., 1974. Vertical movements of the Earth’s crust on the continents. Geophys. Surv., 1: 245—273.
Birger, B.I., 1977. Complex viscosity and Boussinesq approximation in analysis of convective instability of upper mantle. DokI. Akad. Nauk SSSR, 235: 794—797 (in Rus. sian). Birger, B.I., 1978. The finite amplitude thermoconvective waves in the Earth’s mantle. DokI. Akad. Nauk SSSR, 242: 559—562 (in Russian). Birger, B.I., 1980. Thermoconvective waves in the Earth’s mantle. Phys. Earth Planet. Inter., 22: 238—243. Christensen, U., 1984. Convection with pressure- and temperature-dependent non-Newtonian rheology. Geophys. J.R. Astron. Soc., 77: 343—384. Goetze, C. and Brace, W.F., 1972. Laboratory observations of high-temperature rheology of rocks. Tectonophysics, 13: 583—600. Murrell, S.A.F. and Chakravarty, S., 1973. Some new rheological experiments on igneous rocks at temperatures up to 1120°C. Geophys. J.R. Astron. Soc., 34: 211—250. Parmentier, E.M., Turcotte, D.L. and Torrance, K.E., 1976. Studies of finite amplitude non-Newtonian thermal convection with application to convection in the Earth’s mantle. J. Geophys. Res., 81: 1839—1846. Parsons, B. and McKenzie, D.P., 1978. Mantle convection and the thermal structure of the plates. J. Geophys. Res., 83: 4485-4496. Von Doetsch, G., 1967. Anleitung zum praktischen Gebrauch der Laplace-Transformation und der Z-Transformation. Springer, München, 365 pp. Weertman, J., 1978. Creep laws for the mantle of the Earth. Philos. Trans. R. Soc. London A, 288: 9—26. Yuen, D.A., Peltier, W.R. and Schubert, G., 1981. On the existence of a second scale of convection in the upper mantle. Geophys. J.R. Astron. Soc., 65: 171—190.