CC BY
27
GEODYNAMICS & TECTONOPHYSICS
PUBLISHED BY THE INSTITUTE OF THE EARTH'S CRUST SIBERIAN BRANCH OF RUSSIAN ACADEMY OF SCIENCES
ISSN 2078-502X
2018 VOLUME 9 ISSUE 1 PAGES 287-307
https://doi.org/10.5800/GT-2018-9-1-0349
Modeling of moving deformable continents by active
tracers: closing and opening of oceans, recirculation of oceanic crust
A. M. Bobrov1, A. A. Baranov1, 2
1 O.Yu. Schmidt Institute of Physics of the Earth of RAS, Moscow, Russia
2 Institute of Earthquake Prediction Theory and Mathematical Geophysics of RAS, Moscow, Russia
Abstract: The evolution of the 'mantle - moving deformable continents' system has been studied by numerical experiments. The continents move self-consistently with the mantle flows of thermo-compositional convection. Our model (two-dimensional mantle convection, non-Newtonian rheology, the presence of deformable continents) demonstrates the main features of global geodynamics: convergence and divergence of continents; appearance and disappearance of subduction zones; backrolling of subduction zones; restructuring of mantle flows; stretching, breakup and divergence of continents; opening and closing of oceans; oceanic crust recirculation in the mantle, and overriding of hot mantle plumes by continents. In our study, the continental crust is modeled by active markers which transfer additional viscosity and buoyancy, while the continental lithosphere is marked only by increased viscosity with neutral buoyancy. The oceanic crust, in its turn, is modeled by active markers that have only an additional buoyancy. The principal result of our modeling is a consistency between the numerical calculations and the bimodal dynamics of the real Earth: the oceanic crust, despite its positive buoyancy near the surface, submerges in subduction zones and sinks deep into the mantle. (Some part of the oceanic crust remains attached to the continental margins for a long time.) In contrast to the oceanic crust, the continental crust does not sink in subduction zones. The continental lithosphere, despite its neutral buoyancy, also remains on the surface due to its viscosity and coupling with the continental crust. It should be noted that when a continent overrides a subduction zone, the subduction zone disappears, and the flows in the mantle are locally reorganized. The effect of basalt-eclogite transition in the oceanic crust on the mantle flow pattern and on the motion of continents has been studied. Our numerical experiments show that the inclusion of this effect in the model considerably alters the pattern of mantle flows and leads to distinct changes in the evolution of continents. Moreover, a new effect arises - bulging of heavy material (eclogitized former oceanic crust) at the core-mantle boundary, where-from it arises with the mantle plumes on the surface of the Earth.
Key words: numerical modeling; thermo-compositional convection; continent; tracer; oceanic crust; recirculation; backrolling
RESEARCH ARTICLE Received: April 25, 2017
Revised: October 10, 2017
Handling Editor: K.Zh. Seminsky Accepted: December 4, 2017
For citation: Bobrov A.M., Baranov A.A., 2018. Modeling of moving deformable continents by active tracers: closing and opening of oceans, recirculation of oceanic crust. Geodynamics & Tectonophysics 9 (1), 287-307. doi:10.5800/GT-2018-9-1-0349.
Движущиеся деформируемые континенты, смоделированные
активными маркерами: закрытие и раскрытие океанов,
рециркуляция океанической коры
А. М. Бобров1, А. А. Баранов1, 2
1 Институт физики Земли им. О. Ю. Шмидта РАН, Москва, Россия
2 Институт теории прогноза землетрясений и математической геофизики РАН, Москва, Россия
Аннотация: В численных экспериментах изучена эволюция системы мантия - движущиеся деформируемые континенты. Континенты движутся самосогласованно с мантийными течениями термокомпозиционной конвекции. Определены основные черты глобальной геодинамики - схождение и сжатие континентов, возникновение и исчезновение зон субдукции, бэкроллинг зон субдукции, перестройка мантийных течений и растяжение континентов с их последующим расхождением, раскрытие и закрытие океанов, а также рециркуляция океанической коры в мантии. Континентальная кора смоделирована активными маркерами, имеющими добавочную вязкость и плавучесть, тогда как континентальная литосфера - маркерами только с повышенной вязкостью, имеющими нейтральную плавучесть. Океаническая кора, в свою очередь, смоделирована активными маркерами, имеющими плавучесть. Принципиальный результат - соответствие численных расчетов реальной двухмодовости динамики Земли: океаническая кора, несмотря на свою положительную плавучесть у поверхности, погружается в зонах субдукции и уходит глубоко в мантию. (Часть океанической коры остается налипшей на окраинах континентов и не тонет, сохраняясь на длительное время.) В отличие от нее, континентальная кора не субдуцирует в зонах субдукции. Континентальная литосфера, несмотря на свою нейтральную плавучесть, за счет вязкости и сцепления с континентальной корой также остается на поверхности. При этом при «наезде» континента на зону субдукции происходит ее исчезновение и локальная перестройка течений в мантии. Изучено влияние перехода базальт - эклогит в океанической коре на структуру мантийных течений и движение континентов. Установлено, что включение в модель эффекта этого перехода существенно меняет картину мантийных течений и положения континентов. Кроме того, появляется новый эффект - возникновение скоплений остатков океанической коры на дне мантии. Вещество океанической коры накапливается на дне мантии неоднородно, образуя несколько скоплений. Поднимаясь вместе с плюмами, вещество океанической коры вновь оказывается на поверхности Земли.
Ключевые слова: численное моделирование; термокомпозиционная конвекция; континент; маркер; океаническая кора; рециркуляция; бэкроллинг
1. Introduction
The effect of continents on mantle convection has been studied for several decades. Already early works [Trubitsyn et al., 1999; Gurnis, 1988] showed that the presence of continents in the model considerably changes the mantle convection pattern. Continents moved by mantle convection themselves change the convection and create a supercontinental cycle, including the convergence of all continents into a supercontinent and its subsequent disintegration accompanied by radical restructuring of all the mantle flows. In the early studies, continents were modeled by a heat-insulating line on the surface without a specified thickness, then by incompressible rectangles.
In [Bobrov, Trubitsyn, 2008], the movements and mixing of the submerged oceanic lithosphere material in the mantle and its subsequent rise to the surface were determined in the successive stages of the calculated supercontinental cycle. Simultaneously, the fields
of viscous maximum shear stresses and the orientations of their axes were calculated from the obtained fields of mantle material flows. The oceanic lithosphere was modeled by passive tracers; the number of tracers in the model was small. In the works mentioned above, the continents were considered to be undeformable.
Later on, various models were used to model the continents. In [Butler, Jarvis, 2004], an axisymmetric spherical model of a mantle with immobile continent modeled by a high-viscosity region and subduction of an oceanic plate near an active continental margin is considered.
In subsequent years, the continents have been modeled by active tracers, which allowed to take into account the deformation and buoyancy of the continents. In [Trubitsyn et al., 2007], subduction of the oceanic and continental crust into the mantle was simulated by the method of active tracers. Their model shows that the oceanic crust (both normal and thickened) submerges in subduction zones together with the oceanic
lithosphere. In its turn, the buoyant continental crust does not sink in the subduction zones. This model has several shortcomings: the square area of computation, low Rayleigh numbers, and Newtonian rheology in the entire design area (softening of the substance at high stresses is not taken into account). Actually, the principal conclusions published in [Trubitsyn et al., 2007] are correct, although based on the estimations from the parameters that differ significantly from the real parameters of the Earth.
In [Bobrov, Baranov, 2011], a continent simulated by tracers with increased viscosity and zero buoyancy was added to the two-dimensional model of mantle convection. It was assumed that the continent moved self-consistently with mantle currents. Stress fields were also calculated at different stages of the evolution of the system.
Magni et al. [2013] considered various modes of lithosphere subduction, such as delamination of the continental lithosphere and slab detachment. The calculations were performed using the well-tested Citcom program for a two-dimensional Cartesian model (limited to the upper mantle) with the aspect ratio of 1 to 5. The continents were modeled by active tracers with positive buoyancy and increased viscosity. A shortcoming of their model is that it does not consider a phase boundary at a depth of 410 km; moreover, a 660 km phase boundary is also absent since the model is limited to the upper mantle.
The spherical model considered in [Yoshida, 2012] includes deformable continents, which are more rigid than low-viscosity zones, i.e. young orogenic belts and suture continental zones. It is shown that the weak low-viscosity continental margins play an important role, partly protecting the platform areas from stretching by the convecting mantle and so ensuring their longevity. For continental areas, the yield stress was assumed to be infinite, so they did not have the property of plastic softening at high stresses. The yield stress of the oceanic lithosphere was assumed equal to 100 MPa. The transport and recycling of the oceanic crust material in the mantle was not considered in [Yoshida, 2012] (the tracers were introduced only for the continental areas).
In our study, considerable attention is paid to the role of a heavy eclogitized material and its effect on the pattern of mantle convection.
Numerous studies have been aimed at computer modeling of thermochemical convection and considered the role of both the heavy component (eclogite) sinking into the mantle in the subduction zones, and the light fraction that results from differentiation and ascends in the mantle from the core-mantle boundary.
The role of the heavy component was numerically investigated, in particular by Lobkovsky and Kotelkin.
In the series of papers, models of various geometry were considered: the two-dimensional Cartesian model [Lobkovskii et al., 2014]; the model of the thin conic ring in the equatorial plane [Lobkovsky, Kotelkin, 2004]; and the spherical model [Lobkovsky, Kotelkin, 2015]. They consecutively introduced more and more new factors into their model and thus refined it to represent more details.
An important element of the model is the presence of phase boundaries at the depths of 410 and 660 km. The latter boundary plays a particularly important role, since it is an endothermic transition, that is, to some degree (depending on the parameter of the problem -the slope of the phase equilibrium curve) prevents the passage of mantle streams. In the computer experiments, Lobkovsky and Kotelkin obtained the phenomenon of intermittent convection: during some stages the convection in the upper mantle and the convection in the lower mantle took place almost separately, which led to substantial cooling of the upper mantle. Due to its weight, the cooled material passed the 660 km barrier and began to sink quite rapidly into the lower mantle, thus causing the convection in the whole mantle. With the appropriate choice of parameters, Lobkovsky and Kotelkin achieved consistency between their results and the duration of the geologically observable cycles. Their model shows that the heavy component (eclogite) plays a significant role: it participates in the descending thermal flows, intensifying them, so the avalanches in the thermochemical convection occur more often than in case of the purely thermal convection.
In the Lobkovsky model, a powerful global descending flow emerges at the stage of global convection and pulls together the material above the flow to make a supercontinent on one side of the planet, and an ascending superplume occurs on the other side. Based on this process termed 'overturn', Lobkovsky et al. proposed a fundamental explanation for the division of the Earth into the oceanic and continental hemispheres.
The above-mentioned studies should not be regarded as providing the final answers to the issues of convection in the Earth's mantle. As the authors themselves point out, some important factors have not yet been taken into account. In the spherical model, the viscosity of the entire volume of the mantle material is assumed to be constant. For the same reason, the model does not take into account the viscosity jump at the depth of the 660-km boundary. Because of this, the flow velocities in the upper and lower mantle in the model should have close values (as in the animation to [Lobkovsky, Kotelkin, 2004]), which seems unlikely. In addition, the inclusion of this viscosity jump in the model should, apparently, influence the process of 'overturn'. Due to the high viscosity of the lower mantle (in the real Earth, the viscosity jump at the 660-km
boundary is 30:1 or more), an avalanche would progress more slowly, and its penetration would take place for a longer time and less intensively.
In contrast to the two-dimensional Cartesian model [Lobkovskii et al., 2014], in the Lobkovsky spherical model, continents do not have an effective viscosity exceeding the viscosity of a homogeneous mantle. The continents are simulated by numerous (~105) tracers transferring the property of buoyancy (but not of continental viscosity and strength). As a consequence, low-viscosity tracers, in our opinion, too easily combine into continents and the supercontinent with their subsequent dispersion; these aggregations are round-shaped. Nevertheless, despite these simplifications, a fundamental conclusion in the qualitative form was obtained about the reasons of the difference in the continental and oceanic hemispheres of the Earth. Further development of this model seems very promising.
In our present work, we consider a two-dimensional Cartesian model and numerically investigate the evolution of mantle flows, temperatures, viscosity fields, and the displacement of the material during the motion of two continental plates on the surface of the convective mantle. Non-Newtonian rheology of the material is taken into account. Two continents of different sizes are superimposed on the mantle flows calculated with the thermal Rayleigh number Ra=2x107. The continental crust and the continental lithosphere, as well as the oceanic crust, are modeled by active tracers. They play an important role - firstly, they show the movement of the material in the course of evolution. For deformable continents, this makes it possible to take into account such effects as compression, stretching of continents, their local thinning and thickening, that is, those processes that take place in reality. The modeling of continents by rigid rectangles or lines does not allow for these processes to be taken into account. Secondly, tracers are defined here as active, i.e. as transferring the properties of a given material: its viscosity (which determines either the preservation of the compactness of this part of the material (at its high viscosity) or its disintegration); and its density in relation to the surrounding material, what determines its buoyancy.
The viscosity of the material in the model depends on the temperature, and also, in the oceanic area to depths of 150 km, on the stresses at a given point. Namely, the effective viscosity drops abruptly when the critical value of the applied stress is reached. This property of the drop in the viscosity of the material when the critical value of the maximum shear stress is reached (non-Newtonian viscosity, quasi-plastic behavior) plays a very important role. Due to such rheology, bending of thick rigid oceanic plates (in the region of high stresses) and their subduction becomes possible. This property is one of the main factors determining the base features of mantle convection.
The model assumes that the viscosity of the continents is higher than the temperature viscosity of the corresponding oceanic area.
It is, of course, evident that the situation in reality is much more complicated. Very interesting results were obtained, in particular, from the regional models with kinematic boundary conditions [Burov et al., 2014].
The advantage of the global dynamic model under consideration is its coverage of the entire mantle as a whole. Unlike regional models, the global model does not require to formally introduce any closely located vertical side walls of the regions and, accordingly, to impose artificial boundary conditions (which sometimes do not even vary with time). Obviously, these factors can influence the derivable results.
The phase boundaries in the mantle at the depths of 410 and 660 km are taken into account in our modeling, since a number of works have shown their essential role. In our model, the density of material abruptly changes at these boundaries. The upper/lower mantle boundary at 660 km also is considered as the position of the jump in mantle viscosity.
As found in a number of studies, for example, in [Tosi, Yuen, 2011; Bobrov, Baranov, 2016], the presence of these boundaries leads to a number of effects, such as the 'two-storied' form of ascending plumes, with some delay of the material immediately underneath the 660-km boundary; immersed slabs 'lying' on this border, which are of frequent occurrence; and, in general, the impeded passage across the 660-km border in both directions. Further, the difference in the viscosities of the upper and lower mantle leads to the fact that the flow velocities in the upper mantle are about 4 times higher than in the lower mantle, and the flows are oriented (mainly) subhorizontally. These significant subhorizontal upper-mantle flows may transfer the upgoing plumes to thousands of kilometers away from the place where they passed the 660-km boundary. In this case, the forms of the transfer of the plumes may be different. If, for example, a large-scale mantle flow comprises a group of plumes, the lateral plumes are strongly affected by the flows and drawn to sides, while the central ones ascends almost vertically. All these effects, of course, also affect other geophysical fields.
2. The equations and the model
We use the Cartesian 2D model. It is assumed that the mantle is heated from the core and by the decay of the radioactive isotopes uniformly distributed in the mantle.
The 2D fluid convection equations for the coordinates х and z in the dimensionless form are as follows.
The continuity equation is
dVx/dx+dVz/dz = 0, (1)
the momentum transfer equations are
-dp/dx+dTxx/dx+dTxz/dz = 0 (2)
and
-dp / dz+dTxz/dx+dr zz/dz = -RaT+Ra cC+RaphAr, (3) the heat transfer equation is
dT/dt+VxdT/dx+Vz dT/dz = d2T/dx2+d2T/dz2+H, (4)
and the transfer equation for the concentration C of the material is
dC/dt+Vx dC/dx+Vz dC/dz =
= (Kd/K)(d2C/dx2+d2C/dz2). (5)
Here, Ra = ap0gATD3/Kn0 is the thermal Rayleigh number; Ra c = ApgD3/Kn 0 is the compositional Rayleigh number; Ap=p 1-P0 is the density difference of considering material and the mantle material; Raph= =5pgD3/(K^ 0) is the phase Rayleigh number that characterizes the effect of the density jump 5p upon phase transition.
r(x, z) is the volumetric fraction of the second phase (which has a higher density). The gamma function is defined through the hyperbolic tangent. Function Ar=r(x,z)-r 0 (where r 0 is the phase distribution function undisturbed by the convection) is Ar = 1 in the region where the phase boundary is elevated relative to undisturbed level ha, and Ar = -1 in the region where the phase boundary is lowered.
The Earth mantle is described in general by the following parameters [Schubert et al., 2001]: the coefficient of thermal expansion a=2x10-5 K-1; gravitational acceleration g=9.8 m/s2; p=4600 kg/m3; the coefficient of thermal diffusion k=10-6 m/s2; the reference kinematic viscosity V0=n0/p=0.5x1018 m/s2; the dimensionless thermometric density of the heat sources H=15 [Lowman et al., 2004]; and the drop in the superadiabatic (potential) temperature between the core/mantle boundary and the surface, AT, which is assumed to be 2500 K in our calculations. The scaling factors for the mantle are D = 2850 km for the length and k/D=1.08x10-3 cm/yr for velocity; t0=D2/K=266 Ga; and c0 = k^0 /D2 = 0.235x103 Pa. These values give the Rayleigh number Ra= =(agATD3)/(KV0)«2x107. We note that Ra is determined here through the temperature difference AT between the lower and upper boundaries. In our model, y 410=1.6 MPa/K, 5p 410/p 0 = =0.07; y660=-1.3 MPa/K and 5p660/p0=0.09 [Fei et al., 2004]. With these parameters, the Rayleigh numbers are Raph410=5p 410 gD3/(Kpv 0>2.8x107; Raph660= =5p 660 gD3/(Kpv 0)«4x107. The values of dimensionless Y are 0.035 for the transition at the depth of 410 km
and -0.025 for the transition at the depth of 660 km. At Ra=2x107, the dimensionless time interval of 1x10-5 corresponds to 2.7 Ma.
In more details the equations are described in [Bobrov, Baranov, 2011,2016].
The deviatoric viscous stress tensor has the following components:
Txx = 2ndVx/dx, (6)
t zz = 2n\dVz/dz, (7)
Txz = n (dVx/dz+dVz/dx), (8)
where n is the dimensionless dynamic viscosity at a given point.
In the general form, the viscosity depending on temperature and strain rate of the material at a given point is described by the following formula [Simon et al., 2009]:
n = A - i/n [e]d-n)/n exp[(Ea + PV)/(nRT)]. (9)
Here, T and P are dimensional temperature and pressure; R is the universal gas constant; A is the coefficient before the exponent (pre-exponential factor); n is the nonlinearity index; e is the second invariant of the strain rate tensor; Ea and V are the activation energy and activation volume, respectively. The term [e]ti-"V" describes the plastic character of the deformations. In the present work we use a model with a simplified temperature dependence of viscosity (Arrhenius's law):
n t = exp(2E/(T+ To )-2E/(Tref+Tbot)), (10)
where the dimensionless parameter E, which is connected with activation energy, determines the viscosity range in our model; T is the dimensionless superadia-batic temperature; Tref is assumed to be 0.5; and the temperature at the mantle bottom is Tbot = 1.
In the present study, E=ln1045=10.36, which corresponds to the activation energy (431 kJ/mole) [Yoshida, 2010]. This value is approximately the activation energy of wet olivine [Karato, Wu, 1993]. In our model, the jump in viscosity at the boundary with the lower mantle is assumed to be 50. Thus, in our model
n t = exp(20.72/(T+1)-20.72/(0.5+1)) for the upper mantle, and
n t = 50 exp(20.72/(T+1)-20.72/(0.5+1)) (11) for the lower mantle.
Considering the oceanic lithosphere (above 150 km), in addition to Arrhenius' law, the plastic character of the deformations should be also taken into account, by the condition
Пmho = min(x/2e, nt), if nt ^30, (12)
where e is strain-invariant at a given point:
e = (0.5 eij ¿tj )1/2, где ¿ij = dVi/dXj+dVj/dXi, (13)
т is the effective yield stress [Yoshida, 2010].
In our study, we assume т to be 50 MPa in dimension form, which falls in the interval that favors the generation of relatively rigid plates [Trubistyn, 2012; Moresi, Solomatov, 1998]. The strength of the continents in our calculations is assumed to be 10 times greater than the strength of the oceanic lithosphere. Thus, we have the continental strength (yield stress) equal to 500 MPa. Note that in [Yoshida, 2010, 2012] the strength of continents is taken to be infinite; [Magni et al., 2013] assumes a value of up to 400 MPa. With the continental yield stress value used in our study, there is no massive sinking of continents into the mantle in the zones of strong descending flows. At lower values of the yield stress of continents, such a run-off in the model can occur.
In our model, the thickness of the continental crust is 40 km, and the continental lithosphere (including the crust) is 150 km thick. The continental crust in the model has a positive buoyancy (the density difference between the continental crust and the mantle is 2800-3200=-400 kg/m3), which corresponds to the Rayleigh number Ra C1=-4x107.
Below we consider two versions of buoyancy of the oceanic crust. In the first version, an oceanic crust with the thickness of 7 km is modeled by active tracers with positive buoyancy (the density difference of the materials of the oceanic crust and mantle is 3000-3200= =-200 kg/m3), which corresponds to the Rayleigh number RaC3=-2x107. The basalt-eclogite phase transition is not taken into account.
In the real Earth, at a pressure exceeding 2.5 GPa [Sobolev et al., 2007], the oceanic crust basalt transforms into eclogite with a density of approximately 3400 kg/m3. The second version of our calculations takes into account this transition, and the buoyancy of the oceanic crust in our model becomes negative from the depth of 80 km (the density difference is 3400-3200=200 kg/m3), which corresponds to the Rayleigh number Ra c3 = 2x107. Thus, Rac3 in the oceanic crust changes sign at a depth of about 80 km.
The continental lithosphere has neutral buoyancy. Indexes 1, 2 and 3 denote the tracers of the substances of the continental crust, the continental lithosphere, and the oceanic crust, respectively. Concentrations of C1,2,3 can vary during the process of movement and mixing of the material from 1 to 0.
The viscosity of the continents is assumed to be increased in comparison with the temperature viscosity of the corresponding oceanic region n t. For the conti-
nental crust, the viscosity is taken in dimensionless quantities equal to nT+1000xC1, for the continental lithosphere it is nT+500xC2. Thus, the viscosity of the substance is transferred by the continental tracers (indices 1 and 2). The buoyancy of the substance in this version of the model is transferred by the tracers of the continental and oceanic crust (tracers 1 and 3).
The calculations in our study are conducted by the 2D Citcom code [Moresi, Gurnis, 1996; Zhong et al., 2000] with some updates and the automated graphics developed by A. Evseev. This code has been widely used and thoroughly tested. For each time instant, the momentum transfer equation and the heat transfer equation were solved. The thermal convection was calculated with the Rayleigh number Ra=2x107. The numerical finite element solution was computed in a box area with the aspect ratio L:D=5:1 on a uniform 401x201 grid, i.e., with the horizontal resolution of 36 km and the vertical resolution of 15 km. In the course of the computations, the model achieved the regime that is inherent to the assumed parameters (the systematical trend of the solution disappears). The parameter of artificial compressibility (tole_compressibi-lity) was assumed to be 5x10-6, and the accuracy of the Uzawa algorithm was 1x10-5. This combination of the parameters is optimal from the standpoint of the accuracy and performance [Brooks, Hughes, 1982; Hughes, 1987; Moresi, Gurnis, 1996]. The parameters of the algorithm are described in more detail in the above-mentioned papers.
3. Results and discussion
3.1. Oceanic crust with positive buoyancy
Figures 1-11 show the results for eleven successive stages of the convection. In each figure, the following fields are shown from top to bottom: the distributions of temperature and velocity, logarithmic viscosity, and the distribution of the three types of tracers in the calculation area (tracers of the continental crust are shown in red, the continental lithosphere is green, and the oceanic crust is blue.
The dimensionless time of stage 1 is assumed to be zero. The dimensionless times of stages 1-11 are t1=0; 12=2.63x10-4; t3=4.16x10-4; t4=8.88x10-4; t5=1.36x10-3; t6=1.82x10-3; t7=2.34x10-3; ts=3.01x10-3; t9=3.11x10-3; 110=3.25x10-3; tn=3.52x10-3, respectively. In the dimensional units, these values correspond to t1=0 Ma; 12=70 Ma; ts=110 Ma; t4=236 Ma; t5=361 Ma; te=483 Ma; t7=623 Ma; t8=802 Ma; t9=828 Ma; tM=864 Ma; 111=935 Ma, respectively. The calculation of this evolution required 6800 time-steps.
Our modeling shows the following characteristic stages:
^^Г^^Г" | | Logarithmic viscosity z -0.5 0 0.5 1 1.5 2 2.5 3 1-,--
0.5
1 I I I I I I I I I
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Markers у
Fig. 1. Mantle model (viscoplastic rheology), t=0. Spatial distribution patterns (from top to bottom): di-mensionless temperature and flow velocities; dimensionless viscosity (logarithmic viscosity scale); the spatial distribution of the three tracer fields: tracers of the continental crust are shown in red, the continental lithosphere is green, and the oceanic crust is blue.
Рис. 1. Модель мантии с вязко-пластической реологией, начальный момент времени (t=0). Представлены поля (сверху вниз): пространственное распределение безразмерной температуры и скоростей течений; пространственное распределение безразмерной вязкости (шкала вязкости логарифмическая); поля маркеров трех типов: маркеры континентальной коры показаны красным цветом, континентальной литосферы -зеленым, океанической коры -синим.
Fig. 2. Mantle model, t=70 Ma. See Fig. 1 for the legend.
Рис. 2. Модель мантии, второй момент времени (t=70 млн лет). Обозначения те же, что и для рис. 1.
Fig. 3. Mantle model, t=110 Ma. See Fig. 1 for the legend.
Рис. 3. Модель мантии, третий момент времени (t=110 млн лет). Обозначения те же, что и для рис. 1.
Fig. 4. Mantle model, t=236 Ma. See Fig. 1 for the legend.
Рис. 4. Модель мантии, t=236 млн лет. Обозначения те же, что и для рис. 1.
Fig. 5. Mantle model, t=361 Ma. See Fig. 1 for the legend.
Рис. 5. Модель мантии, t=361 млн лет. Обозначения те же, что и для рис. 1.
Fig. 6. Mantle model, t=483 Ma. See Z Fig. 1 for the legend. 1
Рис. 6. Модель мантии, t=483 млн лет. Обозначения те же, что и для рис. 1.
Fig. 7. Mantle model, t=623 Ma. See Fig. 1 for the legend.
Рис. 7. Модель мантии, t=623 млн лет. Обозначения те же, что и для рис. 1.
Initial moment, t=0. Continents are introduced in the model. Due to the impact of the mantle flows, the continents begin to move to the left in the model.
Fig. 2: t=70 Ma. The left slow-moving continent has almost stopped (as it reached the boundary of the calculated region). By this time, an inclined slab has formed under its edge (and the edge itself thickened).
Fig. 3: t=110 Ma. As soon as the left-side continent stopped, the motion of the right-side continent slowed down. Due to this slowing-down, the motion of the oceanic lithosphere material is hindered by the stalled plate, which leads to creation of a descending flow, i.e. a new subduction zone, behind it (x=1.8^2.0), or, in other variants, at some distance behind it. The closure of the ocean is almost maximal.
Fig. 4: t=236 Ma. Maximum convergence of the continents. The slab, formed by the oceanic plate between the continents, broke away under its own gravity (x=0.8^0.9, z=0.5^0.8), and subduction ceased. Thus, the scenario of slab detachment during the collision of two continents is realized, as described in [Davies, von Blanckenburg, 1995; Duretz et al., 2011; van Hunen, Allen, 2011]. The left-side continent experienced a total compression of ~10 % as compared to its initial state.
Then the opening of the ocean begins. The temperature and viscosity fields at this stage demonstrate
that a portion of the cold oceanic lithosphere has already broken off and is now plunging into the lower mantle.
Fig. 5: t=361 Ma. The opening of the ocean continues. The field of the tracers shows the role of the 660-km phase boundary: sinking of the oceanic crust is hindered when it passes this border, and the oceanic crust either lies on the phase boundary or is divided into two segments (right- or left-side subduction, respectively). The upper mantle flows begin to transfer the upper portion of the subducted oceanic crust into the opening ocean area (x=1.0^1.3, z=0.07^0.20).
Fig. 6: t=483 Ma. The portion of the oceanic crust remaining in the upper mantle has been already uplifted to the surface in the opening ocean area (x=1.2^1.3). An intensive descending mantle flow has been formed to the right of the right-side continent (x=3.2^3.8). Large velocity vectors show that it is the main driving force at this stage. The opening of the ocean between the continents is quite intense. Due to the impact of the slab formed at this stage, the right-ward motion of the right-side continent has also increased from 0.6 cm/yr in the previous stage to about 5.4 cm/yr. The slab retreats from the approaching continent to the right (the effect of backrolling) and moves slower that the continent.
| | И^И 11 IHJ mill и MI iij -0.5 0 0.5 1 1.5 2 2.5 3
-..... * - ■■ . 1 ч.
1 ll н У / i } l . ■ ' , V - f ... ry i
1 1 1 1 1 1....... 1 .......1............................
1-n—:-1-1-^-"1-1-1 I-1-—I-
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Markers x
Fig. 8. Mantle model, t=802 Ma. See Fig. 1 for the legend.
Рис. 8. Модель мантии, t=802 млн лет. Обозначения те же, что и для рис. 1.
Fig. 9. Mantle model, t=828 Ma. See Fig. 1 for the legend.
Рис. 9. Модель мантии, t=828 млн лет. Обозначения те же, что и для рис. 1.
Fig. 10. Mantle model, t=864 Ma. See Fig. 1 for the legend.
Рис. 10. Модель мантии, t=864 млн лет. Обозначения те же, что и для рис. 1.
Fig. 11. Mantle model, t=935 Ma. See Fig. 1 for the legend.
Рис. 11. Модель мантии, t=935 млн лет. Обозначения те же, что и для рис. 1.
Fig. 7: t=623 Ma. Opening of the ocean stops. Its size reached 7500 km (analogue of the Atlantic Ocean). During the entire period of opening, the left-side continent slowly stretches. The oceanic crust remaining in the upper mantle noticeably rises up due to convection in the upper mantle. The right-side continent overrided the slab and broke the upper-mantle portion of this slab, and the velocity of the continent decreased by an order of magnitude. The remnant slabs are located deep in the lower mantle (x=3.0^3.6, z=0.1^0.4). Simultaneously, at the continental sides, descending mantle flows formed, while the continent itself was shortened. The right-side descending flow is smaller and lies on the 660-km boundary for some time. The left-side descending flow (x=3.5^3.7) is very intense and continues to develop. It should be noted that an intensive descending flow in our model occurred thrice - it appears behind the moving continent at the moments when the continent motion is abruptly decelerated: t=110 Ma, t=623 Ma, and t=935 Ma.
Fig. 8: t=802 Ma.The left-side continent was separated from the left-side border of the calculation area 80 Ma before this stage. The right-side continent surrounded by subduction zones is under compression. The left-side ocean is opening; practically in the middle (x=0.7) the mid-oceanic ridge appeared, and here the recycled oceanic crust is transferred to the surface. At the front edge of the left-side continent, there is a fragment of the oceanic crust (x=2.08^2.15, z=0.9^1). This fragment joined the continent at an early stage of the evolution. Thus, the oceanic crust, attached to the continent during subduction, can remain at the surface for a long time (in the model, the full time of the numerical experiment is ~1 billion years). The right-side ocean is closing. The subduction zone (x=3.1^3.2) located in the closing ocean retreats to the right-side continent that is surrounded by the subduction zones and experiences compression.
Fig. 9: t=828 Ma. Active opening of the left-side ocean continues, and the velocity of the left continent reaches 9 cm/yr.
The left-side continent with the subduction zone reaches the mantle plume located in front of it (x=3.0^3.1). In the continental convergence region, there are three subduction zones, which are later combined into one.
Fig. 10: t=864 Ma. The moment of the continental collision. The left-side continent with the subduction zone destroyed the mantle plume in the closing ocean.
In this case, a group of three subduction zones in continental convergence region is combined into a single intensive mantle flow. Note that a similar phenomenon occurs in case of the Philippine oceanic plate that is surrounded by subduction zones.
In the region (x=0.8^1.0), there is a mid-oceanic ridge on the surface, where, as can be seen in the
figures, the recycled oceanic crust is transferred to the surface.
The left-side continent sharply slowed down during the collision. The subduction starts at its left edge. This effect occurs in this model for the third time.
Fig. 11: t=935 Ma. The rightward movement of the continents as a whole slows down as they approach the wall of the computational area. Now the continents are surrounded by subduction zones. The third subduction zone is located at the junction of the two continents.
Stretching in the left-side continent is initiated by a mantle flows underneath this continent. Note that in the upper part of the plunging slab, the values of the maximum shear stresses can reach values up to 200000 in dimensionless form, i.e. up to 50 MPa in dimensional form. Thus, the yield stress limit adopted in our model is reached in these region, which results in quasiplastic deformation of the material.
Summarizing the above-described results, we conclude that the model demonstrates all the basic elements of the modern theory of floating continents.
The positive buoyancy of the oceanic crust slows down its plunging, as well as the subduction of the oceanic slab (due to the viscous cohesion), if the crust does not detach from the slab.
In the areas where the oceanic crust thickness is comparable with the continental one (up to 30 km), the subduction of the slab slows down sharply. This effect is observed for basaltic plateaus in the ocean, such as, for example, the Ontong-Java plateau [Taira et al., 2004]. In detail, the dependence of crustal subduction on the crustal thickness and buoyancy was studied in [Trubitsyn et al., 2007].
The small thickness and low buoyancy of the oceanic crust as compared to the continental one explains the fundamental distinction in its evolution. The oceanic crust sinks in subduction zones, and only partially remains at the sides of the continents. The continental crust, on the contrary, practically does not subduct due to its large thickness and higher buoyancy (the density difference of the continental crust and the mantle is 2800-3200=-400 kg/m3, and for the oceanic crust it is -200 kg/m3). Thus, the oceanic crust exists on the surface for no more than 200 Ma (before the closure of the ocean), while the continental crust can exist for billions of years.
Our results show the constantly arising oscillations of the process (fluctuations in the shape of the flows, in their velocities, in heat flow, etc.), which vary in frequency. Using the experimental data on liquids [Dob-retsov et al., 1998; Kirdyashkin et al., 2000] obtained the (scaled) results which were correlated with the characteristic durations of the Wilson supercontinental cycle. For this purpose, the hydrodynamic homochronicity criterion Ho=ut/l was introduced, where u is the
0.5-
-*— I
1 Л \ i
0.5
1.5
2 2.5 Markers
3.5
4.5
5 X
Fig. 12. Mantle model (viscoplastic rheology of the oceanic lithosphere and eclogitization of the oceanic crust), t=236 Ma. Spatial distribution patterns (from top to bottom): dimensionless temperature and flow velocities; dimensionless viscosity (logarithmic viscosity scale); the spatial distribution of the three tracer fields: tracers of the continental crust are shown in red, the continental lithosphere is green, and the oceanic crust is blue.
Рис. 12. Модель мантии с вязкопластической реологией и эклогитизацией, t=236 млн лет. Представлены поля (сверху вниз): пространственное распределение безразмерной температуры и скоростей течений; пространственное распределение безразмерной вязкости (шкала вязкости логарифмическая); поля маркеров трех типов: маркеры континентальной коры показаны красным цветом, континентальной литосферы - зеленым, океанической коры - синим.
characteristic flow velocity, t is the period of temperature pulsations, and l is the thickness of the layer. For experiments with viscous fluids at Rayleigh numbers (6.8x105^1.0xl06), the estimate Ho=7.2 was obtained [Kirdyashkin et al., 2000].
Let us consider our results from this point of view. Taking the characteristic values from our calculations: u=4 cm/yr; tchar=500 million years - time comparable with the duration of the Wilson cycle; l=2900 km - the thickness of the mantle, we obtain Ho=7; the evaluation shows a good stability of this value for a wide range of models. Of course, there are shorter-period oscillations, in particular, with quasiperiodicity of 10-40 million years, associated with the emergence of individual plumes and the immersion of slabs.
3.2. Transition of basalt-eclogite in the oceanic crust
In this version of our calculations, the basalt-eclogite phase transition is included, i.e. RaC3 of the oceanic crust changes sign at a depth of about 80 km, and so oceanic crust becomes heavier than the surrounding mantle.
The calculations show a number of differences from the previous version. Figures 12-19 show some characteristic stages.
The initial stages during the first 250 Ma develop almost identically; then significant discrepancies begin to occur.
Unlike the first version of the model, where after 250 million years, a new ocean begins to open between
^^Z^^T | I Logarithmic viscosity
1-0.5 0 0.5 1 1.5 2 2.5 3 3.5
-- 7
--J ) У J с ---- /1.
1-1-1-1— .....T-1-1-1-1-f"""—\
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
Markers у
Fig. 13. Mantle model (viscoplastic rheology and eclogitization), t=490 Ma. See Fig. 12 for the legend.
Рис. 13. Модель мантии, t=490 млн лет. Обозначения те же, что и для предыдущих рисунков.
Fig. 14. Mantle model (viscoplastic rheology and eclogitization), t=635 Ma. See Fig. 12 for the legend.
Рис. 14. Модель мантии, t=635 млн лет. Обозначения те же, что и для предыдущих рисунков.
Fig. 15. Mantle model (viscoplastic rheology and eclogitization), t=710 Ma. See Fig. 12 for the legend.
Рис. 15. Модель мантии, t=710 млн лет. Обозначения те же, что и для предыдущих рисунков.
Fig. 16. Mantle model (viscoplastic rheology and eclogitization), t=845 Ma. See Fig. 12 for the legend.
Рис. 16. Модель мантии, t=845 млн лет. Обозначения те же, что и для предыдущих рисунков.
.^-^Ш^ИИ^^Е Logarithmic viscosity 1-0.5 0 0.5 1 1.5 2 2.5 3 3.5
7 2
//• k - г /А v^. / # ■ —A. -— А У J is* ttf* S ' у i ( 1 1 1 1
Markers x
Fig. 17. Mantle model (viscoplastic rheology and eclogitization), t=875 Ma. See Fig. 12 for the legend.
Рис. 17. Модель мантии, t=875 млн лет. Обозначения те же, что и для предыдущих рисунков.
Fig. 18. Mantle model (viscoplastic rheology and eclogitization), t=935 Ma. See Fig. 12 for the legend.
Рис. 18. Модель мантии, t=935 млн лет. Обозначения те же, что и для предыдущих рисунков.
Markers
| Fig. 19. Mantle model (viscoplastic rheology and eclogitization), t=990 Ma. See Fig. 12 for the legend. | Рис. 19. Модель мантии, t=990 млн лет. Обозначения те же, что и для предыдущих рисунков.
almost converged continents begins, here it does not occur (Fig. 12, t=236 Ma).
Until t=490 Ma, the continents remain almost converged and inactive. The gradual extension of the leftside continent by ~10 % takes place later on. Then a new ocean begins to open to the left of both continents, similar to the previous version.
Thus, the model without eclogitization shows the opening of the first ocean between continents at t=240 Ma and the second ocean to the left of both continents at t=720 Ma. According to the model with account of eclogitization, only the second ocean occurs: its opening begins at t=635 Ma.
It seems that the opening of the (first) ocean between the continents is impeded by a subvertically plunging strip of the eclogitized heavy oceanic crust beneath the suture of these continents (see Fig. 13: t=490 Ma, x=0.5^1.1, z=0^0.7). Hence, the sensitivity of the results to eclogitization is evident.
In the process of opening of the ocean, the gap between the two continents is completely closed (t=710 Ma), and the oceanic crust becomes sandwiched between continents.
Fig. 16: t=845 Ma. The recycled oceanic crust from the lower mantle begins to transfer to the surface in the opening ocean (blue tracers). Thus, the phenomenon of recirculation of the oceanic crust in the mantle [Christensen, Hofmann, 1994; Sobolev et al., 2007] is numerically confirmed.
Simultaneously, an important new effect arises: at the bottom of the lower mantle, two small clusters formed, containing the heavy material of the eclo-gitized oceanic crust. The upward movement of this heavy material occurs with difficulty (as seen in all subsequent stages of the calculation). The reason is the negative buoyancy of eclogite. In the previous version, the basal t-eclogite phase transition is not taken into account, so that there is no bulging of heavy material at the mantle bottom. Due to its positive buoyancy, the oceanic crust material is easily transferred to the surface of the Earth.
Fig. 17: t=875 Ma. The right-side continent reaches the mantle plume located in front of it. The descending flows fall sub-horizontally to the 660-km boundary.
The continent finally destroys the mantle plume (t=935 Ma), overriding it by the subduction zone.
When the movement of the united continents slows down, a slab appears behind them (t=990 Ma).
Summarizing, we can conclude that the above-discussed versions of the model taking into account eclogitization of the oceanic crust and without eclogiti-zation, show both different and common features of development. Common phenomena are: the opening of a new ocean; backrolling of the subduction zone (t=845 Ma 4 875 Ma); overriding the mantle plume (see t = 875 Ma 4 935 Ma); the emergence of slabs behind the continents experiencing decelerated motions (t=990 Ma).
Some authors (e.g. [Yoshida, 2010,2012; Magni et al., 2013]) do not take into account the buoyancy of the oceanic crust and the basalt-eclogite phase transition and consider this effect insignificant. One may assume that the influence of buoyancy of such a thin (7 km) layer of the oceanic crust material is insignificant. However, our calculations show that the inclusion of the basalt-eclogite transition effect in the oceanic crust significantly alters the pattern of mantle flows and the positions of continents. Moreover, the model with eclo-gitization shows a new effect - bulging of heavy material (eclogites, i.e. portions of the former oceanic crust) at the mantle bottom. These bulges move along the core-mantle boundary under the influence of the mantle convection, in the direction from slabs to mantle plumes. The process of their formation is connected with local plumes: the ascending plumes collect the material located on the bottom to the left and to the right of them and causes bulging. In the models and calculations that do not take the 'basalt-eclogite' transition into account, this phenomenon is absent.
Thus, our numerical experiments show that the inclusion of basalt-eclogite transition in the model considerably alters the evolution of mantle convection with continents. Moreover, a new effect arises - bulging of heavy material (eclogites, i.e. portions of the former oceanic crust) at the mantle bottom. The bulging forms several piles that are gradually brought upward by the plumes and thus appear on the surface of the Earth.
4. Conclusion
The main results of our numerical experiments include the following:
1. Our model of the two-dimensional mantle convection with non-Newtonian rheology in the presence of deformable continents demonstrates and allows to analyze all the basic elements of the modern theory of floating continents: convergence and compression of continents; the emergence of subduction zones; stretching, breakup and divergence of continents; opening and closing of the oceans; backrolling of sub-
duction zones; overriding of hot plumes by moving continents, and the recirculation of the oceanic crust in the mantle.
2. The model with the parameters accepted for our study shows considerable shortening (up to ~10 %) and thickening of the continents during their collision, as well as extension and thinning at other stages of their evolution. Stretching takes place only in the continental region located above the ascending mantle flow. The other regions of this continent may experience shortening in the same period of time.
3. When the moving continent slows down due to continental collision, it becomes an obstacle to the subhorizontal mantle flow behind it, and a descending flow (i.e. a new subduction zone) begins to form behind the continent. The new subduction zone may occur either at the edge of the continent or at some distance from it. Thus, the model predicts the moment when a new subduction zone arises.
4. The calculations based on the model with real parameters show that the continental crust does not sink, while the oceanic crust, despite its positive buoyancy at the surface, plunges into the mantle. This fundamental difference between the continental crust and the oceanic crust is the basis of the dual modes (the difference between continents and oceans) in the Earth geody-namics.
5. The phenomenon of recirculation of the oceanic crust in the mantle has been numerically confirmed. At the same time, some part of the oceanic crust remains attached to the continental margins for a long time.
6. The effect of buoyancy of the thin (7 km) oceanic crust layer is significant. According to our calculations, the effect of eclogitization is important for convection in the entire mantle. The inclusion of basalt-eclogite transition in the model show an additional new effect -bulging of heavy material (eclogitized remnants of the oceanic crust) at the core-mantle boundary, wherefrom it can be transferred upward by the mantle plumes and thus appear on the surface of the Earth.
Of course, the processes in reality are much more complicated. In the future numerical experiments, we plan to include new additional factors in the model.
5. Acknowledgments
The work was supported by the Russian Foundation for Basic Research (projects 16-55-12033 and 13-0501123). The authors would like to thank the two anonymous reviewers, whose thorough reviews helped to improve the manuscript significantly. We are grateful to Louis Moresi, Shijie Zhong, Michael Gurnis, and other authors of the 2D CITCOM code for providing the possibility of using this software.
6. References
Bobrov A.M., Baranov A.A., 2011. Horizontal stresses in the mantle and in the moving continent for the model of two dimensional convection with varying viscosity. lzvestiya, Physics of the Solid Earth 47 (9), 801-815. https:// doi.org/10.1134/S1069351311090023.
Bobrov A.M., Baranov A.A., 2016. The mantle convection model with non-Newtonian rheology and phase transitions: the flow structure and stress fields. lzvestiya, Physics of the Solid Earth 52 (1), 129-143. https:// doi.org/10.1134/S1069351316010031.
Bobrov A.M., Trubitsyn A.P., 2008. Numerical model of the supercontinental cycle stages: integral transfer of the oceanic crust material and mantle viscous shear stresses. Studia Geophysica et Geodaetica 52 (1), 87-100. https:// doi.org/10.1007/s11200-008-0007-1.
Brooks A.N., Hughes T.J.R., 1982. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Computer Methods in Applied Mechanics and Engineering 32 (1-3), 199-259. https://doi.org/10.1016/0045-7825(82)90071-8.
Burov E., Francois T., Agarda P., Le Pourhiet L., Meyer B., Tirel C., Lebedev S., Yamato P., Brun J.-P., 2014. Rheological and geodynamic controls on the mechanisms of subduction and HP/UHP exhumation of crustal rocks during continental collision: Insights from numerical models. Tectonophysics 631, 212-250. https://doi.org/10.1016/j.tecto.2014. 04.033.
Butler S.L., Jarvis G.T., 2004. Stresses induced in continental lithospheres by axisymmetric spherical convection. Geophysical Journal International 157 (3), 1359-1376. https://doi.org/10.1111/j.1365-246X.2004.02257.x.
Christensen U., Hofmann A., 1994. Segregation of subducted oceanic crust in the converting mantle. Journal of Geophysical Research: Solid Earth 99 (B10), 19867-19884. https://doi.org/10.1029/93JB03403.
DaviesJ.H., von Blanckenburg F., 1995. Slab breakoff: A model of lithosphere detachment and its test in the magmatism and deformation of collisional orogens. Earth and Planetary Science Letters 129 (1-4), 85-102. https://doi.org/ 10.1016/0012-821X(94)00237-S.
Dobretsov N.L., Kirdyashkin A.A., Kirdyashkin A.G., Popov S.P., 1998. Temporal characteristics of nonstationary free-convective flows in the horizontal layer and time scales of lower mantle convection. Doklady Earth Sciences 363 (8), 1132-1135.
Duretz T., Gerya T.V., May D.A., 2011. Numerical modelling of spontaneous slab breakoff and subsequent topographic response. Tectonophysics 502 (1-2), 244-256. https://doi.org/10.1016/j.tecto.2010.05.024.
Fei Y., Orman J.V., Li J., van Westrenen W., Sanloup C., Minarik W., Hirose K., Komabayashi T., Walter M., Funakoshi K., 2004. Experimentally determined postspinel transformation boundary in Mg2SiO4 using MgO as an internal pressure standard and its geophysical implications. Journal of Geophysical Research: Solid Earth 109 (B2), B02305. https://doi.org/10.1029/2003JB002562.
Gurnis M., 1988. Large-scale mantle convection and aggregation and dispersal of supercontinents. Nature 332 (6166), 696-699. https://doi.org/10.1038/332695a0.
Hughes T.J.R., 1987. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Prentice-Hall, Englewood Cliffs, New Jersey, 803 p.
Karato S., Wu P., 1993. Rheology of the upper mantle: a synthesis. Science 260 (5109), 771-778. https://doi.org/ 10.1126/science.260.5109.771.
Kirdyashkin A.A., Dobretsov N.L., Kirdyashkin A.G., 2000. Experimental modeling of the influence of subduction zones on the spatial structure of lower mantle convection and characteristic periods of heat flow fluctuations in the mantle. Doklady Earth Sciences 371A (3), 565-568.
Lobkovskii L.l., Inyukhin A.V., Kotelkin V.D., 2014. Subduction and cyclic processes in the upper mantle. Doklady Earth Sciences 459 (1), 1348-1352. https://doi.org/10.1134/S1028334X14110233.
Lobkovsky L.l., Kotelkin V.D., 2004. Numerical analysis of geodynamic evolution of the Earth based on a thermochemi-cal model of the mantle convection. Russian Journal of Earth Sciences 6 (1), 49-58.
Lobkovsky L., Kotelkin V., 2015. The history of supercontinents and oceans from the standpoint of thermochemical mantle convection. Precambrian Research 259, 262-277. https://doi.org/10.1016Zj.precamres.2015.01.005.
Lowman J.P., King S.D., Gable C.W., 2004. Steady plumes in viscously stratified, vigorously converting, three-dimensional numerical mantle convection models with mobile plates. Geochemistry, Geophysics, Geosystems 5 (1), Q01L01. https://doi.org/10.1029/2003GC000583.
Magni V., Faccenna C., van Hunen J., Funiciello F, 2013. Delamination vs. break-off: the fate of continental collision. Geophysical Research Letters 40 (2), 285-289. https://doi.org/10.1002/grl.50090.
Moresi L.N., Gurnis M., 1996. Constraints on the lateral strength of slabs from three dimensional dynamic flow models. Earth and Planetary Science Letters 138 (1-4), 15-28. https://doi.org/10.1016/0012-821X(95)00221-W.
Moresi L.N., Solomatov V., 1998. Mantle convection with a brittle lithosphere: Thoughts on the global tectonic styles of the Earth and Venus. Geophysical Journal International 133 (3), 669-682. https://doi.org/10.1046/j.1365-246X. 1998.00521.x.
Schubert G., Turcotte D.L., Olson P., 2001. Mantle Convection in the Earth and Planets. Cambridge University Press, New York, 940 p.
Simon K., Huismans R.S., Beaumont C, 2009. Dynamical modeling of lithospheric extension and small-scale convection: implications for magmatism during the formation of volcanic rifted margins. Geophysical Journal International 176 (1), 327-350. https://doi.org/10.1111/j.1365-246X.2008.03891.x.
Sobolev A.V., Hoffman A.W., Kuzmin D.V., Yaxley G.M., Arndt N.T., Chung S-L, Danyushevsky L.V., Elliott T, Frey F.A., Garcia M.O., Gurenko A.A., Kamenetsky V.S., Kerr A.C., Krivolutskaya N.A., Matvienkov V.V., Nikogosian I.K., Rocholl A., Sigurdsson I.A., Sushchevskaya N.N., Teklay M., 2007. The amount of recycled crust in sources of mantle-derived melts. Science 316 (5823), 412-417. https://doi.org/10.1126/science.%201138113.
Taira A., Mann P., Rahardiavan R., 2004. Incipent subduction of the Ontong Java plateau along the North Colomon trench. Tectonophysics 389 (3-4), 247-266. https://doi.org/10.1016Zj.tecto.2004.07.052.
Tosi N., Yuen D.A., 2011. Bent-shaped plumes and horizontal channel flow beneath the 660 km discontinuity. Earth and Planetary Science Letters 312 (3-4), 348-359. https://doi.org/10.1016/j.epsl.2011.10.015.
Trubitsyn V., 2012. Rheology of the mantle and tectonics of the oceanic lithospheric plates. Izvestiya, Physics of the Solid Earth 48 (6), 467-485. https://doi.org/10.1134/S1069351312060079.
Trubitsyn V., Baranov A., Kharybin E., 2007. Numerical models of subduction of the oceanic crust with basaltic plateaus. Izvestiya, Physics of the Solid Earth 43 (7), 533-542. https://doi.org/10.1134/S1069351307070014.
Trubitsyn V.P., Rykov V.V., Jacoby W.R., 1999. A self-consistent 2-D model for the dip angle of mantle downflow beneath an overriding continent. Journal of Geodynamics 28 (2-3), 215-224. https://doi.org/10.1016/S0264-3707(98) 00038-6.
Van Hunen J., Allen M.B., 2011. Continental collision and slab break-off: A comparison of 3-D numerical models with observations. Earth and Planetary Science Letters 302 (1-2), 27-37. https://doi.org/10.1016/j.epsl.2010.11.035.
Yoshida M., 2010. Preliminary three-dimensional model of mantle convection with deformable, mobile continental lithosphere. Earth and Planetary Science Letters 295 (1-2), 205-218. https://doi.org/10.1016/j.epsl.2010.04.001.
Yoshida M., 2012. Dynamic role of the rheological contrast between cratonic and oceanic lithospheres in the longevity of cratonic lithosphere: A three-dimensional numerical study. Tectonophysics 532-535, 156-166. https://doi.org/ 10.1016/j.tecto.2012.01.029.
Zhong S., Zuber M.T., Moresi L.N., Gurnis M., 2000. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection. Journal of Geophysical Research: Solid Earth 105 (B5), 11063-11082. https://doi.org/10.1029/2000JB900003.
Alexandr M. Bobrov, Candidate of Physics and Mathematics, Lead Researcher O.Yu. Schmidt Institute of Physics of the Earth of RAS 10 Bol'shaya Gruzinskaya street, Moscow D-242 123242, GSP-5, Russia И e-mail: a_m_bobrov@yahoo.com
Александр Марович Бобров, канд. физ.-мат. наук, в.н.с. Институт физики Земли им. О.Ю. Шмидта РАН 123242, ГСП-5, Москва Д-242, ул. Большая Грузинская, 10, Россия И e-mail: a_m_bobrov@yahoo.com
Aleksei A. Baranov, Candidate of Physics and Mathematics, Lead Researcher O.Yu. Schmidt Institute of Physics of the Earth of RAS
10 Bol'shaya Gruzinskaya street, Moscow D-242 123242, GSP-5, Russia Institute of Earthquake Prediction Theory and Mathematical Geophysics RAS
84/32 Profsoyuznaya street, Moscow 117485, Russia e-mail: Baranov@ifz.ru
Алексей Андреевич Баранов, канд. физ.-мат. наук, в.н.с. Институт физики Земли им. О.Ю. Шмидта РАН
123242, ГСП-5, Москва Д-242, ул. Большая Грузинская, 10, Россия Институт теории прогноза землетрясений и математической геофизики РАН
117485, Москва, ул. Профсоюзная, 84/32, Россия e-mail: Baranov@ifz.ru
