(2020) 72:70
Miyagoshi et al. Earth, Planets and Space
https://doi.org/10.1186/s40623-020-01195-1

Open Access

EXPRESS LETTER

Tectonic plates in 3D mantle convection
model with stress‑history‑dependent rheology
Takehiro Miyagoshi1, Masanori Kameyama1,2*

and Masaki Ogawa1,3

Abstract
Plate tectonics is a key feature of the dynamics of the Earth’s mantle. By taking into account the stress-history-dependent rheology of mantle materials, we succeeded in realistically producing tectonic plates in our numerical model ofmantle convection in a three-dimensional rectangular box. The calculated lithosphere is separated into several pieces(tectonic plates) that rigidly move. Deformation of the lithosphere caused by the relative motion of adjacent plates isconcentrated in narrow bands (plate margins) where the viscosity is substantially reduced. The plate margins develop
when the stress exceeds a threshold and the lithosphere is ruptured. Once formed, the plate margins persist, evenafter the stress is reduced below the threshold, allowing the plates to stably move over geologic time. The verticalcomponent of vorticity takes a large value in the narrow plate margins. Secondary convection occurs beneath old
tectonic plates as two-dimensional rolls with their axes aligned to the direction of plate motion. The surface heat flow
decreases with increasing distance from divergent plate margins (ridges) in their vicinity in the way the cooling halfspace model predicts, but it tends towards a constant value away from ridges as observed for the Earth because of
the heat transport by the secondary convection.
Keywords: Mantle convection, 3D numerical model, Plate tectonics, Stress-history-dependent rheology
Introduction
One of the most challenging issues in studies of the
Earth’s mantle dynamics is to self-consistently produce
tectonic plates in three-dimensional numerical models
of mantle convection. The lithosphere develops along thesurface boundary, when the viscosity strongly depends
on temperature (Weertman 1970). The lithosphere produced in this way, however, behaves as a stagnant lid thatdevelops on top of the convecting mantle (Ogawa et al.
1991; Moresi and Solomatov 1995). To allow the lithosphere to be further divided into several pieces (tectonic
plates) that rigidly move, it is necessary to introduce platemargins, where the mechanical strength is much lower
than that of the surrounding tectonic plates in the lithosphere. In the Earth, plate margins develop when a stresshigher than the rupture strength of the lithosphere is*Correspondence: kameyama@sci.ehime‑u.ac.jp
2
Geodynamics Research Center (GRC), Ehime University, 2‑5 Bunkyo‑cho,Matsuyama, Ehime 790‑8577, Japan
Full list of author information is available at the end of the article

induced there by, for example, Large Igneous Provinces(Coffin and Eldholm 1994). Once formed, plate margins
persist, even after the stress is reduced below the rupture strength. Indeed, various geophysical and tectonic
observations (Kanamori 1980; Zhong and Watts 2013;
Gao and Wang 2014) show that mechanically weak platemargins persist in the lithosphere at stress levels lower
than those for mechanically strong plate interiors. This
stress-history-dependent behaviour of the lithosphere isthought to play a crucial role in plate tectonics (Bercovici
1998; Zhong et al. 1998; Ogawa 2003). Here, we present a
three-dimensional model of mantle convection to examine how a stress-history-dependent rheology exerts control over the dynamics of the lithosphere.

Methods
The model of mantle convection with tectonic plates
was calculated with the ACuTEMan code (Kameyama
et al. 2005; Kameyama 2005) as thermal convection of an
incompressible Newtonian fluid with an infinite Prandtl
number in an internally heated three-dimensional© The Author(s) 2020. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) andthe source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material
in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material
is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds thepermitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat​iveco​
mmons​.org/licen​ses/by/4.0/.

Miyagoshi et al. Earth, Planets and Space(2020) 72:70

