Evolution of stress fields during the supercontinent cycle

Alexnder Bobrov , Alexey Brnov ,b,*, Robert Tenzer

a Schmidt Institute of Physics of the Earth, Russian Academy of Sciences, Moscow, Russia

b Institute of Earthquake Prediction Theory and Mathematical Geophysics, Russian Academy of Sciences, Moscow, Russia

c Department of Land Surveying and Geo-Informatics, Hong Kong Polytechnic University, Hong Kong, China

Keywords:Supercontinent cycle Floating deformable continents Thermochemical convection Horizontal stresses Dynamic topography

ABSTRACT We investigate the evolution of stress fields during the supercontinent cycle using the 2D Cartesian geometry model of thermochemical convection with the non-Newtonian rheology in the presence of floating deformable continents.In the course of the simulation,the supercontinent cycle is implemented several times.The number of continents considered in our model as a function of time oscillates around 3. The lifetime of a supercontinent depends on its dimension. Our results suggest that immediately before a supercontinent breakup,the over-lithostatic horizontal stresses in it(referring to the mean value by the computational area)are tensile and can reach-250 MPa.At the same time,a vast area beneath a supercontinent with an upward flow exhibits clearly the over-lithostatic compressive horizontal stresses of 50-100 МРа.The reason for the difference in stresses in the supercontinent and the underlying mantle is a sharp difference in their viscosity.In large parts of the mantle,the over-lithostatic horizontal stresses are in the range of ±25 MPa, while the horizontal stresses along subduction zones and continental margins are significantly larger.During the process of continent-to-continent collisions,the compressive stresses can approximately reach 130 MPa, while within the subcontinental mantle, the tensile overlithostatic stresses are about -50 MPa. The dynamic topography reflects the main features of the supercontinent cycle and correlates with real ones. Before the breakup and immediately after the disintegration of the supercontinent,continents experience maximum uplift.During the supercontinent cycle,topographic heights of continents typically vary within the interval of about±1.5 km,relatively to a mean value. Topographic maxima of orogenic formations to about 2-4 km are detected along continent-tocontinent collisions as well as when adjacent subduction zones interact with continental margins.

1. Introduction

Earth is currently the only planet in our solar system with the supercontinent cycle [1,2]. The aggregation and breakup of supercontinents is a fundamental process in the Earth's evolution. This cycle controls the growth of global mountain belts[3],continental uplift and sea level [4], flood basalt volcanism [5,6], global heat balance,and climate[7-9].Nowadays,the role of continents in the Earth's evolution is rather well understood. A detailed review of supercontinents was given, for instance, in [2,10].

A number of authors used numerical modeling techniques to study the nonlinear influence of mantle convection on a supercontinent cycle.2D models were used by Gurnis[11],Lowman and Jarvis [12], Lowman and Jarvis [13], Trubitsyn and Bobrov [14],Lowman and Jarvis [15], Honda et al. [16], Butler and Jarvis [17],Bobrov and Trubitsyn[18],Neil et al.[19],Heron and Lowman[20],Lobkovsky et al. [21], Rolf et al. [22], Trim and Lowman [23], Dal Zilio et al. [24], Kameyama and Harada [25], Bobrov and Baranov[26,27], Jain et al. [28], Mao et al. [29], and others. 3D mantle convection models were used to investigate the supercontinent cycle by Rykov and Trubitsyn[30],Lowman and Gable[31],Yoshida et al. [32], Honda et al. [16], Phillips and Bunge [33], Phillips and Bunge [34], Zhong et al. [35], Li and Zhong [36], Zhang et al. [37],Yoshida [38], Yoshida [39], Heron and Lowman [20], Yoshida [40],Yoshida [41], Yoshida and Santosh [42], Lobkovsky and Kotelkin[43], Zhang et al. [44], Mao et al. [29], Yoshida [45] and others. In addition to theoretical models, experimental studies were also conducted to investigate the mantle convection pattern with moving continents. For instance, Guillou and Jaupart [46] used silicon oil while involving a temperature-dependent viscosity. The classification of existing studies of supercontinent cycle depending on applied models(2D,3D),stresses,type of supercontinent cycle,and main results is given in Table 1.

