Additional Materials for “Tectonic plates in 3D mantle convection model with stress-history-dependent rheology”
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 Bunkyocho, Matsuyama, Ehime 790-8577, Japan, kameyama@sci.ehime-u.ac.jp
Masaki Ogawa, Department of Earth Sciences and Astronomy, University of Tokyo atKomaba, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan, cmaogawa@mail.ecc.u-tokyo.ac.jp

Preprint submitted to Earth, Planets and Space on March 23, 20201

Details of numerical model
Thermal convection is calculated in a rectangular box with an aspect ratio 4×4 by ACuTEMan numericalcode under the Boussinesq approximation (Kameyama et al. 2005; Kameyama 2005). The convecting
ﬂuid is incompressible and Newtonian, and the Prandtl number is inﬁnite. The box is internally heated,
and the internal heating rate is 12.8 in a dimensionless unit (equivalent to the surface heat ﬂow of about
23 mW m−2 ). The viscosity depends on temperature T , height z measured from the bottom boundary,and “degree of damage” ω, discussed below as[]Fωη = exp E(Tref − T ) + V (1 − z) −1+ω(1)(see Equation (1) of our earlier paper by Ogawa (2003)), where Tref , E, V , and F are constants. InEquation (1) and below, quantities without asterisks are non-dimensional; the length is normalised by
the height of the box d∗ = 3000 km, the temperature by the temperature scale ∆T ∗ = 1650◦ C, andthe viscosity by η0∗ = 2.4 × 1020 Pa s, a typical value for the uppermost mantle. We assume that theRayleigh number in the momentum equation Ra =

ρ∗0 α∗ ∆T ∗ g ∗ d∗ 3is 1.8 × 108 , where ρ∗0 = 3300 kg m−3η0∗ κ∗

is the reference density, α∗ = 3 × 10−5 K−1 is the thermal expansivity, κ∗ = 10−6 m2 s−1 is the thermal
diﬀusivity, and g ∗ = 10 m s−2 is the gravitational acceleration. The degree of damage ω evolves as(Bercovici 1998; Bercovici et al. 2001)∂ω+ V · ∇ω = Γσij ε̇ij − λω + ν∇2 ω,∂t(2)where Γ is a constant, σij is stress, ε̇ij is strain rate, andλ = λ0 exp [E(T − Tsrfc )] .

(3)

(See Equations (2) and (3) of Ogawa (2003).) Here, Tsrfc is the temperature at the surface whose∗
dimensional value is Tsrfc= 0◦ C, and the time t∗ is normalised by the thermal diﬀusion time d∗ 2 /κ∗ . An

artiﬁcial diﬀusion term ν∇2 ω is added on the right-hand side of Equation (2) to avoid spurious ﬂuctuationsin the ω-distributions that cause too wild variations of the viscosity in the lithosphere. In this study, thedimensionless value of ν is chosen to be ν = c(δxmin )2 /δt, where δxmin is the minimum mesh spacing, δt
is a time increment which is determined from the numerical stability of heat-advection-calculation, andc is a constant on the order of 10−1 .
The assumed parameter values are E = 19.3, V = 4.6, F = 8.06, Γ = 4 × 10−4 , λ0 = 105 , and∗
Tref= 1300◦ C, which are motivated by the plate-like regime of an earlier two-dimensional model (Ogawa22003). The calculation is started from the initial state taken from the ﬁnal state of the weak plate
regime shown in Figure A2, and is continued over the period of dimensionless time on the order of 10−2 .
The continuity and momentum equations are iteratively solved for the velocity V and pressure at eachtimestep, until the L2 -norm of residual becomes smaller than that of driving forces by a factor of 10−6 .

Notes on stress-history-dependence in rheology in our model
The formulation (1)-(3) is not a model of a speciﬁc physical process taking place in the Earth’s lithosphericmaterials, despite its similarity to that of void-volatile mechanism (Bercovici 1998; Bercovici et al.
2001). We employed this formulation together with the values of parameters described above so as
to phenomenologically introduce a stress-history-dependence into rheology, which is deﬁnitely importantin understanding the Earth’s plate motion (Zhong et al. 1998; Tackley 2000; Bercovici et al. 2015) in theway illustrated in Figure A1.
Let us consider the steady solution ωs of (2) calculated from
0 = Γσij ε̇ij − λωs = ΓλωsσII 2− λωs = −ηη(η−ΓσII 2λωs),(4)where σII is the second invariant of the stress tensor. The large value of λ0 adopted here implies that
ω is always kept very close to ωs . When F > 4, the steady solution ωs depends on σII , as illustrated inFigure A1a. Here,[σm,p = σref√where σref ≡F−1±2√

F2−F4] 12[Fexp −4(1±√41−F)],(5)

