RUSSIAN JOURNAL OF EARTH SCIENCES, English Translation, VOL 1, NO 1, DECEMBER 1998
Russian Edition: JULY 1998
A self-consistent 2D model of mantle convection with a floating continent
V. P. Trubitsyn and V. V. Rykov
Schmidt United Institute of Physics of the Earth, Russian Academy of Sciences, B. Gruzinskaya ul. 10, Moscow, 123810 Russia
Abstract. The mantle is modeled by a viscous fluid filling a horizontally elongated 2D region with an aspect ratio of 10:1. A model with Ra = 106 is constructed on a 200x80 mesh. Developed nonsteady-state thermal convection including narrow downwellings and upwellings sets in, with mantle flow velocities ranging from 1 to 10 cm/yr. Then, at a certain moment, a continent floating on the mantle is introduced into the model. The continent is modeled by a thin long plate of a thickness of 0.03 and a length of 2.0 relative units with respect to the mantle thickness, which corresponds to dimensional values of 90 and 6000 km, respectively. To demonstrate mantle heating beneath continent, the latter is positioned at the coldest place of the mantle where downgoing flows dominate at the moment chosen. The evolution of the mantle-continent system is found from numerical solution of equations governing the momentum, mass, and heat transfer in viscous fluid and rigid continent. The problem is rigorously formulated, a self-consistent method is given for the solution of coupled integrodifferential equations, and a technique of their numerical implementation is described. The continent remains virtually immobile during a long time (about 500 Ma), but the mantle flow pattern dramatically changes, which results in suppression of cold mantle downwellings under the continent and their gradual replacement by hot upwellings. Afterwards the viscous drag of mantle flows begins to move the continent at a variable velocity averaging about 1 cm/yr. The mantle flow pattern and continent velocity constantly changes under the action of mechanical coupling and thermal interaction between the mantle and moving continent. After a time of about 1.5 X 109, when the continent has traveled over a distance of about 15000 km, it arrives at a place where several cold mantle downwellings concentrate. Then the continent velocity sharply decreases, and the continent continues its motion in either primary or opposite direction, depending on the general mantle flow pattern. The results of the numerical experiment can be used for the analysis of mechanism responsible for the motion of Eurasia-type continents, origination and ascent of plumes, and geodynamic processes in the subcontinental mantle.
1. Introduction
Presently, various numerical experiments have been used for the study of thermal convection in the Earths mantle. The role of variable, temperature-dependent viscosity, pressure and temperature, as well as the effect of phase transformations on the mantle flow pattern, have been examined in terms of 2D and 3D models
©1998 Russian Journal of Earth Sciences.
Paper No. TJE98001.
Online version of this paper was published on July 20 1998. URL: http://eos.wdcb.rssi.ru/tjes/TJE98001 /TJE98001.htm
in Cartesian and spherical coordinates (see e.g. [Christensen, 1983, 1984; Christensen and Yuen, 1984; Glatz-maier, 1988; Machetel and Weber, 1991; Lenardic and Kaula, 1994; Solheim and Peltier, 1994; Tackley et al., 1994; Parmentier et al., 1994]). Numerical modeling has provided an insight into many basic regularities of global geodynamics. In the nearest decade, a numerical modeling problem of vital importance is the construction of self-consistent models in which oceanic high-viscosity lithosphere breaks into separate plates.
Another important problem is the exploration of the effect of continents on the mantle convection structure. Continents cover more than a quarter of the Earth’s surface. The mantle heat flow in continental areas is
about 30 mW/m2, which is one-third as much as the heat flow in oceanic areas, averaging 90 mW/m2. This is explained by the fact that continents are thermal insulators of the mantle, because only the conductive mechanism of heat transfer is operative within continents. Various estimates of the Nusselt number, characterizing the whole-mantle convective heat loss, give a value of about 20. The oceanic lithosphere, involved in the convective circulation of matter, virtually remains within the convective thermal boundary layer and, despite its high viscosity, cannot serve as an effective thermal screen for the mantle. Thus, even immobile continents should considerably affect the structure of mantle convection.
Based on a simplified model, Christensen, [1983] compared the mantle flow patterns beneath continents and oceans. Trubitsyn and Fradkov, [1985] showed that thermal convection in the upper mantle is suppressed, resulting in a threefold decrease in the continental heat flow. As was shown in later works [Trubitsyn et al., 1993a, 1993b, 1994; Trubitsyn and Bobrov, 1993; 1996, 1997; Lowman and Jarvis, 1993, 1995, 1996; Bobrov and Trubitsyn, 1996; Nakakuki et al., 1997], an immobile continent initially suppresses the underlying mantle convection and broadens the convective cell; afterwards, when the subcontinental mantle had been heated during the following few hundreds of millions of years, a hot upwelling mantle flow develops beneath the continent.
Since continents are not fixed but float on the mantle, their effect on the mantle convection structure is even greater. In early studies, the effect of continental drift was included in a simplified form, as an effective boundary condition. In continental areas, the free upper boundary condition was replaced by fixed values of horizontal velocity [Lux et al., 1979]. In a similar manner, Dom M.-P., [1997] simulated the effect of moving continents by a time-dependent velocity specified at the upper boundary.
Gurms, [1988] presented the results of a consistent numerical 2D model of mantle convection with free floating continents. For the first time, a numerical 3D model of mantle convection with two floating continents was constructed by Trubitsyn and Rykov, [1995] and Rykov and Trubitsyn, 1996a, 1996b]. This model reconstructs general regularities of the formation and breakup of Pangea [Trubitsyn and Rykov, 1995]. Structures similar to the Atlantic and Pacific oceans form upon the Pangea breakup. In terms of this model, a very steep subduc-tion zone (of the Kuril-Kamchatka type) develops at one of the Pacific margin, and a very gently dipping subduction zone (of the South American type) develops at the other. Assuming an alternative initial position of continents, Rykov and Trubitsyn, [1996b] obtained a
structure consisting of two coupled continents (similar to North and South Americas), which developed upon the breakup of a supercontinent.
Trubitsyn and Rykov, [1997] developed a self-consistent numerical model of low-Ra convection in the upper mantle having a variable viscosity and interacting with a moving continent modeled by a thick rigid plate. The cases with a free plate floating over the mantle and moving at a prescribed fixed velocity were considered. The continent approaching a downgoing cold mantle flow was found to deflect this flow and form structures similar to inclined subduction zones. The dip angle of the downgoing mantle flow decreases with the velocity of the approaching continent.
This work describes in detail the mathematical formulation of the problem and the method of its solution. A Ra = 10® convection model close to the real Earth is constructed on a 200x80 mesh, including a thin continent of thickness d = 90 km and horizontal size I = 6000 km. The long-term evolution of the mantle-continent system is calculated. Comparison of the evolution of unsteady-state convection in the mantle with and without a continent shows that the moving continent dramatically changes the mantle convection structure.
2. Model
The mantle is modeled by an incompressible constant-viscosity fluid occupying the volume of an elongated rectangular 2D layer of a thickness D and length L, with the aspect ratio L : D = 10 : 1. The upper boundary is assumed to be free and other boundaries are impermeable, with a slip condition. The lower and upper boundaries have fixed temperatures T = To and T = 0, respectively; i.e. the layer is heated from below. The side walls are thermally insulated. The coordinate origin coincides with the left-hand lower corner, the z axis is directed upward, and the x axis is directed to the right.
Continents are represented as light rigid rectangular thick plates that are heat conductive, float in the mantle, and have the length I and thickness d + do, where d is the depth to which the continent is immersed in the mantle and do is the continent elevation above the mantle. The viscous coupling with mantle flows drives the continent toward downgoing mantle flows. Because the heat transfer mechanism within continents is purely conductive, the continents produce a heat insulating effect which prevents heat to escape from the mantle. Also, the continents affect mechanically the mantle, because the no-slip condition tends to attenuate velocity contrasts in mantle flows adjacent to the lateral and lower boundaries of the continents.
3. The System of Equations of Mantle Convection with Floating Continents
3.1. Equations of thermal convection
Thermal convection in a viscous mantle is described by the distributions of convective vector velocities Vl(x, y, z), temperature T(x,y,z), and pressure p(x,y,z). These unknown functions are found from the solution of three equations of momentum, heat, and mass transfer
prdVi
dt
dp
dx;
dxj
■ pgSi;
dp
dt
d^ = d(k-^)dx-dt 1 dxi ’ %
d(ViP)
dxi
= 0, i = 1,2,3,
(1)
(2)
(3)
viscous stresses
_ dVi dVj % + d^}
(4)
dp
dx
,d2Vx
dp
d~z + dT ~dt
dz2 dxdz d2Vz d2K
dx2
= 0,
dx2 dzdx
vxdT , vzdr
dx
dz
dVx dVz
dx dz
where Ra is Rayleigh number,
,52v;
dz2 _ d2T dx2
= 0,
+ RaT = 0,
d2T dz2 ’
Ra =
apogTpD
kp
(5)
(6)
(7)
(8)
(9)
3.2. Equations of motion of a continent
In a 2D model, velocities u = (ux,uy) at all points of a rigid continent floating in the mantle (along its horizontal surface) have the same values and are equal to the velocity of its center of gravity:
,(x, y) = M0, uy(x, y) = 0
(10)
The continents move under the action of viscous forces of mantle flows. The Euler equations of horizontal movement of a rigid continent have the form
Mduo
dt
J f {-pSxj + Sxj)nj
(11)
where p is density, g is gravity, T is temperature measured from the adiabatic distribution, k is diffusivity, Sij is the Kronecker delta, and Sa is the deviator tensor of
where df is the elementary area of the surface, rij is the unit vector of outward normal to the area, and M is the dimensionless mass of the continent, referred to the unit length of the y axis in the 2D model considered. The integral is taken over the whole surface of the continent portion immersed in the mantle, which includes its base (z = 1 — d) and side walls (x = x\ and x = x\ + /):
Here p is viscosity. The relative values of inertial terms in the left-hand side of mantle momentum transfer equation (1) are of the order of kp/p ~ 10-23 with respect to the right-hand terms. Therefore, these inertial terms may be neglected.
We put p = po(l — a) in the buoyancy term of equation (1) and p = po in other terms of equations (l)-(2). Hereinafter the pressure is measured from its hydrostatic distribution po(z) determined from the condition Vp = —pog. Also, we introduce dimensionless variables taking, as measurements units, the mantle thickness D for length, k/D for velocity, D2/k for time, To for temperature, p for viscosity, and pk/D2 for pressure and stresses.
Using these variables, 2D model equations of convection (2)-(4) take the form
Mduo
dt
/ [p{x = Xi) - Sxx(x = Xi)]dz-
Jl-d
/ [p{x = x2) - Sxx(x = X2)\dz- (12)
Jl-d
- Sxz(z = 1 - d)dx
J X i
where d and I are thickness and length of the immersed portion of continent; x\(t) and X2(t) = x\(t) + I are instantaneous coordinates of the left- and right-hand walls of the continent, which obey the conditions
dx i dt
= u 0,
(13)
The equation governing the temperature distribution Tc within the rigid continent in the initial fixed system of coordinates is reduced to the equation of heat conductivity with advective heat transfer controlled by the continent velocity Mo along x axis,
dT dT
~dt +U°~dx = d(kdTc/dxi)dxi (i4)
Similar to the mantle, the relative value of inertial terms in the left-hand side of the equation of continent motion (12) is of the order of kp/p Pd 10-23 compared to the terms in its right-hand side. Substituting deviator tensor of viscous stresses (4) into (12) and omitting inertial terms, we obtain the equations of continent motion in dimensionless variables
u
2 fi[dVx(x — x\) [p(x = xi) - p(x = x2)j---------
' l — d
dx
dVx(x = x2)
dx
T = 1 at z = 0;
^ = 0,
dx ’
T = Tr
, <9T
and =
dz
'T 'T A 9T ^
T = TC and — = —
ox ox
Vz = 0, and
dVI dx
= 0 when z = 0,
Similar conditions are imposed on the lateral boundaries,
Vx = 0,
(15)
dVx
dz
= 0 when x = 0 and x = 10.
r-X i + i
- fj,[dVx(z = 1 - d)/dz+
J X\
+dVz(z = 1 — d)/dx]dx = 0,
where the coordinate x\ in the integration limits is determined from (13).
3.3. Boundary conditions
Equations of mantle convection (6)—(10), continent motion (15) and (13), and heat transfer within the continent (14) are interconnected through boundary conditions.
Temperature boundary conditions at the lower and lateral boundaries of the whole area studied have the form
(2°)
Condition (19) is also imposed on the free upper surface z=1 outside the continent:
Vz = o,
dVz
dx
= 0
(21)
at x < x\(t) and x > x\(t) + I.
The non-slip condition on the interface between the viscous mantle and moving continent implies that the velocities of mantle flow and continent motion coincide on this boundary:
Vx = uq and Vz = 0
(22)
(16)
at x = 0 and x = 10.
The mantle temperature at the free upper, z = 1, surface vanishes ((T = 0) only in the oceanic area, outside the continent, namely (x < x\(t) and x > x\(t) + /).
At the continental surface immersed in the mantle, the temperature and heat flow are set to be continuous between the mantle and continent; i.e.,
at the base of the continent for x\(t) < x < xi(t) + I and z = 1 — d and at its wall sections for 1 — d < z < 1, x = x\ and x = x2.
Conditions (22) simplify the equation of continent motion. Since the impermeability condition is satisfied all along the base of the continent irrespective of x, we also have Vz = 0 at the base. Consequently, the last term under the integral sign in (15) vanishes. Likewise, since the no-slip condition Vz = 0 is satisfied all over the walls of the continent irrespective of z, we have dVz/dz = 0. Then, incompressibility condition (8) implies that dVx/dx = 0 on the walls. As a result, the equations of continent motion assume a simple form
at the base of the continent, z = 1 — d and x\(t) < x < xi(t) + I, and
at its lateral surface sections, 1 — d < z < 1 and x = x\ or x = x2.
Zero temperature is set on the upper surface of the continent, i.e., Tc = 0 at x\(t) < x < x\(t) + I and z = 1 + do, where do is the elevation of the continent above mantle surface (commonly, do <C d).
As noted above, impermeability and slip conditions are imposed on mantle flows on the lower and lateral boundaries of the entire calculated region. Therefore, the lower boundary condition is = 0, Sxz = 0 at z = 0. Because it is valid at all x, we have Vz/x = 0 at z = 0. In view of (4), we obtain
fl
/ \p{x = x\) — p(x = x2 )\dz—
Jl — d
fXl+l rdvx(z = i-d)1J n
- nl---------—----------]dx = 0,
dx i
= u0, where xi(t = 0) = xio
(23)
(24)
(19)
Finally, the resulting mathematical problem may be stated as follows. In all, there are seven unknown functions. These are four position- and time-dependent functions for the mantle convection: two velocity components of mantle flows Vx(x,z,t) and Vz(x,z,t), temperature distribution T(x,z,t), and pressure distribution p(x,z,t), and three functions for the continent: temperature distribution within the continent Tc(x, z,t), instantaneous translational velocity of the continent as a whole uo(t), and the coordinate of its left-hand edge x\(t). These functions can be found from a closed system of seven coupled equations: four differential equations of convection (5)—(8), thermal conductivity equation for the continent (14), and the equation of motion
of the rigid continent, which is reduced to the condition imposed on mantle flow velocity derivatives (23) and to relation (24) between the velocity of the continent and its position. The constants of integration of the differential equations are found from boundary conditions
(18H22).
The difference between our problem with a free floating continent and the well-known problem with a fixed continent consists in the fact that the impermeability and no-slip conditions on the upper surface must be satisfied at the place currently occupied by the floating continent whose velocity and position are not a priori known but should be defined by solving the system of coupled differential equations at each time step.
4. Numerical Method
Two basically different methods are applicable to the solution of the system of thermal convection equations including floating continents. In the two-region method [Trubitsyn and Bobrov, 1996, 1997], equations of heat and mass transfer are solved at each time step separately outside and within the continent, and continuity conditions are then set at the interface between continent and mantle. In the one-region region [Trubitsyn and Rykov, 1995; Rykov and Trubitsyn, 1996a, 1996b], computations are conducted in a single coherent domain and explicitly incorporate the jump in material properties between mantle and continents on the surface of the continent.
The numerical algorithm solving the system of thermal convection equations including floating continents can be summarized as follows. Let convective flow velocities and mantle distributions of temperature T(ti) and pressure p{t\), as well as position xi(ti) and velocity uo(ti) of the continent, be known at a time moment t\. We have to find the solution of system (5)-(8), (14), (23) and (24) at the next time moment t2 = ti+At. The new position of the continent *1(^2) at the moment^ is readily found from (24): xi(t2) = xi(ti) + uo(t)l)At. If the continent velocity at this moment uo(t2) were also known, one could solve thermal convection equations (5)-(8) with boundary conditions for temperature (18) and velocities (19)—(22), consistent with the new position of the continent, and find mantle flow velocities Vx{t2) H Vz(t2), temperature T(t2) and pressure p(t2) at the time moment t2. However, the complexity of the problem lies just in the fact that the velocity of the continent uo(t2) is unknown. Its value must comply with mantle flow velocities Vx(t2) and Vz(t2) that obey equation of continent motion (24). Therefore, an iterative technique should be developed for determining this velocity of the continent. In principle, based on continent velocity values chosen in accordance with a specified
enumeration procedure, convective velocity fields and integrals in (23) can be calculated until a value of uo(t2) is found for which the right-hand side of (23) deviates from zero by a quantity e complying with desired accuracy. Since the right-hand side physically represents the force that acts on the continent, we have e > 0 if the chosen velocity of the continent uo(t2) is underestimated and e < 0 if it is overestimated.
A finite difference method was used for solving the system of equations of thermal convection including floating continents [Rykov and Trubitsyn, 1996b]. Numerical solution of temperature transfer equation (7) or (14) employed the flux-corrected transport method of Zalesak, [1979]. Equations for velocities and pressure (5), (6) and (8) were reduced to elliptic equations with variable coefficients (generalized Poisson equations). They were solved with the help of the three-layer triangular method with a conjugate-gradient choice of iteration parameters [Samarskii and Nikolaev, 1978].
5. Numerical Results
5.1. Free nonsteady-state mantle convection
The mantle was modeled by a viscous fluid that occupies an elongated 2D region. Given the mantle thickness D k, 3000 km, its outer circumference 2nR & 40000 km, and its inner circumference 2ir(R — D) K, 21000 km, the region approximating the mantle had an aspect ratio of L : D = 10 : 1. The simplest, Ra=106 model with heating from below was computed on a 200 x 80.
Figures la, 2a, 3a, and 4a present model patterns of the mantle thermal convection developed after the dimensionless time moment f=1.0995. As is known, the whole region, even with no convection, is heated throughout at t > 1. Therefore, we may consider the convection as well-developed by that moment. However, at Ra > 2 x 105 thermal convection is nonstation-ary and quasi-turbulent, because the nonlinear terms VxdT/dx and VzdT/dz in (7) generate new, smaller-degree harmonics at each time step. With increasing Rayleigh number, these nonlinear terms become larger, and smaller-scale time-variable inhomogeneities arise as a result of self-organization.
The dimensionless temperature equal to T = 0 at the upper surface and to To = 1 at the lower one is shown in the figure by variable shades. The magnitude and direction of the mantle flow velocity at each point are indicated by an arrow. Velocity scales are shown in the lower right-hand corners of the figures. Maximum dimensionless values of the velocities are Vmax 1500.
The thin (red in the on-line color figure version) line in the upper part of the figures shows the inferred distribution of the Nusselt number (dimensionless heat flow q = —dT/dz) over the external surface of the region.
Figure 1. The state of developed free thermal convection, adopted as an initial moment at which a continent is introduced into the mantle: (a) free mantle convection; (b) mantle convection with the continent. Velocity vectors are shown by arrows; dimensionless temperature is quantified by shades ranging from black (T = 0) to light (T = 1). The upper (a) plots are the inferred relative heat flow distribution (thin curve, left-hand axis) and ocean bottom topography (thick curve, right-hand axis). The upper (b) plot is the relative heat flow distribution.
The mean Nusselt. number is Nu Pd 17. The amplitude of relative variations in the heat flow is about 100%).
The thick (green in the on-line color figure version) line shows the inferred topography, i.e., deformations of the upper free surface arising under the action of viscous convection stresses. The scale of dimensionless relief elevations multiplied by a factor of 10-3 is shown on the right-hand axis. The elevations were estimated from the formula h = Szz/pg. Taking into account the stress unit defined in Section 3, the relief elevation unit is pk/(pgD2). As seen from Figures la, 2a, 3a, and 4a, mean amplitudes of the dimensionless elevations are of an order of 25 x 103.
Setting, for the Earth, D Pd 3 • 10b m, a Pd 3 • 10-5 K-\ AT Pd 2 • 103 K, k Pd 10-6 mV1, p Pd 1022 Pa s, p Pd 5 • 103 kg/m3, and g Pd 10 m s-2, the time and velocity units can be found. The velocity unit is k/D Pd 10-5 m/yr, and the time unit is r = D2/k Pd 3 x 1011 years. According to (9), the Rayleigh number, which characterizes the intensity of mantle thermal convection, is Ra Pd 8 • 106.
Since the above parameters of the Earth are not accurate enough, we used an approximate Rayleigh num-
ber Ra = 10b. Since mantle flow velocities obey the proportionality relation V' ~ Ra2/3 [Turcotte and Schubert, 1982], the velocities and times can easily obtained for somewhat different Rayleigh numbers. With the Earth parameters specified above, the inferred results and observed data may be compared if the velocity values are multiplied by 82/3 = 4 and the time is divided by 4. As a result, maximum mantle flows have Vmax ^ 1500 x 4 x 10-5 m/yr Pd 6 cm/yr.
Figures la, 2a, 3a, and 4a illustrate the inferred evolution of mantle convection over a large time interval. The dimensionless times are indicated in the left-hand parts of the figures. As noted above, the unit time is t Pd 3 x 1011 years. In order to calculate the time interval between the convection states shown in Figures, the dimensionless times indicated in the figures should be multiplied by 3 x 1011 years and divide by 4 (due to Rayleigh number rescaling); i.e., the dimensionless times should be multiplied by 75 x 109 years. Then, the overall evolution time of mantle convection shown in Figures is (1.1495 — 1.0995) x75x 109 years Pd 3.6x 109 years.
The relief elevation unit is pk / (pgD2) Pd 0.02 m so
Figure 2. Same as in Figure 1, the time moment i = 109 years. At the unsteady-state convection, the mantle flow pattern and temperature distribution are not constant in time. Due to the heat insulating effect, the subcontinental mantle is heated, and a system of hot plumes arises under the continent.
that dimensionless amplitudes of inferred relief elevations of 2.5 x 104 correspond to 0.5 km.
5.2. Evolution of the mantle convection with a floating continent
Figures lb, 2b, 3b, and 4b present the inferred mantle convection patterns including thermal and mechanical interaction effects with a moving continent. For comparison, Figures la, 2a, 3a, 4a and Figures lb, 2b, 3b, 4b show the patterns calculated at the same time moments. A continent represented by a rigid plate of the thickness d = 0.03 D = 90 km and length I = '2 D = 6000 km is introduced at the time moment t = 1.0995. In Figures lb, 2b, 3b, and 4b, the continent is shown by gray. It is intentionally chosen to initially occupy the coldest place in the mantle where outgoing mantle flows are most abundant. The figures show that, due to the heat-insulating effect, the mantle under the continent is gradually heated. By the timet = 1.1151 Pd 1.2x 109 years, a well-defined system of hot mantle upwellings has formed under the continent, and the continent starts drifting to the right under the action of mantle flow viscous drag force. Since the continent at each moment interacts with mantle flows, the mantle convection structure changes as the continent moves. The dimensionless continental drift velocity first increases from 100 to 400 and then
decreases. According to the definitions made above, dimensional velocities range from 0.4 to 1.6 cm/yr. As seen from Figures lb, 2b, 3b, and 4b, after two billion years, when the continent had traveled a distance of about 5 D Pd 15000 km, it nearly stops.
In the on-line version of this paper, Figures la, 2a, 3a, 4a and lb, 2b, 3b, 4b are complimented by movies (Figures 5 and 6, respectively), each including 30 frames.
6. Conclusion
The evolution of free unstable-state mantle convection and evolution of mantle convection coupled with a free floating continent were calculated. The comparison of concurrent distributions of temperature and mantle flow velocities derived from these two models has demonstrated that floating continents dramatically change the structure and evolution of mantle convection. Detailed evolution results presented in this paper are useful for investigations into the driving mechanism of Eurasia-type continent motion, plume origination and upwelling, and thermal regime of subcontinental mantle.
Figure 3. Same as in Figure 1, the time moment t = 1.5 x 109 years. Due to self-organization of nonlinear processes, the free floating continent starts drifting. It traveled a dimensionless distance of 1.9 over 0.5 x 109 years, which yields a velocity of 1 cm/yr.
Acknowledgements. This work was supported by the
International Center of Research and Technology, Grant
No. 415-96 and the Russian Foundation for Basic Research,
Project No. 96-05-66069.
References
Bobrov A. M., Trubitsyn V. P., Physics of the Earth, N. 7,
1995, p. 5-13 (in Russian).
Bobrov. A. M., Trubitsyn V. P., Times of rebuilding of mantle flows beneath continent, Izvestiya, Physics of the Solid Earth, 31, 1996, pp. 551-559.
Christensen U., A numerical model of coupled subcontinental and oceanic convection, Tectonophysics, 95, 1983, pp.
1 - 23.
Christensen U. R., Convection will pressure- and temperature-dependent non-Newtonian rheology, Geophys. J. Astr. Soc., 77, 1984, pp. 343-384.
Christensen U. R. and Yuen D. A., The interaction of a subducting lithospheric slab with a chemical or phase boundary, .7. Geophys. Res., 89, 1984, pp. 4389-4402.
Doin M.-P., Fleitout L. and Christensen U., Mantle convection and stability of depleted and undepleted continental lithosphere, .7. Geophys. Res., 102, 1997, pp. 2771-2787.
Glatzmaier G. A., Numerical simulation of mantle convection: Time-dependent, three-dimensional, compressible,
spherical shell, Geophys. astrophys. Fluid Dyn., 43, 1988, pp. 223-264.
Guillou, L. and Jaupart G, On the effect of continents on mantle convection, .7. Geophys. Res., 100, 1995, pp. 24217-24238.
Gurnis M., Large-scale mantle convection and aggregation and dispersal of supercontinents, Nature, 332, 1988, pp. 696-699.
Gurnis M. and Zhong S., Generation of long wavelengh het-erogeneitiey in the mantle dynamics interaction between plates and convection, Geophys. Res. Lett., 18, 1991, pp. 581-584.
Lenardic A. and Kaula W. M., Tectonic plates, D thermal structure, and the nature of mantle plumes, .7. Geophys. Res., 99, 1994, pp. 15697-15708.
Lowman J. P. and Jarvis J. T., Mantle convection flow reversals due to continental collisions, Geophys. Res. Lett., 20, 1993, pp. 2091-2094.
Lowman J. P. and Jarvis J. T., Mantle convection models of continental collision and breakup incorporating finite thickness plates, Phys. Earth Planet. Inter., 88, 1995, pp. 53068.
Lowman J. P. and Jarvis J. T., Continental collisions in wide aspect ratio and high Rayleigh number two-dimensional mantle convection models, .7. Geophys. Res., 101, Bll,
1996, pp. 25485-25497.
Lux R. A. Davies G. F. and Thomas J. H., Moving litho-
Figure 4. Same as in Figure 1, the time moment t = 2.0 x 109 years. The no-continent, convection pattern changes but repeats itself quasi-periodically. The continent is still moving. It traveled a distance of 5 over t = 109 years, which yields 3 cm/yr.
spheric plates and mantle convection, Geophys. J. R. As-tron. Soc., 57, 1979, pp. 209-228.
Machetel P. and Weber P., Intermittent layered convection in amodel mantle with an endothermic phase change at 670 km, Nature, 350, 1991, pp. 55-57.
Nakakuki T. Yuen D. A. and Honda S., The interaction of plumes with the transition zone under continents and oceans, Earth Planet. Sci. Lett., F{6, 1997, pp. 379-391.
Parmentier E. M. Sotin C. and Travis B. J., Turbulent 3-D thermal convection in an infinite Prandl number, volu-metrically heated fluid; Implication for mantle dynamics, Geophys. J. Int., 116, 1994, pp. 241-254.
Rykov V. V., Trubitsyn V. P., Computational seismology, v. 26, Geodynamics and Earthquake Prediction, 1994, p. 94-102 (in Russian).
Rykov V. V., Trubitsyn V. P., Computational seismology, v. 27, Theoretical Problems of Geodynamics and Seismology, 1994, p. 21-41 (in Russian).
Rykov V. V. and Trubitsyn V. P., Numerical technique for calculation of three-dimensional mantle convection and tectonics of continental plates, in Computational Seismology and Geodynamics, v. 3ed. by D. K. Chowdhury. Am. Geophys. Un., Washington D.C. 1996a. pp. 17-22. (English translation of Vychisliteljnaya seismologiya, Nauka, M., No. 26, 1994, pp. 94-102)
Rykov V. V and Trubitsyn V. P., 3-D model of mantle convection incorporating moving continents, in Compu-
tational Seismology and Geodynamics, v. 3 ed. by D. K. Chowdhury. Am. Geophys. Un., Washington D.C., 1996b, pp. 23-32. (English translation of Vychisliteljnaya seismologiya, Nauka, M., No. 27, 1994, pp. 21-41)
Samarsky A. A., Nikolayev E. C., Methods of Solving Residual Equations, Nauka, M., 1978, 591 pp. (in Russian).
Samarskii A. A. and Nikolaev E. S., Method of solving finite-difference equations (in Russian), Nauka, M., 1978, 591 p.
Solheim L. P. and Peltier W. R., Phase boundary deflections at 660-km depth and episodically layered isochemical convection in the mantle, .7. Geophys. Res., 99, 1994, pp. 15861-15875.
Tackley P. J., Stevenson D. J., Glatzmaier G. A., Schubert G., Effect of multiple phase transitions in three dimension spherical model of convection in Earth’s mantle, .7. Geophys. Res., 99, 1994, pp. 15877-15901.
Trubitsyn V. P., Belavina Yu. F., Rykov V. V., Physics of the Earth, N. 11, 1993, p. 3-13 (in Russian).
Trubitsyn V. V., Belavina Yu. F., Rykov V. V., Thermal and mechanical ineraction of mantle and continental lithosphere, Izvestiya, Physics of the solid Earth, 29, 1993, pp. 933-945. (English translation of Fizika Zemli, No. 11, 1993, pp. 3-13. Published by American Geophysical Union and Geological Society of America).
Trubitsyn V. P., Belavina Yu. F., Rykov V. V., Physics of the Earth, N. 7/8, 1994, p. 5-17 (in Russian).
Trubitsyn V. V., Belavina Yu. F., Rykov V. V., Ther-
mal mantle convection in a varying viscosity mantle with a finite-sized continental plate, Izvestiya, Physics of the solid Earth, 30, 1995, pp. 587-599. (English translation of Fizika Zemli, No. 7/8, 1994, pp. 5-17. Published by American Geophysical Union and Geological Society of America).
Trubitsyn V. P., Bobrov A. M., Physics of the Earth, N. 9, 1993, p. 27-37 (in Russian).
Trubitsyn V. P., Bobrov A. M., Kubyshkin V. V. Physics of the Earth, N. 5, 1993, p. 3-11 (in Russian).
Trubitsyn V. P., Bobrov A. M., and Kubyshkin V. V., Influence of continental lithosphere on structure of mantle thermal convection, Izvestiya, Physics of the solid Earth, 29, 1993, pp. 377-385. (English translation of Fizika Zemli, No. 9, 1993, pp. 27-37. Published by American Geophysical Union and Geological Society of America).
Trubitsyn V. P., Bobrov A. M., Computational seismology, v. 21, Theoretical Problems of Geodynamics and Seismology, 1994, p. 3-20 (in Russian).
Trubitsyn V. P., Bobrov A. M., Computational seismology, v. 28, Modern Problems of Seismicity and Earth Dynamics, 1996, p. 22-31 (in Russian).
Trubitsyn V. P. and Bobrov A. M., Structure evolution of mantle convection after breakup of supercontinent, Izvestiya, Physics of the solid Earth, 29, 1994, pp. 768-778. (English translation of Fizika Zemli, No. 9, 1993, pp. 27-37. Published by American Geophysical Union and Geological Society of America).
Trubitsyn V. P. and Bobrov A. M., Thermal and mechanical interaction of continents with the mantle, in Computational Seismology and Geodynamics, v. 3, ed. by D. K. Chowdhury. Am. Geophys. Un., Washington D.C,
1996, pp. 33-41. (English translation of Vychisliteljnaya seismologiya, Nauka, M., No. 27, 1994, pp. 3-20).
Trubitsyn V. P and A. M. Bobrov A. M., Structure of mantle convection beneath stationary continents, in
Computational Seismology and Geodynamics, v. Jt, ed. by D. K. Chowdhury . Am. Geophys. Un., Washington D.C.,
1997, pp. 42-53. (English translation of Vychisliteljnaya seismologiya, Nauka, M., No. 28, 1996, pp. 22-31).
Trubitsyn V. P., Fradkov A. C., Physics of the Earth, N. 7, 1985, p. 3-14 (in Russian).
Trubitsyn V. P. and Fradkov A. S., Convection under con-tinnts and oceans, Izvestiya, Physics of the solid Earth, 21 , No. 7, 1985, pp. 491-498. (English translation of Fizika Zemli, No. 7, 1985, pp. 3-13. Published by American Geophysical Union and Geological Society of America).
Trubitsyn V. P. and Rykov V. V., A 3-D numerical model of the Wilson cycle, I. Geodynamics, 20, 1995, pp. 63-75.
Trubitsyn V. P., Rykov V. V., Physics of the Earth, N. 6,
1997, p. 1-12 (in Russian).
Trubitsyn V. P. and Rykov V. V., Mechanism of formation of an inclined subduction zone, Izvestiya, Physics of the Solid Earth, 33, No. 6, 1997, pp. 427-437. (English translation of Fizika Zemli, No. 6, Interperiodica Publishing, 1997, pp. 1-12).
Turcotte D. L. and Schubert G., Geodynamics: Applications of Continuum Physics to Geological Problems, John Wiley, New York, 1982, 449 pp.
Zalesak S. T., Fully multidimensional flux-correced transport algorithms for fluids, I. Comp. Phys., 31, 1979, pp. 335-361.
Zhong S. and Gurnis M., Dynamic feedback between a continentlike raft and thermal convection, I. Geophys Res., 98, 1993, pp. 12219-12232.
Zhong Sh. and Gurnis M., Role of plates and temperature-dependent viscosity in phase change dynamics, I. Geophys. Res., 99, 1994, pp. 15,903-15,917.
(Received May 15 1998.)