One of the most important parameters of the supercontinent cycle is the stress field. When the ultimate strength of the supercontinent is exceeded,supercontinents break up.Moreover,it would be impossible for the continents to converge and assemble without exceeding the ultimate strength of the oceanic plate and its subduction. Other parameters (including various stress fields) provide not only a piece of additional information about a supercontinent cycle but also allow excluding models that are plausible based only on commonly used fields and,conversely,select acceptable models.

Stresses during the supercontinent cycle have been investigated in numerous studies.Gurnis[11],using a 2D model,estimated that the supercontinent tensile strength was about 50-70 MPa. Lowman and Jarvis [13,15], using a 2D model, suggested that critical stress for supercontinent breakup is 80 MPa,and this value could be reached with or without considering internal heating.Lowman and Jarvis [12] used a more moderate value of the supercontinent tensile strength of 48 MPa with two rigid continents and oceanic plates.Bobrov and Baranov[47]suggested that continental stresses strongly depend on the position of the continent with respect to subduction zones. Under this assumption, the over-lithostatic pressure in continents ranges from -5 to 15 MPa; the horizontal over-lithostatic tensile stress amounts up to -4 MPa. Butler and Jarvis [17] investigated various axisymmetric spherical mantle models with immobile supercontinent represented by a highviscosity area. According to their results, the tensional stresses in the supercontinent change within the range from 25 to 72 MPa.Yoshida [39] applied a spherical model with a temperaturedependent viscosity in the mantle, superimposed by an immobile supercontinent. He found out that, depending on the model parameters,the maximum deviatoric tensile stresses generated in the supercontinent were 30-90 MPa.

Rolf et al.[22],Phillips and Bunge[34],and Bobrov and Baranov[27] demonstrated the existence of irregularities in the supercontinent cycle. For the supercontinent, typical shear stresses change in a wide range of 50-200 MPa. Before the breakup, maximum shear stresses generated in the supercontinent reach 200 MPa[27].In a more recent spherical study, Yoshida [45] suggested that stresses under continents could reach 100 MPa.

Although numerous studies have investigated stresses in the mantle and continents during some stages of the supercontinent cycle,we argue that the stress distribution in time and space is not yet fully understood, particularly the evolution of stresses during the supercontinent cycle. The stresses required to break up a supercontinent should be generated by mantle convection while considering various parameters, such as rheology, the intensity of mantle convection, and internal heat. Open questions also remain about the contribution of ascending flows to a supercontinent breakup and the necessity of pre-existing weak zones for the breakup.To address these aspects,we investigated the stresses and dynamic topography during the supercontinent cycle by using a 2D Cartesian model and considering a non-Newtonian rheology. The continents and oceanic crust were modelled by active tracers(chemical heterogeneities)with their own viscosities and densities.Phase transitions at 410 and 660 km, internal and bottom heating were included in the model. Moreover, beneath the oceanic crust,we added the phase transition basalt-eclogite at the depth of about 80 km. We also assumed that motion of continents is consistent with a mantle flow pattern, meaning that the velocity of continental motions varies in the course of their motion according to time-dependent forces that act from the viscous mantle.The study is organized into three sections. Theoretical model is explained in section 2.Results are presented in section 3.The main findings are discussed, and the study is concluded in section 4.

2. Theoretical model

In our model, the Earth's mantle is modeled as a Boussinesq fluid with an infinite Prandtl number.The mantle is assumed to be heated from the core,and due to the decay of radioactive elements in the mantle itself(internal heating).Under these assumptions,the thermochemical convection with different compositions is governed by the dimensionless equations for conservation of mass,momentum, energy, and advection. For a 2D model, these equations are defined in the following form [26,27]:

where Vand Vare velocities; p is the dynamic pressure; τ, τand τare components of the deviatoric tensor of viscous stresses;T is the temperature; C is the concentration; Γ is the phase function; H is the internal heat production; Ra is the thermal Rayleigh number;Rais the phase Rayleigh number and Rais the chemical Rayleigh number.

Since the mantle is considered to be incompressible, the expression in Eq. (1) describes a volume conservation. The momentum equations in Eqs. (2) and (3) postulate the Newton's second law as a force balance between viscous forces, dynamic pressure, and the thermal buoyancy force. The internal heat production H in our model is set to be 15.

The thermal Rayleigh number in Eq. (3) is defined by

Table 1 A classifciation of studies on the supercontinent cycle.

Table 1 (continued)

mantle. Themodelshowsshortening (upto ～10%)and thickeningof thecontinents during theircollision.The influence of thesize of thesupercontinent on convectionpatternhasbeen studied.Itwasshownan irregularity of supercontinent cycle,typical velocitiesforcontinents before collisionare3-10 cm/year,for supercontinent 0.5-1.5cm/yearandafterthebreakup4-8cm/year Decreasinginfluence of continentalinsulation on subcontinentalwarmingwith increasing Ra andtimedependent correlationbetweencontinents andthe underlyingmantle arefound.The extroversion process(i.e., theclosureof thePacific)is dominantwhen viscosityincreasesby afactor of 30 from the upperto thelowermantle,whereastheintroversion process (i.e., theclosureof theAtlantic)is dominant when factor=100.Immobile supercontinent Incompletecycle,fivecontinents Incompletecycle,two continents or one supercontinent Incompletecycle,formation of anew Earth's supercontinent Beforethebreakup,maximum shear stressgeneratedin the supercontinent can reach200MPa Tensional and compressional stressacting under the moving continentsreach 100 MPa continents. Continents were modelled by theactive markerswith increased viscosityandlowdensity.2Dand3D mantle isoviscous models,Ra=1×106.Rigidcontinents are modeledby rectangles.2D(1:10) modelwith p, T, andstress dependentviscosity,Ra=2.5×107 withphasetransitions,oceaniccrustanddeformable continents. Continents were modelled by theactive markerswith increased viscosityandlowdensity.Spherical annulus(2D),sphericalshell, with p, T, stress-dependent viscosity.The continents aremodelled as acompositionallybuoyant and viscoustracers.Spherical modelwith p,T,stress dependentviscosityin themantle,Ra=5.89 ×107 withphasetransitionsat 400and670km.Continentsaremodelled by acompositionallybuoyantandviscous tracers.Mao et al.[29]BobrovandBaranov[27]Jainet al.[28]Yoshida [45]29 30 31 32

where Δρ is the density contrast between two chemical components.

We have considered three chemical components,namely for the continental crust, the continental lithosphere, and the oceanic crust, which are modeled by active buoyant and viscous tracers with thickness values of 40, 90, and 7 km and density values of ρ= 2800, ρ= 3200 and ρ= 3000 kg·mrespectively. The continental lithosphere has a neutral buoyancy (ρ= ρ=3200 kg·m),whereas the oceanic crustal density changes from 3000 to 3400 kg·mdue to the basalt-eclogite phase transition at a depth of about 80 km[48].In this case,Ra= -4 ×10,Ra= 0, Ra= -2×10(above 80 km) and Ra= 2×10(deeper than 80 km). The phase Rayleigh number Rafor phase transitions at 410 km and 660 km is computed from