rectangular box with an aspect ratio of 4 × 4 and a heightof 3000 km. The temperature on the surface boundary isfixed at 0◦ C, whereas the other boundaries are insulating.
All of the boundaries are impermeable and shear-stressfree. The viscosity strongly depends on the temperature,and the lithosphere develops as the coldest and most viscous part of the cold thermal boundary layer along thesurface. The viscosity is also a two-valued function ofstress in a range from σm (the strength of mechanical coupling at the plate margins) to σp (the rupture strength ofthe lithosphere) (Ogawa 2003). The viscosity takes a high(or “intact”) value for the plate interior when the stress issufficiently low, but it drops to a low (or “damaged”) value
for the plate margins as the stress increases and exceeds
σp; the viscosity remains low, even when the stress is
reduced below σp as long as it is higher than σm (see theAdditional file 1 for more detail). Which of the two values the viscosity takes in the stress range is determinedby whether or not the material has experienced a stresshigher than σp in the past, and this memory of stresshistory persists for an indefinitely long period of time
as long as σ stays in the range from σm to σp. We stored
this information of stress history, which is transported
by convection, in a scalar quantity called damage parameter (Bercovici et al. 2001; Ogawa 2003). The persistentmemory of stress-history of our model is different fromthe memory assumed in some earlier models (e.g., Bercovici 1998; Foley et al. 2014) that fades away with a decay
time of several hundred million years. This persistentmemory is the reason why a more plate-like behaviour ofthe lithosphere is produced in our model. The viscosity
for plate interiors is ηp = 1027 Pa s, whereas that for platemargins is 3.2 × 1023 Pa s along the surface boundary. We
also assume σp = 200 MPa, a value somewhat lower thanearlier estimates (Kohlstedt et al. 1995; Watts and Burov
2003; Jain et al. 2017), and σm = 67 MPa.

Results
The lithosphere simulated in our model is rifted into
several highly viscous pieces (tectonic plates) separated
by narrow plate margins where the viscosity is substantially lower than that in the plate interiors (Fig. 1a).
Each of the plates rigidly moves (see the arrows forvelocity vector V in Fig. 1a), and the deformation of thelithosphere caused by relative motion between adjacentplates is accommodated in the narrow plate margins.
The cold plates sink into the mantle at convergent plate∂Vy∂Vxmargins (subduction zones) where
+
<0
∂x
∂y
holds (see the temperature distribution on the verticalsection presented in Fig. 1a). The vertical component of
the vorticity field (∇ × V )z takes a large value in theplate margins except in the subduction zones (Fig. 1b);Page 2 of 6abFig. 1 The overall three-dimensional structure of thermal convection
obtained in our numerical model with stress-history-dependentrheology. a The distribution of temperature T in a vertical crosssection and that of viscosity η (colour) and velocity V (arrows) on
the surface. The value of η is normalised by ηp, the viscosity of intact
material at the surface. b The distribution of the vertical componentof vorticity (∇ × V)z on the surface(∇ × V )z is somewhat large over the entire small plate
along the west side in Fig. 1b too, implying that thisplate rigidly rotates. As far as the authors know, Fig. 1b
is the first example where (∇ × V )z takes a significant
value in narrow plate margins in 3D numerical modelsof thermal convection where tectonic plates spontaneously develop (Tackley 2000; Bercovici et al. 2015).
The plate motion induces lateral variation in the surface heat flow (Fig. 2a). The heat flow decreases withincreasing distance L from the divergent plate margins(ridges) in their vicinity. The almost linear dependenceof 1/(heat flow)2 on L observed in Fig. 2b is consistent
with a prediction from the cooling half-space model oftectonic plates (Sclater et al. 1980) because L is almost
proportional to the surface age because of the calculated stable plate motion in our models (see below).

Miyagoshi et al. Earth, Planets and Spacea(2020) 72:70

Page 3 of 6bcFig. 2 The lateral variations in the surface heat flow and convecting flow structure. a The distribution of the inverse square of the surface heat flow
shown in both colour and white contour lines. b The plot of the inverse square of heat flow at the surface as a function of the distance L from thedivergent plate boundary (ridge) along a horizontal track indicated by the dashed line in a. Also shown by the dash-dotted line is the theoreticalprediction from the half-space cooling for small L, while by the dashed line is the deviation from it for large L. c The distribution of verticalcomponent of velocity Vz on the horizontal plane at the height of z = 0.672