λ0 ηcoldand ηcold = exp [E(Tref − Tsrfc ) + V (1 − z)] is the viscosity at T = Tsrfc . (In theΓmain text, ηp stands for ηcold at z = 1.) The plus/minus sign takes + for σm and − for σp . In thelithosphere, σp corresponds to the rupture strength of tectonic plates, whereas σm corresponds to thestrength of the mechanical coupling at plate margins. As illustrated in Figure A1, ωs is a two-valued
function of stress in the range of σm ≤ σII ≤ σp , and whether a material takes the intact (ωs < 1) or thedamaged (ωs > 1) branch in that range is determined by whether or not the material has experienced astress higher than σp in the past. Namely, the viscosity is stress-history dependent.
The viscosity in our model depends on stress history in a way diﬀerent from that assumed in earlier
void-volatile models (Bercovici 1998; Foley et al. 2014). In the earlier models, ω is in the range from0 to 1, and the viscosity depends on ω as 1 − ω. The steady solution ωs is, therefore, always a singlevalued function of the stress σII ; there is no hysteresis like the one illustrated in Figure A1b. Instead, a3stress-history-dependence arises in these models by assigning a small value to λ: ω retains the memoryof its past value for a period of λ−1 (see Equation (2)), which is assumed to be on the order of 100 Myr(Bercovici 1998). These earlier works, however, suggest that it is still diﬃcult to produce the plate-likebehaviour of the lithosphere in this way. In contrast, the lithosphere retains the memory of past σ foran indeﬁnitely long period of time, as long as σm < σII < σp holds in our model (Figure A1). This morepersistent memory is the reason why the lithosphere shows a more plate-like behaviour in our model.

Inﬂuences of rheological parameters on the occurrence of plate-like motions
As discussed in our earlier paper (Ogawa 2003), rigid motion of tectonic plates separated by weak platemargins is expected to occur along the surface boundary, when the stress σcTBL in the cold thermalboundary layer (TBL) satisﬁesσm < σcTBL < σp .

(6)

For the case shown in Figures 1-3 (E = 19.3, F = 8.06), the dimensionless values of σm and σp are2.5 × 106 and 7.4 × 106 , respectively. On the other hand, assuming that ridge push force is the major
source of stress in the cold TBL, the dimensionless stress σcTBL is∫ z∗ =d∗d∗ 2∗
ρ∗0 α∗ (∆TcTBL− T ∗ )g ∗ dz ∗η0∗ κ∗ z∗ =0∫ z=1∆TcTBL ℓ= Ra,(∆TcTBL − T )dz ≈ Ra
2z=0σcTBL ≈(7)

where ℓ is the dimensionless thickness of the cold TBL, and ∆TcTBL is temperature contrast across thecold TBL. By taking ∆TcTBL ∼ 0.6, ℓ ∼ 0.1 (see Figure 1) and Ra = 1.8 × 108 (the value assumed here),
the dimensionless value of σcTBL is about 5.4 × 106 , which satisﬁes the inequality (6).
When the criterion of σcTBL < σp is violated, we obtain a convection of the “weak plate” regime (seeFigure A2). The assumed parameter values are E = 14.9, F = 8.06, and Ra = 4.6 × 108 . The resulting
σcTBL and σp are 1.86 × 107 and 5.6 × 106 , respectively. In this regime, the “tectonic plates” are fragile
and are highly damaged by the stress induced by the convective ﬂow in the mantle, and are involved into
the overall convection. The occurrence of the “weak plate” regime is consistent with the regime diagram
summarised in our earlier two-dimensional model (Ogawa 2003). In Figure A3, on the other hand, weshow convection in the “stagnant lid” regime calculated at the same values of Ra, E and F as those forthe case shown in Figures 1-3 (Ra = 1.8 × 108 , E = 19.3, F = 8.06) but started from a diﬀerent initialcondition where ω is around 0 everywhere in the top cold TBLs. The convection occurs in the stagnant
4lid regime where a stiﬀ and stagnant lid of highly viscous ﬂuid is formed along the top cold surface andthe convection occurs beneath the lid. The occurrence of both the plate-like regime and “stagnant lid”
regime at the same values of rheological parameters is consistent with our earlier ﬁnding (Ogawa 2003),
and is due to the hysteresis illustrated in Figure A1.