where ηis the compositional prefactor, the activation energy is E=10.36,Tis equal to 0.6,and the temperature at the lower bound of the mantle is T= 1 [38]. The viscosity jump between the upper and lower mantle is 30 according to Zhong et al.[35];see also Zhang et al. [37]. For this model, the viscosity at the surface becomes very high,and the oceanic lithosphere ceases to subduct.To model the rheology close to the real Earth,it is necessary to take into account the dependence of the viscosity of the oceanic lithosphere on the stress in it. For the oceanic lithosphere (above 150 km),we took into account the plastic character of deformations so that

where τ is the effective yield stress [38]. For more details about theoretical definitions summarized above, we refer readers to the study by Bobrov and Baranov [27]. We assume that the strength(i.e., the yield stress) of the oceanic lithosphere is 50 MPa [49],whereas the strength of continents is accepted to be infinite [42].Such rheology is sufficient to generate rigid oceanic plates separated by narrow weak plate boundaries.Moreover,the subduction of rigid oceanic plates is permissible.

The compositional prefactor ηis defined for the continental crust and continental lithosphere to maintain the integrity of the continent.The dimensionless viscosity of the continental crust and lithosphere depends on the composition. We then write

We used the CitcomS code [50] with some improvements(active tracers) for solving the expressions in Eqs. (1)-(5). The numerical finite element solution was computed in a box area with the aspect ratio L:D=10:1 on a uniform 801×201 grid,i.e.,with the horizontal resolution of 36 km and the vertical resolution of 15 km.

In initial computations, the pure thermal convection model reaches the regime that is inherent to the assumed parameters(the systematical trend of the solution disappears).Then,we introduced five continents into the model, with the initial time set to be t=0(Fig. 1). The positions of continent are: the first left continent x=0.05-1,z=0.0-0.045;the second x=4.1-5.2,z=0.0-0.045;the third x = 6.2-6.4, z = 0.0-0.045; the fourth x = 7.3-7.7,z = 0.0-0.045; the fifth x = 9.3-9.6, z = 0.0-0.045. Under the action of mantle flows, continents began to drift and unite in groups and after disintegrate.Here and further by supercontinent,we mean the united group of continents.

3. Results

3.1. Evolution of fields arising in the process of convection

Mantle flows, temperature, viscosity as well as the distribution of normal vertical stresses on the surface of the computational domain are shown at various stages of the supercontinent cycle(Figs.2-7).We also computed the normal horizontal stresses(σ)within the continental crust and the underlying mantle.The values of σare given relative to the average value. Figs. 2-4 show the process of asymmetric breakup of the supercontinent, while at stages shown in Figs. 5-7, the breakup proceeds symmetrically.Simultaneously in other areas shown in Figs. 2-4, the process of continental aggregation occurs.

The temperature field possesses significant thermal inertia and shows shapes of plumes, slabs, their delay or passage through the 660-km phase boundary, their gradually warming up remnants at the bottom of the lower mantle, etc. The viscosity field shows a jump at the 660-km phase boundary as well as the temperaturedependent nature of viscosity, i.e., the increase of viscosity with depth.The field also manifests small areas in the uppermost part of subduction zones, where the viscosity drops because of the softening effect caused by high stresses.Black sub-horizontal lines with small elevations and depressions indicate the position of the 410-km and 660-km phase boundaries. The analysis of the presented results is done next.

3.1.1. Evolution of normal horizontal stresses

Repeating assemblies and breakups of continents in the model give the possibility to reveal some typical features of the supercontinent cycle.

Fig.1. The mantle model, with the initial time t = 0. Spatial distributions of dimensionless temperature and flow velocities (top); maximum velocities are about 8 cm per year.Spatial distributions of dimensionless viscosity(logarithmic viscosity scale)and flow velocities(bottom).The continents are presented as a purple high-viscous area on the surface.

Fig. 2. The mantle model, the stage t = 1.440 Ga. (a) Spatial distributions of dynamic topography on the surface; (b) Spatial distributions of dimensionless temperature and flow velocities; (c) Spatial distributions of dimensionless logarithmic viscosity and flow velocities (the continents are presented as a purple high-viscous area on the surface); (d)Dimensionless normal horizontal stresses σxx. The black lines depict the locations of phase boundaries.