Away from ridges, the value of 1/(heat flow)2 deviates from the linear dependence and tends towards aconstant value (Fig. 2b). Figures 2a and c show that thedeviation is due to heat transport by the secondary convection that occurs beneath old plates: The distributionof the vertical component of velocity on the horizontal
plane at z = 0.672 plotted in Fig. 2c shows that the secondary convection indeed occurs in the form of twodimensional rolls whose axes are aligned in the directionof plate motion. The convection rolls extend to the depthof about 1000 km. The secondary convection also induces
a stripe pattern in the surface heat flow (Fig. 2a). The
onset and pattern of secondary convection beneath oldoceans are consistent with a prediction from earlier laboratory (Parsons and McKenzie 1978; Richter and Parsons

1975) and numerical (van Hunen et al. 2003; Huang andZhong 2005) experiments of thermal convection beneathmoving plates. Such stripe patterns may be identified in
the future by high-resolution tomographic studies of theupper mantle beneath old parts of large oceanic plates
such as the Pacific and Indian plates.
The time evolution of the velocity distribution along
the surface presented in Fig. 3 shows that the tectonic
plates produced in our model move rather steadily on
the non-dimensional time scale of 0.01 (or about 3 billion years). Plate A, for example, moves to the southeast, whereas plate B moves to the north, and the ridge
between them exists throughout the presented time span.
The velocities of these plates do not change so wildly,
although the rigid rotation of plate A observed in Fig. 3aMiyagoshi et al. Earth, Planets and Space(2020) 72:70

Page 4 of 6

a

b

c

d

Fig. 3 The time evolution of the surface plate motions. Snapshots of the distributions of viscosity η (colour) and velocity V (arrows) in a horizontal
plane along the surface at the elapsed times indicated in the figure. The values of η are normalised by ηp, the viscosity of intact material at thesurfacestops by t = 5.328 × 10−2. The plate margins move inaccordance with the plate motion. Ridge R (Fig. 3b), forexample, moves to the south in frames (b) to (d) because
the southward motion of plate C is much faster than thenorthward motion of plate B.
Plate margins persist in our calculation because thelithosphere remembers the location of existing plate margins owing to the assumed stress-history dependenceof the viscosity. The margins do not spontaneously disappear; they disappear only when they merge with eachother or with the sidewalls of the convecting box owing
to the advection by tectonic plates.
An interesting feature of the modelled tectonic plates
shown in Fig. 3 is that the ridges are just simple narrow
bands but the subduction zones are more diffuse. Besides,
fragments of plates (micro-plates) are often formed along
subduction zones (see the area around the southwest andsoutheast corners of the convecting box). This featurearises because high stress tends to be induced aroundsubduction zones by plate motion. This feature is alsoMiyagoshi et al. Earth, Planets and Space(2020) 72:70

observed for the Earth (Gordon et al. 1998). In Fig. 3 wecould not find a clear relationship between formationof micro-plates and the curvature of plate margins, as
opposed to the suggestion by Mallard et al. (2016).

