%%
%% This is a LaTeX template file for EPS
%%
%% modified by M. M. on 17 March 2014
%% modified by Y.O. on 7 September 2014
%% modified by Y.O on 12 December, 2018
%% modified by Y.O and M. M. on 6 December, 2019

% \documentclass[submit]{eps}
\documentclass{EPS}

\usepackage{lineno}
\linenumbers*[1]

% MCK needs
\usepackage{amsmath,bm}
%

\setcounter{page}{1}
\title{Tectonic plates in 3D mantle convection model with
stress-history-dependent rheology}
\author{%
Takehiro Miyagoshi,
 Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 3173-25
 Showa-machi, Kanazawa-ku, Yokohama 236-0001, Japan,
 miyagoshi@jamstec.go.jp\\
%
 Masanori Kameyama,
 Geodynamics Research Center (GRC), Ehime University, 2-5 Bunkyo-cho,
 Matsuyama, Ehime 790-8577, Japan,
 kameyama@sci.ehime-u.ac.jp\\
 %
 Masaki Ogawa,
 Department of Earth Sciences and Astronomy, University of Tokyo at Komaba, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan,
 cmaogawa@mail.ecc.u-tokyo.ac.jp
 }

\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 of mantle 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 is concentrated 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, even after the stress is
 reduced below the threshold, allowing the plates to stably move over
 geologic time.
 The vertical component 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
 half-space 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}

\begin{document}

\maketitle