As seen in Fig. 2, left and right, two supercontinents are at different stages.The large left supercontinent,x=3.75-5.85,is in a state before the breakup. The stresses in it are tensile. The smaller right-hand supercontinent,x=6.95-7.60,was formed recently and is in a state of compression. At stages before the beginning of a supercontinent breakup, we find the state of pronounced horizontal tensile (over-lithostatic) stresses (violet and blue colors in the lower panel of figures) throughout a supercontinent (stages t= 1.440 Ga,Fig.2; and t= 2.882 Ga,Fig.5).At the same time,an extensive area of subcontinental mantle surrounding the mantle upward flow exhibits clearly the horizontal compressive stresses(red and orange colors).In the area beneath the supercontinent in Fig. 2, for instance, the over-lithostatic compressive stresses range from 60 to 100 МРа. We remind here that the value of 400 000 in the dimensionless scale shown in the figures corresponds to 113 МРа. At the same time, the over-lithostatic stresses in the supercontinent are tensile and can reach -250 МPа just before its breakup. The cause of this stress differences in the supercontinent and in the underlying mantle is a sharp difference in the viscosities of the continental crustal and the underlying mantle materials. In other words, a supercontinental material cannot move like the underlying mantle material because of its high viscosity.The effect is visible in Figs. 2 and 5, although the mantle flows under the supercontinent at the stage shown in Fig. 2 are significantly more asymmetric than at the stage shown in Fig. 5. However, immediately with the beginning of a supercontinent breakup, the overlithostatic tensile stresses in the supercontinent almost everywhere change to compressive ones (Fig.3, stage t = 1.505 Ga; and Fig.6,t=2.980 Ga for the next cycle).This process occurs under the influence of actively intruding hot plume material into the opening gap. At the stages of continental assembly, an opposite contrast of the stresses in continents and underlying mantle takes place. For example,at stage t=1.701 Ga(see Fig.4),a supercontinent that just was formed by the collision of two continents above a downwardmoving slab exhibits the horizontal compressive stresses,while the tensile (relative to the mean value) stresses take place in the vast mantle area below it. An analogous case is presented at the aforementioned stage t = 1.440 Ga (see Fig. 2) where a small (compound) continent is located in the right side of the computational domain above the remnants of two descending mantle flows.Beneath this small supercontinent, the over-lithostatic tensile stresses аге about -60 МРа (about -213 000 in dimensionless units), while the supercontinent is characterized by compressive stresses of about 130 МРа. The boundary of a sharp change in the values of σ, which appears in these areas and at these stages,outlines the position of effective viscosity jump,meaning that this interface marks the base of the high-viscosity lithosphere.

The ascending flows that rise not in the subcontinental region but the oceanic zone show significantly lower stresses(of the order of 100 000 in dimensionless form, that is, of the order of 28 МРа),occupying much smaller areas.It is important to note that the same phenomenon of sharp contrast of stresses is also presented in oceanic areas and two variants.In the first case,the over-lithostatic tensile stresses arise in a relatively thin oceanic lithosphere above the layer with upward-moving flows where the horizontal stresses are compressive. In particular, this can be seen in the uppermost part of the computational domain in Fig. 5 at x = 4.9-6.8 at the stage t=3.053 Ga,and in Fig.7 at x=2.5-3.5 and x=4.5-5.3.In the second case, an opposite effect takes place. The horizontal compressive stresses in the oceanic lithosphere along with tensile or close-to-zero horizontal stresses in the mantle are seen at the stage shown in Fig. 6: the uppermost part of the computational domain, x = 4.1-5.1 where the rigid oceanic lithosphere moves onto a subduction zone at x = 4.1; and at the stage of Fig. 3:x=6.05-6.7.Thus,the difference in viscosities of 3 orders and even 2 orders is sufficient for the appearance of this phenomenon. The general type of horizontal stresses in the continents (tensile or compressive), as well as the distribution of stresses along the continent, is determined by the resulting effect of three factors -the normal horizontal stresses acting on the left- and right-end faces of the continent, and shear stresses distributed along the lower boundary. Therefore, their combined contribution may not be obvious.