Discussion
Our model of tectonic plates shown in Figs. 1 to 3 indicates that a stress-history-dependent rheology that allows
the lithosphere to “remember” the location of plate margins is crucial for plate tectonics to stably operate overgeologic time (Zhong et al. 1998; Gurnis et al. 2000). It
is difficult to produce stable plate motion together witha large value of (∇ × V )z in narrow plate margins innumerical models where the viscosity depends only on
the instantaneous stress (Moresi and Solomatov 1998;
Tackley 2000). An example is the yielding model; the lithosphere is intact and its viscosity is high when the stressis less than a threshold, called the yield strength, but it isruptured and its viscosity is reduced as the stress reaches
the threshold. In the yielding model, however, plate margins often spontaneously disappear or migrate at a velocity much higher than the plate velocity as the stress fieldevolves with time in the lithosphere. Stable motion of
tectonic plates separated by narrow plate margins is difficult to produce, even in some models where the viscosity
in the lithosphere does depend on its stress history overthe past several hundred million years (Foley et al. 2014).
See the Additional file 1 for more discussion.
The stable plate motion produced in our model isimportant for understanding the global-scale heterogeneity of the Earth’s mantle revealed by tomographic
studies. The large low shear velocity provinces (LLSVPs)
develop on the core-mantle boundary away from the slab
graveyards where subducted slabs have accumulated overthe past several hundred million years (Richards andEngebretson 1992; Cottaar and Lekic 2016). The samefeature is observed for numerical models of the evolutionof thermo-chemical piles on the core-mantle boundary
where the Earth’s plate motion in the past is given as themechanical boundary condition on the surface (Zhang
et al. 2010). If plate margins jump and/or migrate sowildly and plate velocity varies with time as strongly as
observed in some numerical models where the rheology
depends on instantaneous stress rather than stress-history, regions such as LLSVPs and slab graveyards would
not have developed (Nakagawa and Tackley 2014). Stable
plate motion produced in our model opens a way to selfconsistent numerical modelling of the development andevolution of large-scale structures in the Earth’s mantle
throughout its history. It is important to further extend
our model to three-dimensional spherical geometry of
the mantle in future exploration.

Page 5 of 6

Conclusions
In this study we have realistically produced for the firsttime tectonic plates in our model of thermal convection of three-dimensional mantle, by assuming a stresshistory-dependent rheology for mantle materials. The
tectonic plates calculated in our model show the following five features: (i) rigid motion of cold and stiff
plates (both translation and rotation), (ii) accommodation of relative motion of adjacent plates in narrowplate margins, (iii) stable tectonic motions over several
billion years, (iv) decrease in surface heat flow in thevicinity of divergent plate boundaries in a way consistent with the prediction from the half-space coolingmodel, and (v) deviation of the surface heat flow fromthe half-space cooling model away from divergent plateboundaries. We believe that our model opens a way to
self-consistent numerical modelling of development
and evolution of large-scale structure in the Earth’s
mantle in its history.
Supplementary information
Supplementary information accompanies this paper at https​://doi.
org/10.1186/s4062​3-020-01195​-1.
Additional file 1. Additional Materials for “Tectonic plates in 3D mantle
convectionmodel with stress-history-dependent rheology”.
Abbreviation
LLSVPs: The large low shear velocity provinces.
Acknowledgements
This research in part used computational resources of the K computer provided by the RIKEN Center for Computational Science through the HPCI System Research project (Project ID: hp160254, hp170225, hp180199, hp190170),
and those of the Earth Simulator at Japan Agency for Marine-Earth Scienceand Technology.
Authors’ contributions
T.M. performed numerical simulations, analysis of simulation data, andprepared the manuscript. M.K. developed the numerical simulation code
ACuTEMan and prepared the manuscript. M.O. prepared the manuscript. All
authors discussed numerical simulation results. All authors read and approved
the final manuscript.
Funding
This work was supported by JSPS KAKENHI Grant Numbers JP25287110,
JP15H05834, 18K03724, a joint research program at Geodynamics ResearchCenter, Ehime University, and “Exploratory Challenge on Post-K computer”(Elucidation of the Birth of Exoplanets [Second Earth] and the Environmental
Variations of Planets in the Solar System).
Availability of data and materials
The data that support the findings of this study are available upon reasonable
request to T.M. (email: miyagoshi@jamstec.go.jp).
The numerical code ACuTEMan is available upon request to M.K. (email:kameyama@sci.ehime-u.ac.jp).
Competing interests
The authors declare no competing interests.

Miyagoshi et al. Earth, Planets and Space(2020) 72:70