\section{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 the surface boundary, when the viscosity
strongly depends on temperature \citep{weertman70}.
The lithosphere produced in this way, however, behaves as a stagnant lid
that develops on top of the convecting mantle \citep{ogawa91,moresi95a}.
To allow the lithosphere to be further divided into several pieces
(tectonic plates) that rigidly move, it is necessary to introduce plate
margins, 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 stress higher than the
rupture strength of the lithosphere is induced there by, for example,
Large Igneous Provinces \citep{coffin94}.
Once formed, plate margins persist, even after the stress is reduced
below the rupture strength.
Indeed, various geophysical and tectonic observations
\citep{kanamori80sos,zhong13jgrb,xgao14sci} show that mechanically weak
plate margins persist in the lithosphere at stress levels lower than
those for mechanically strong plate interiors.
This stress-history-dependent behaviour of the lithosphere is thought to
play a crucial role in plate tectonics
\citep{bercovici98a,zhong98,ogawa03a}.
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.

\section{Methods}

The model of mantle convection with tectonic plates was calculated with
the ACuTEMan code \citep{kameyama05jcp,kameyama05jes} as thermal
convection of an incompressible Newtonian fluid with an infinite Prandtl
number in an internally heated three-dimensional rectangular box with an
aspect ratio of $4\times4$ and a height of $3000$ km.
The temperature on the surface boundary is fixed at $0^\circ$~C, whereas
the other boundaries are insulating.
All of the boundaries are impermeable and shear-stress free.
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 the surface.
The viscosity is also a two-valued function of stress in a range from
$\sigma_m$ (the strength of mechanical coupling at the plate margins) to
$\sigma_p$ (the rupture strength of the lithosphere) \citep{ogawa03a}.
The viscosity takes a high value for the plate interior when the stress
is sufficiently low, but it decreases for the plate margins as the
stress increases and exceeds $\sigma_p$;
the viscosity remains low, even when the stress is reduced below
$\sigma_p$ as long as it is higher than $\sigma_m$.
Which of the two values the viscosity takes in the stress range is
determined by whether or not the material has experienced a stress
higher than $\sigma_p$ in the past.
We stored this information of stress history, which is transported by
convection, in a scalar quantity called damage parameter
\citep{bercovici01jgrb2,ogawa03a};
see the Additional Materials for more details.
The viscosity for plate interiors is $\eta_p=10^{27}$ \mbox{Pa\,s},
whereas that for plate margins is $3.2\times10^{23}$ \mbox{Pa\,s} along
the surface boundary.
We also assume $\sigma_p=200$ \mbox{MPa}, a value somewhat lower than
earlier estimates \citep{kohlstedt95,watts03epsl,cjain17grl}, and
$\sigma_m=67$ \mbox{MPa}.

\section{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
(Figure \ref{zu1}a).
Each of the plates rigidly moves (see the arrows for velocity vector
$\bm{V}$ in Figure \ref{zu1}a), and the deformation of the lithosphere
caused by relative motion between adjacent plates is accommodated in the
narrow plate margins.
The cold plates sink into the mantle at convergent plate margins
(subduction zones) where
$\displaystyle\frac{\partial{V_x}}{\partial{x}}+\frac{\partial{V_y}}{\partial{y}}<0$
holds (see the temperature distribution on the vertical section
presented in Figure \ref{zu1}a).
The vertical component of the vorticity field
$\left(\bm{\nabla}\times\bm{V}\right)_z$ takes a large value in the
plate margins except in the subduction zones (Figure \ref{zu1}b);
$\left(\bm{\nabla}\times\bm{V}\right)_z$ is somewhat large over the
entire small plate along the west side in Figure \ref{zu1}b too,
implying that this plate rigidly rotates.
As far as the authors know, Figure \ref{zu1}b is the first example where
$(\bm{\nabla}\times\bm{V})_z$ takes a significant value in narrow plate
margins in 3D numerical models of thermal convection where tectonic
plates spontaneously develop \citep{tackley00c,bercovici15treatise}.

The plate motion induces lateral variation in the surface heat flow
(Figure \ref{zu2}a).
The heat flow decreases with increasing distance $L$ from the divergent
plate margins (ridges) in their vicinity.
The almost linear dependence of $1/(\text{heat~flow})^2$ on $L$ observed
in Figure \ref{zu2}b is consistent with a prediction from the cooling
half-space model of tectonic plates \citep{sclater80rg} because $L$ is
almost proportional to the surface age because of the calculated stable
plate motion in our models (see below).

Away from ridges, the value of $1/(\text{heat~flow})^2$ deviates from
the linear dependence and tends towards a constant value (Figure
\ref{zu2}b).
Figures \ref{zu2}a and \ref{zu2}c show that the deviation is due to heat
transport by the secondary convection that occurs beneath old plates:
The distribution of the vertical component of velocity on the horizontal
plane at $z=0.672$ plotted in Figure \ref{zu2}c shows that the secondary
convection indeed occurs in the form of two-dimensional rolls whose axes
are aligned in the direction of plate motion.
The convection rolls extend to the depth of about $1000$ km.
The secondary convection also induces a stripe pattern in the surface
heat flow (Figure \ref{zu2}a).
This pattern of secondary convection is consistent with a prediction
from earlier laboratory experiments of thermal convection beneath moving
plates \citep{parsons78jgrb,richter75jgrb}.
Such stripe patterns may be identified in the future by high-resolution
tomographic studies of the upper 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 Figure \ref{zu3} 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 Figure \ref{zu3}a stops by
$t=5.328\times10^{-2}$.
The plate margins move in accordance with the plate motion.
Ridge R (Figure \ref{zu3}b), for example, moves to the south in frames
(b) to (d) because the southward motion of plate C is much faster than
the northward motion of plate B.

Plate margins persist in our calculation because the lithosphere
remembers the location of existing plate margins owing to the assumed
stress-history dependence of the viscosity.
The margins do not spontaneously disappear;
they disappear only when they merge with each other 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 Figure
\ref{zu3} 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 and southeast
corners of the convecting box).
This feature arises because high stress tends to be induced around
subduction zones by plate motion.
This feature is also observed for the Earth \citep{gordon98nature}.
In Figure \ref{zu3} we could not find a clear relationship between
formation of micro-plates and the curvature of plate margins, as opposed
to the suggestion by \cite{mallard16nature}.

\section{Discussion}

Our model of tectonic plates shown in Figures \ref{zu1} to \ref{zu3}
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 over geologic time
\citep{zhong98,gurnis00gm}.
It is difficult to produce stable plate motion together with a large
value of $(\bm{\nabla}\times\bm{V})_z$ in narrow plate margins in
numerical models where the viscosity depends only on the instantaneous
stress \citep{moresi98,tackley00c}.
An example is the yielding model;
the lithosphere is intact and its viscosity is high when the stress is
less than a threshold, called the yield strength, but it is ruptured 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 field evolves 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 over the past several
hundred million years \citep{foley14jgrb}.
See the Additional Materials for more discussion.

The stable plate motion produced in our model is important 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 over the past several hundred million years
\citep{richards92nature,cottaar16gji}.
The same feature is observed for numerical models of the evolution of
thermo-chemical piles on the core-mantle boundary where the Earth's
plate motion in the past is given as the mechanical boundary condition
on the surface \citep{nzhang10jgrb}.
If plate margins jump and/or migrate so wildly 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 \citep{nakagawa14gc}.
Stable plate motion produced in our model opens a way to
self-consistent numerical modelling of the development and evolution 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.

\section{Conclusions}

 In this study we have realistically reproduced for the first time
 tectonic plates in our model of thermal convection of three-dimensional
 mantle, by assuming a stress-history-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 narrow
 plate margins,
 (iii) stable tectonic motions over several billion years,
 (iv) decrease in surface heat flow in the vicinity of divergent plate
 boundaries in a way consistent with the prediction from the half-space
 cooling model,
 and
 (v) deviation of the surface heat flow from the half-space cooling
 model away from divergent plate boundaries.
 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, and hence draws attention from broad
 disciplines in geoscience, ranging from geodynamics, geophysics to
 studies of Earth's evolution.

\section{List of abbreviations}
The large low shear velocity provinces (LLSVPs)

\section{Declarations}

\section{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).

\section{Competing interests}

The authors declare no competing interests.

\section{Funding}
This work was supported by
JSPS KAKENHI Grant Numbers JP25287110, JP15H05834, 18K03724,
a joint research program at Geodynamics Research Center, 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).