References
Bercovici D (1998) Generation of plate tectonics from lithosphere-mantle ﬂow and void-volatile selflubrication. Earth Planet Sci Lett 154(1-4):139–151, DOI 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, DOI
10.1029/2000JB900432
Bercovici D, Tackley PJ, Ricard Y (2015) The generation of plate tectonics from mantle dynamics. In:Schubert G (ed) Treatise on Geophysics, vol 7, 2nd edn, Elsevier, DOI 10.1016/B978-0-444-538024.00135-4
Foley BJ, Bercovici D, Elkins-Tanton LT (2014) Initiation of plate tectonics from post-magma ocean
thermochemical convection. J Geophys Res Solid Earth 119(11):8538–8561, DOI 10.1002/2014JB011121
Kameyama M (2005) ACuTEMan: A multigrid-based mantle convection simulation code and its optimization to the Earth Simulator. J Earth Simulator 4:2–10
Kameyama M, Kageyama A, Sato T (2005) Multigrid iterative algorithm using pseudo-compressibility forthree-dimensional mantle convection with strongly variable viscosity. J Comput Phys 206(1):162–181,
DOI 10.1016/j.jcp.2004.11.030
Ogawa M (2003) Plate-like regime of a numerically modeled thermal convection in a ﬂuid withtemperature-, pressure-, and stress-history-dependent viscosity. J Geophys Res 108(B2):2067, DOI
10.1029/2000JB000069
Tackley PJ (2000) Self-consistent generation of tectonic plates in time-dependent, three-dimensional
mantle convection simulations 2. strain weakening and asthenosphere. Geochem Geophys Geosyst
1(8):1026, DOI 10.1029/2000GC000043
Zhong S, Gurnis M, Moresi L (1998) Role of faults, nonlinear rheology, and viscosity structure ingenerating plates from instantaneous mantle ﬂow models. J Geophys Res 103(B7):15525–15268, DOI
10.1029/98JB00605
5(a) ωs versus σIIωs21028d”
ge=aF am“d11010100=4F=

ωs = 1−1−110tac“in−2−210Ftσm0t”
σt p−1−100
10σII /σref10(b) rω ≡ η(ω = ωs )/η(ω = 0) versus σIIr1ω0100F=0“intact”−1−110F=4−2−210−310−3−410−4Figure A1.

tσmσptF=8“damaged”−1−110σII /σref0010The schematic illustrations of stress-history-dependent rheology employed in this study.

(a) The plots of the relations between the second invariant σII of the stress tensor and the steadysolution of damage ωs calculated from Equation (4) for several values of parameter F , which determinesthe ω-dependence of viscosity F . (b) The plots against σII of the viscosity drop induced by damagerω ≡ η(ω = ωs )/η(ω = 0) for several values of F . The values of σII are normalised by σref . When F ≤ 4,
the value of ωs monotonously increases with increasing σII . Similarly, when 0 < F ≤ 4, the value of rω
monotonously decreases with σII . For F > 4, both ωs and rω are two-valued functions of stress in therange of σm ≤ σII ≤ σp with two branches (the intact branch with ωs < 1 and the damaged branchwith ωs > 1).

6

(a) Temperature T0T(b) Viscosity η and velocity V at top|V | = 10000

Figure A2.

0.75(c) Damage ω and velocity V at top−4 log10 (η/ηwp ) −0.5 |V | = 10000−2log10 ω

2The overall structure of thermal convection in the “weak plate” regime obtained for

Ra = 4.6 × 108 and E = 14.9. (a) The three-dimensional thermal structure. Shown are the distributionof temperature T at the two vertical side walls and at the bottom boundary (colour) and the isosurface
of T = 0.35. (b) The distribution of viscosity η (colour) and velocity V (arrows) in a horizontal plane
along the surface. The values of η are normalised by ηwp , the viscosity of intact material at the surface(note that ηwp < ηp , due to larger Ra and smaller E than those used in Figure 3). (c) The distributionof damage ω (colour) and velocity V (arrows) in a horizontal plane along the surface.

7

(a) Temperature T0T(b) Viscosity η and velocity V at top|V | = 1Figure A3.

0.75(c) Damage ω and velocity V at top−4 log10 (η/η p ) −0.5 |V | = 1−2log10 ω2The overall structure of thermal convection in the “stagnant lid” regime obtained for

Ra = 1.8 × 108 and E = 19.3 but for the diﬀerent initial condition without signiﬁcant rupture in thecold thermal boundary layer. (a) The three-dimensional thermal structure. Shown are the distributionof temperature T at the two vertical side walls and at the bottom boundary (colour) and the isosurfaceof T = 0.5. (b) The distribution of viscosity η (colour) and velocity V (arrows) in a horizontal plane
along the surface. The values of η are normalised by ηp , the viscosity of intact material at the surface.
(c) The distribution of damage ω (colour) and velocity V (arrows) in a horizontal plane along thesurface.

8