Author details
1
Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 3173‑25
Showa‑machi, Kanazawa‑ku, Yokohama 236‑0001, Japan. 2 Geodynamics
Research Center (GRC), Ehime University, 2‑5 Bunkyo‑cho, Matsuyama, Ehime
790‑8577, Japan. 3 Department of Earth Sciences and Astronomy, Universityof Tokyo at Komaba, 3‑8‑1 Komaba, Meguro, Tokyo 153‑8902, Japan.
Received: 23 March 2020 Accepted: 12 May 2020

References
Bercovici D (1998) Generation of plate tectonics from lithosphere-mantle flow
and void-volatile self-lubrication. Earth Planet Sci Lett 154(1–4):139–151.
https​://doi.org/10.1016/S0012​-821X(97)00182​-9
Bercovici D, Ricard Y, Schubert G (2001) A two-phase model for compaction
and damage 3 applications to shear localization and plate boundary formation. J Geophys Res 106(B5):8925–8939. https​://doi.org/10.1029/2000J​
B9004​32
Bercovici D, Tackley PJ, Ricard Y (2015) The generation of plate tectonics frommantle dynamics. In: Schubert G (ed) Treatise on Geophysics, vol 7, 2ndedn, Elsevier, https​://doi.org/10.1016/B978-0-444-53802​-4.00135​-4
Coffin MF, Eldholm O (1994) Large igneous provinces: crustal structure, dimensions, and external consequences. Rev Geophys 32(1):1–36. https​://doi.
org/10.1029/93RG0​2508
Cottaar S, Lekic V (2016) Morphology of seismically slow lower-mantle
structures. Geophys J Int 207(2):1122–1136. https​://doi.org/10.1093/gji/
ggw32​4
Foley BJ, Bercovici D, Elkins-Tanton LT (2014) Initiation of plate tectonics frompost-magma ocean thermochemical convection. J Geophys Res Solid
Earth 119(11):8538–8561. https​://doi.org/10.1002/2014J​B0111​21
Gao X, Wang K (2014) Strength of stick-slip and creeping subduction megathrusts from heat flow observations. Science 345(6200):1038–1041. https​://doi.org/10.1126/scien​ce.12554​87
Gordon RG, DeMets C, Royer JY (1998) Evidence for long-term diffuse
deformation of the lithosphere of the equatorial Indian Ocean. Nature
395(6700):370–374. https​://doi.org/10.1038/26463​
Gurnis M, Zhong S, Toth J (2000) On the competing roles of fault reactivation
and brittle failure in generating plate tectonics from mantle convection. In: Richards MA, Gordon RG, van der Hilst RD, eds. The History andDynamics of Global Plate Motions, Geophysical Monograph, vol 121,
American Geophysical Union, pp 73–94, https​://doi.org/10.1029/GM121​
p0073​
Huang J, Zhong S (2005) Sublithospheric small-scale convection and its implications for the residual topography at old ocean basins and the platemodel. J Geophys Res 110:B05404. https​://doi.org/10.1029/2004J​B0031​53
Jain C, Korenaga J, Karato S (2017) On the yield strength of oceanic
lithosphere. Geophys Res Lett 44(19):9716–9722. https​://doi.
org/10.1002/2017G​L0750​43
Kameyama M (2005) ACuTEMan: A multigrid-based mantle convection simulation code and its optimization to the Earth Simulator. J Earth Simulator
4:2–10Kameyama M, Kageyama A, Sato T (2005) Multigrid iterative algorithm using
pseudo-compressibility for three-dimensional mantle convection withstrongly variable viscosity. J Comput Phys 206(1):162–181. https​://doi.
org/10.1016/j.jcp.2004.11.030
Kanamori H (1980) The state of stress in the Earth’s lithosphere. In: Dziewonski
A, Boschi E (eds) Physics of the earth’s interior. North-Holland, Amsterdam, pp 531–554
Kohlstedt DL, Evans B, Mackwell SJ (1995) Strength of the lithosphere:
constraints imposed by laboratory experiments. J Geophys Res
100(B9):17587–17602. https​://doi.org/10.1029/95JB0​1460