Fig. 3. The mantle model, the stage t = 1.505 Ga. See Fig. 2 for the legend.

Fig. 4. The mantle model, the stage t = 1.701 Ga. See Fig. 2 for the legend.

Fig. 5. The mantle model, the stage t = 2.882 Ga. See Fig. 2 for the legend.

Fig. 6. The mantle model, the stage t = 2.980 Ga. See Fig. 2 for the legend.

Fig. 7. The mantle model, the stage t = 3.053 Ga. See Fig. 2 for the legend.

3.1.2. Normal vertical stresses and surface topography

We calculated the vertical stress field σin the computational domain, particularly on its surface. The stresses on the surface make it possible to determine its topography, assuming that this dynamic relief is supported by the found stresses.In the upper parts of Figs. 2-7, the scale of surface stresses σis shown in dimensionless units as well as the scale of the corresponding dynamic topography in dimensional form (in km).

Our results demonstrate a systematic trend of changes in the topography of continents during the supercontinent cycle. Fig. 2 reveals that before the breakup, when ascending mantle flows take place under the supercontinent, the left half of the supercontinent is on average about 2 km higher than the adjacent oceanic region,and its right half is also partially uplifted.After the beginning of breakup(Fig.3),the height of the left continent(with respect to the adjacent oceanic areas) decreases by about 1 km on average, and the subsidence of the right continent occurs with a delay because of the asymmetric disintegration of the supercontinent in this case. Later, in the course of further divergence of continents, their heights continue to gradually decrease (the stage t = 1.701 Ga,Fig. 4).

Descending mantle flows create an opposite effect. At the stage t=1.440 Ga(Fig.2),the small continent(x=7.0-7.6)located above the descending mantle flows is lowered under their influence and,on average does not differ by height from the neighboring oceanic region.Here we have in mind that it is necessary to add the inherent height of the continents above sea level to the values given here.

The region of a new opening ocean at the stage t = 1.505 Ga(Fig.3)raised as much as 2-4 km concerning to old oceanic areas.Thus, the simulation demonstrates that a new expanding ocean should inevitably be shallow. A similar evolutionary trend takes place during the next disintegration of a supercontinent and the opening of a new oceanic region(Figs. 5-7). Heights of the supercontinent also reach maxima before the beginning of its breakup.Later, the heights of diverging continents gradually decrease(especially near subduction zones). The second case differs from the first case of supercontinent splitting in a more symmetrical structure, while in the first case, the ascending mantle jet is deflected to the left, and the opening occurs mainly in the same direction.As seen in Fig.2,a narrow deep depression appears in the initial stage of the rift opening.The reason is that the divergence of parts of the supercontinent (caused by forces on the base and,possibly,downward flows at the outer ends)has already begun,but filling the opening gap with the mantle material is hindered due to frictional forces on margins of this still narrow channel.

During the continental collision,an orogenic area with a height of several kilometers appears between them (see Fig. 4). The orogenic areas with the heights of 2-4 km also occur when adjacent subduction zones interact with the continental margins (see Figs.2-4).The presented plots also show that the local extrema of vertical stresses σand topography of the continents coincide with the subcontinental upward and downward mantle flows, being located directly above them. Above the uppermost part of the ascending flow, the relief of the continent is elevated, while lowered above the downdrafts. We note that the distribution of horizontal stresses σalso marks the locations of subcontinental downward and upward flows (the latter can be seen, for example,in the bottom panel of Fig. 2 in the left continent at x = 3.9-4.2).

3.2. Cyclicity and characteristic model times of the continental aggregation and dispersal