\section{Authors' contributions}

 T.M.\ performed numerical simulations,
 analysis of simulation data, and prepared 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.

\acknowledgments{%
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).
%https://www.hpci-office.jp/pages/procedure_achievements?tab=publication
%HPCI 「京」課題 hp160254, hp170225, hp180199, hp190170
}

%
\bibliographystyle{spbasic}
\bibliography{manuscript4eps}

\begin{figure}[p]
 \centering
 %\includegraphics[clip]{zu1-crop.pdf}
 \caption{The overall three-dimensional structure of thermal convection
 obtained in our numerical model with stress-history-dependent rheology.
 (a) The distribution of temperature $T$ in a vertical cross section and
 that of viscosity $\eta$ (colour) and velocity $\bm{V}$ (arrows) on the
 surface.
 The value of $\eta$ is normalised by $\eta_p$, the viscosity of intact
 material at the surface.
 (b) The distribution of the vertical component of vorticity
 $(\bm{\nabla}\times\bm{V})_z$ on the surface.
 }
 \label{zu1}
\end{figure}

\begin{figure}[p]
 \centering
 %\includegraphics[clip]{zu2-crop.pdf}
 \caption{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 the divergent plate boundary (ridge)
 along a horizontal track indicated by the dashed line in (a).
 (c) The distribution of vertical component of velocity $V_z$ on the
 horizontal plane at the height of $z=0.672$.
 }
 \label{zu2}
\end{figure}

\begin{figure}[p]
 \centering
 %\includegraphics[clip]{zu3-crop.pdf}
 \caption{The time evolution of the surface plate motions.
 Snapshots of the distributions of viscosity $\eta$ (colour) and
 velocity $\bm{V}$ (arrows) in a horizontal plane along the surface at
 the elapsed times indicated in the figure.
 The values of $\eta$ are normalised by $\eta_p$, the viscosity of
 intact material at the surface.
 }
 \label{zu3}
\end{figure}

\end{document}