Page 6 of 6

Mallard C, Coltice N, Seton M, Müller RD, Tackley PJ (2016) Subduction controls
the distribution and fragmentation of Earth’s tectonic plates. Nature
535(7610):140–143. https​://doi.org/10.1038/natur​e1799​2
Moresi LN, Solomatov V (1995) Numerical investigation of 2D convection withextremely large viscosity variations. Phys Fluids 7(9):2154–2162. https​://doi.org/10.1063/1.86846​5
Moresi LN, Solomatov V (1998) Mantle convection with a brittle lithosphere:
thoughts on the global tectonic styles of the Earth and Venus. Geophys J
Int 133(3):669–682. https​://doi.org/10.1046/j.1365-246X.1998.00521​.x
Nakagawa T, Tackley PJ (2014) Influence of combined primordial layeringand recycled MORB on the coupled thermal evolution of Earth’s mantle
and core. Geochem Geophys Geosyst 15(3):619–633. https​://doi.
org/10.1002/2013G​C0051​28
Ogawa M (2003) Plate-like regime of a numerically modeled thermal convection in a fluid with temperature-, pressure-, and stress-history-dependent
viscosity. J Geophys Res 108(B2):2067. https​://doi.org/10.1029/2000J​
B0000​69
Ogawa M, Schubert G, Zebib A (1991) Numerical simulations of threedimensional thermal convection in a fluid with strongly temperaturedependent viscosity. J Fluid Mech 233:299–328. https​://doi.org/10.1017/
S0022​11209​10004​96
Parsons B, McKenzie D (1978) Mantle convection and the thermal structureof the plates. J Geophys Res 83(B9):4485–4496. https​://doi.org/10.1029/
JB083​iB09p​04485​
Richards MA, Engebretson DC (1992) Large-scale mantle convection andthe history of subduction. Nature 355(6359):437–440. https​://doi.
org/10.1038/35543​7a0
Richter FM, Parsons B (1975) On the interaction of two scales of convection in
the mantle. J Geophys Res 80(17):2529–2541. https​://doi.org/10.1029/
JB080​i017p​02529​
Sclater JG, Jaupart C, Galson D (1980) The heat flow through oceanic and continental crust and the heat loss of the Earth. Rev Geophys 18(1):269–311.
https​://doi.org/10.1029/RG018​i001p​00269​
Tackley PJ (2000) Self-consistent generation of tectonic plates in timedependent, three-dimensional mantle convection simulations 2 strain
weakening and asthenosphere. Geochem Geophys Geosyst 1(8):1026.
https​://doi.org/10.1029/2000G​C0000​43
van Hunen J, Huang J, Zhong S (2003) The effect of shearing on the onset andvigor of small-scale convection in a Newtonian rheology. Geophys Res
Lett 30(19):1991. https​://doi.org/10.1029/2003G​L0181​01
Watts AB, Burov EB (2003) Lithospheric strength and its relationship to theelastic and seismogenic layer thickness. Earth Planet Sci Lett 213(1–2):113–131. https​://doi.org/10.1016/S0012​-821X(03)00289​-9
Weertman J (1970) The creep strength of the Earth’s mantle. Rev Geophys
8(1):145–168. https​://doi.org/10.1029/RG008​i001p​00145​
Zhang N, Zhong S, Leng W, Li ZX (2010) A model for the evolution of the
Earth’s mantle structure since the Early Paleozoic. J Geophys Res
115:B06401. https​://doi.org/10.1029/2009J​B0068​96
Zhong S, Watts AB (2013) Lithospheric deformation induced by loading of the
Hawaiian Islands and its implications for mantle rheology. J Geophys Res
118(11):6025–6048. https​://doi.org/10.1002/2013J​B0104​08
Zhong S, Gurnis M, Moresi L (1998) Role of faults, nonlinear rheology, andviscosity structure in generating plates from instantaneous mantle flow models. J Geophys Res 103(B7):15525–15268. https​://doi.
org/10.1029/98JB0​0605

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.