To reveal characteristic features of oscillations in the number of continents, we formally carried out evolutionary calculations covering a time interval of almost twice as long as the Earth's existence. In particular, these calculations show the repeating assembly of continents and their dispersal. In most cases, the assembly and dispersal of continents occurred in the central area instead of boundaries of the computational domain which could influence the process.

The number of separate continents as a function of time in the model is plotted in Fig. 8. We see the presence of four separate continents or their assemblies during 11.7% of the total time of calculation, the presence of three continents during 67.7% of the total calculation time,and two continents during 17.2%of the total time. Thus, the number of continents on the surface of the model oscillates around 3.In fact,as the third(small)continent almost all the time is immobile at the right-side wall of the computational area, there are two floating supercontinents at the surface; this value may diminish or increase by 1.None of these configurations is stable. As an aspect ratio of the computational box is 10:1, its surface corresponds to the global Earth's model.

All supercontinents undergo the process of disintegration during their evolution.In our calculations,their lifetime is:in the case of a single supercontinent, consisting of the four continents with summary dimension of 2.7 D, 620 Ma as average; in the case of floating supercontinents of smaller dimensions(1.8-2.1 D),904 Ma on average. Here D is the whole mantle thickness (2850 km). The continent with a dimension of about 1 D (consisting of two small continents) experienced disintegration during the entire time of evolution only once (moreover, being at the side wall, i.e., being motionless for a long time).The rest of the time,it only joined and separated from other continents as a whole. On the contrary, the continent having a dimension of 0.7 D does not experience the breakup during total computation time at all. One reason for this difference is the inability of a small continent to create,owing to the heat-insulating effect, an upward mantle flow beneath itself. Such behavior agrees with the results presented by Bobrov et al.[51]for simple isoviscous convection, where the size of the continent decreasing from about 2 D to 0.7 D sharply reduced the effect of thermal insulation by a low-conductivity continent on the underlying mantle.

Another very important factor preventing the formation of ascending flows beneath a supercontinent is its displacement onto a new mantle area.Thus,a supercontinent may disintegrate either in the case that it remained almost motionless for a long time or it turns out to be above the structure of an already existing upwelling mantle flow beneath it, and two downward flows near its edges.

The lifetime of supercontinents in the model is possibly overestimated.The reason is that the breakup in the model occurs only along the already existing sutures (weak zones), which consist of the material of the oceanic lithosphere of reduced strength.Possibly,the division in other places of the supercontinent in reality could have occurred earlier, under favorable conditions.

4. Discussion and conclusions

We have used the 2D Cartesian model of thermochemical mantle convection with the non-Newtonian rheology and floating continents. Despite a relative simplicity of our model, the fundamental properties of mantle convection and supercontinent cycle were preserved. During the total calculation time, the supercontinent cycle was implemented several times. The number of individual continents on the surface as a function of time oscillates around 3.These results are in good agreement with the number of separated continents obtained for a spherical model [34].

The lifetime of supercontinents depends on their dimension and increases with decreasing of their size. According to our results,their lifetime is about 620 Ma for a large supercontinent(about 2.7 D), about 900 Ma for smaller dimensions (1.8-2.1 D), and several billion years for a micro-supercontinent with a dimension less than about 1 D. Excluding a single case of breakup, this small supercontinent only joined and separated from other continents as a whole. These values correlate with the results of Heron and Lowman[20],where the minimum width of a supercontinent at which its breakup can still occur was found to be about 1.5-1.6 D. On contrary, the continent having a dimension of 0.7 D does not experience any breakup during the whole period. The irregularity of the supercontinent cycle in our model is realized as in studies by Phillips and Bunge [34], Rolf et al. [22], and Bobrov and Baranov[27].

We find that the horizontal tensile(over-lithostatic)stresses σthroughout the supercontinent are about-200 to-250 MPa before the breakup of the supercontinent, while the subcontinental mantle with mantle plumes exhibits the horizontal compressive stresses of 50-100 MPa. Stresses generated in the supercontinent are large enough to break it up (see e.g., Gurnis [11]; Lowman and Jarvis[15];Yoshida,[39];Bobrov and Baranov[27]).The reason for this phenomenon is a sharp difference in the viscosities of the supercontinent and the underlying mantle materials (i.e., the highviscosity continental material cannot flow together with the underlying mantle). With the beginning of the supercontinent breakup, tensile stresses in the supercontinent change to compressive ones under the influence of actively intruding plume material into the opening rift. At stages of the assembly of continents,the opposite contrast of stresses take place.A supercontinent emerging as a result of a continental collision exhibits the horizontal compressive stresses(80-130 MPa),while a subcontinental mantle with plunging slab demonstrates tensile (relative to the average) stresses (-40 to -60 MPa). The boundary of a sharp change in values of σ, outlines the position of effective viscosity jump; i.e., this interface marks the base of the high-viscosity lithosphere. A similar phenomenon of a sharp spatial change in stresses also occurs in oceanic areas. For the occurrence of this phenomenon, differences in viscosities of 3 and even 2 orders of magnitude are sufficient. It should also be noted that for the real Earth, the sphericity effect reduces stresses by 20-30% [17]. However,recent spherical mantle models show rather high stresses(up to 100 MPa)at the base of continents [45].

Fig. 8. The number of continents as a function of time.

The calculated dynamic topography on the surface reflects the main features of the supercontinent cycle.At the initial stages,the bottom of the newly-opened ocean is dynamically raised by 2-4 km relative to old oceanic regions. When the continents collide, an orogen with a height of several kilometers appears between them. The same situation takes place along active continental margins.Here the height of rising orogens is about 2-4 km.During and immediately after the breakup, the supercontinent is dynamically uplifted by 1-2 km due to mantle plumes. Later, the heights of continents relative to the adjacent oceanic regions decrease on average by almost 1 km.Then,the heights of diverging continents continue to decrease gradually.When the continent is in the area of the mantle downwelling, it is lowered under its influence and actually does not differ in height from the neighboring oceanic areas.Note that the inherent height of continents above sea level should be added to these heights.

As a result, our model shows that during the supercontinental cycle, the dynamic topography of continents changes within the interval of about ±1.5 km relative to the mean level, which is in good accordance with real Earth values [52]. Previous mantle models, e.g., Lowman and Jarvis [15], provide the maximum amplitude of dynamic topography of about 0.5 km only.Such small value is explained by the fact that they used a constant viscosity in the mantle in their model.

Our results revealed that localized extreme values of vertical stresses σand topography are located in continents directly above subcontinental upward,and downward mantle flows,thus marking their positions.The continental relief above the uppermost part of the ascending stream is uplifted, while lowered over the downdrafts. Our modeling results thus demonstrate that a new expanding ocean must inevitably be shallow.

Our model adopted in this study took into account a set of features inherent in the real Earth. The considered stress fields supplemented the picture of commonly used fields of velocity,temperature, and viscosity. Computations were performed with a realistic Rayleigh number, which significantly changes the pattern of convection. Typical for lower Ra values, all fields become more unstable.Consequently,the feedback effect of continents on mantle convection weakens.As a result,the supercontinent cycle becomes irregular. Despite these findings, it is important to note that some other essential factors, such as fluids, should be taken into consideration to model the continental cycle more realistically.Moreover,the spherical model instead of 2D(planar)model will be used in our forthcoming study.

Author contribution

Alexander Bobrov provided a main framework of the research.

Alexey Baranov designed and performed calculations.

Robert Tenzer corrected the manuscript.

Con

fl

icts of interest

The authors declare that there is no conflicts of interest.

Acknowledgments

The authors would like to thank the anonymous reviewer and the editor, whose review 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.