ARTICLE
Received 12 May 2016 | Accepted 5 Dec 2016 | Published 18 Jan 2017
DOI: 10.1038/ncomms14164
OPEN
Origin and evolution of the deep thermochemical
structure beneath Eurasia
N. Flament1,w, S. Williams1, R.D. Müller1, M. Gurnis2 & D.J. Bower3
A unique structure in the Earth’s lowermost mantle, the Perm Anomaly, was recently identified beneath Eurasia. It seismologically resembles the large low-shear velocity provinces
(LLSVPs) under Africa and the Pacific, but is much smaller. This challenges the current
understanding of the evolution of the plate–mantle system in which plumes rise from the
edges of the two LLSVPs, spatially fixed in time. New models of mantle flow over the last 230
million years reproduce the present-day structure of the lower mantle, and show a Perm-like
anomaly. The anomaly formed in isolation within a closed subduction network B22,000 km
in circumference prior to 150 million years ago before migrating B1,500 km westward at an
average rate of 1 cm year 1, indicating a greater mobility of deep mantle structures than
previously recognized. We hypothesize that the mobile Perm Anomaly could be linked to the
Emeishan volcanics, in contrast to the previously proposed Siberian Traps.
1 EarthByte Group, School of Geosciences, Madsen Building F09, University of Sydney, Sydney, New South Wales 2006, Australia. 2 Seismological Laboratory,
California Institute of Technology, Pasadena, California 91125, USA. 3 Institute of Geophysics, Department of Earth Sciences, ETH Zürich, Sonneggstrasse 5,
8092 Zürich, Switzerland. w Present address: School of Earth and Environmental Sciences, University of Wollongong, Northfields Avenue, Wollongong, New
South Wales 2522, Australia. Correspondence and requests for materials should be addressed to N.F. (email: nflament@uow.edu.au).
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
1
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
T
he long-wavelength structure of Earth’s lowermost mantle
is characterized by two large low-shear velocity provinces
(LLSVPs) under Africa and the Pacific, B15,000 km
(Fig. 1a) in diameter and 500–1,000 km high1,2. In addition, a
single, spatially small (Bo1,000 km in diameter, B500 km high)
deep mantle structure named the ‘Perm Anomaly’ was recently
identified through seismic tomography3. The discovery of the
Perm Anomaly poses fundamental questions about its dynamic
relationship with the much larger LLSVPs, its uniqueness, its age
and the formation of lower mantle structures in general. The
structure of the lower mantle is important for reconstructions of
the plate–mantle system4–8 in deep geological time because the
reconstructed locations of most large igneous provinces (LIPs)
and kimberlites over the past 320 million years (Myr) correlate
with the edges of present-day LLSVPs, leading to the concept of a
plume generation zone at LLSVP boundaries6. This concept has
been used to build a model of absolute motion of tectonic plates
over the Phanerozoic (past 540 Myr)8, under the assumptions
that LLSVPs are fixed9 and non-deforming10 through time.
However, numerical models suggest that the influence of sinking
slabs11,12 should result in LLSVP deformation and motion13–16
over hundreds of million years. Additionally, measurements of
differential splitting of SKS and SKKS seismic phases reveal
anisotropy along the boundary of the African LLSVP17 and along
the eastern boundary of the Perm Anomaly18, suggesting
deformation is occurring on the edges of these structures.
Finally, while SKS phases passing along the edge of
LLSVPs suggest gradients in shear velocity that are too
large to be explained by thermal variations alone19,20, how
chemically distinct LLSVPs are from the rest of the mantle
remains unclear21. Previous studies have shown that the largestscale structure of the lower mantle results from past subduction
history11–13,15,16,22,23, without quantifying the geographic match
between predicted and tomographic structures or the motion of
individual thermochemical structures. Although some geodynamic models23,24 produce a Perm-like anomaly linked to the
African LLSVP, its tectonic origin remains to be explained.
In this study, we report paleogeographically constrained
forward global mantle flow models predicting a discrete structure
similar in scale and geographical location to the recently
discovered Perm Anomaly, and quantify the match between
predicted and seismically inferred lower mantle structure across a
series of mantle flow and tomography models. In the flow models,
the Perm-like anomaly forms in isolation before 150 million years
ago (Myr ago), within a long-lived, B22,000 km-long, closed
subduction network consisting of the Mongol-Okhotsk subduction zone to the west, northern Tethys subduction zone to the
south, and east Asian subduction zone to the east. The models
predict that the discrete Perm-like anomaly has coherently
migrated westward at a rate of 1 cm year 1 over the last 150 Myr,
which is incompatible with the hypothesis that lower mantle
structures can be considered fixed and rigid over time. Because of
its past mobility, the Perm Anomaly may not be linked to the
Siberian Traps, but rather to the Emeishan volcanics.
Results
Predicted lowermost mantle temperature. We address the
questions raised by the discovery of the Perm Anomaly through
comparison of the lowermost mantle thermal structure predicted
by forward global mantle flow models constrained by tectonic
reconstructions25 (Methods) to tomography images. In dynamic
models, slabs subducting deep into the mantle deform a basal
layer, initially uniform, which is either thermal or thermochemical (Methods, Table 1). In the reference model (case 1,
Table 1), predicted present-day temperature B200 km above the
2
core–mantle boundary (CMB) is characterized by two large hightemperature regions under Africa and the Pacific and one smaller,
spatially distinct high-temperature region north-east of the
African Anomaly, under eastern Europe and western Russia
(Fig. 1b). For case 1, the predicted present-day CMB and surface
heat flow are respectively 10.4 TW and 40.3 TW, which is
consistent with constraints26. The spatial extent of the large hightemperature regions is in first-order agreement with the position
and shape of the African and Pacific LLSVPs in individual
tomography models (for example, S40RTS (ref. 27), Fig. 1a) and
in a vote map of tomography models3 (Fig. 1), and the predicted
smaller structure under eastern Europe and western Russia
matches the location of the Perm Anomaly in tomography
(Fig. 1). Visual comparison suggests the reference model better
fits the long-wavelength shape than similar models16. The edges
of the model LLSVPs tend to be hotter than their interior (Fig. 1),
which is consistent with a plume generation zone6.
Cluster analysis of mantle flow and tomography models. To
make a more quantitative comparison between lower mantle
temperature predicted by global geodynamic models with seismic
velocity anomalies of selected tomography images, we use cluster
analysis (Methods). Cluster analysis objectively classifies a set of
points into groups of points with similar variations in a given
property with depth. Following Lekic et al.3 we consider two
clusters and depths between 1,000 and 2,800 km (Fig. 2). For
S40RTS (ref. 27), the procedure reveals a low-velocity cluster
below B2,400 km depth, in which seismic velocity anomalies are
reduced to 0.9%, and a high-velocity cluster in which seismic
velocity anomalies are increased to þ 0.4% (Fig. 2b). For case 1f,
which was seismically filtered28 (Methods) for direct comparison
to tomography, a low-velocity and high-velocity cluster are also
distinct below B2,400 km depth, although predicted anomalies
are larger (down to 1.2% and up to þ 0.6%, Fig. 2d). The
geographic distribution of low-velocity clusters shows two large
LLSVPs and a Perm-like anomaly in both S40RTS and case 1f
(Fig. 2a,c), confirming that the extent and location of predicted
deep mantle structures is compatible with seismic images. Given
the small influence of seismic filtering, including on the extent
and location of the Perm-like anomaly (Fig. 2c,e), we do not
seismically filter other cases for which the clustering procedure
reveals a high-temperature and a low-temperature cluster below
B2,400 km depth (Fig. 2f).
Sensitivity of model success to parameters. We test the sensitivity to model parameters by considering 27 cases with varying
Rayleigh number, initial model age, initial slab depth, viscosity,
relative and absolute4,5,29 plate motions, and basal layer density
(Table 1). To assess model success, we introduce a ‘Perm Score’
PS (Table 1) that visually characterizes model clusters; the
method scores whether a predicted Perm-like anomaly is present
and separate from the African LLSVP (PS ¼ 2), present and
linked to the African LLSVP (PS ¼ 1), or absent (PS ¼ 0). A
Perm-like anomaly is present and separate from the African
LLSVP in 15 out of 27 cases. Moreover, this Perm-like anomaly is
the only isolated, small anomaly that forms in all of these 15
cases; consequently there must be a specific cause for the
generation of this unique feature. The Perm-like anomaly is
separate at present-day for initial model age 4200 Myr ago,
initial slab depth 4800 km, and when absolute plate motions are
based on hotspot tracks4,29 and paleomagnetic data5 as opposed
to mapping slab remnants from seismic tomography7 (Table 1).
These results confirm the influence of initial conditions and
subduction history on model results12, and we verified
elsewhere23 that models initiated with LLSVPs in the initial
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
a
60°
60°
30°
30°
0°
0°
–30°
–30°
–60°
–60°
S40RTS
dln VS (%)
–2.0
–1.6
–1.2
–0.8
–0.4
0.0
0.4
0.8
1.2
1.6
2.0
b
60°
60°
30°
30°
0°
0°
–30°
–30°
–60°
–60°
Case 1
0.0
0
0.1
350
0.2
0.3
700
0.4
1,050
0.5
1,400
0.6
1,750
0.7
2,100
0.8
0.9
2,450
1.0 T
2,800 T (°C)
Figure 1 | Lower mantle structure inferred from seismic tomography and predicted by a mantle flow model. (a) Seismic velocity anomalies at 2,677 km
depth for tomography model S40RTS (ref. 27). (b) Predicted present-day mantle temperature at 2,677 km depth for case 1. The solid gray contour indicates
a value of five, and the dashed gray contour a value of one in a vote map for tomography models3. Present-day coastlines are shown in black.
condition are consistent with the present-day mantle structure.
No Perm-like anomaly is predicted when the Rayleigh number Ra
is ten times larger or 100 times smaller than in the reference case
(Ra ¼ 7.8 107).
To assess model success beyond the prediction of a Perm-like
anomaly, we calculate the accuracy with which the global
geographic distribution of predicted model clusters reproduces
that of tomography clusters. This accuracy, defined as the ratio of
successfully predicted areas to total area (Methods, Fig. 3), is
calculated for distinct mantle flow and tomography models,
each of which is based on different assumptions and delivers
non-unique inferences of the true pattern of mantle structure.
The accuracy varies between 0.54 and 0.81 across 27 model cases
and seven tomography models (Methods, Fig. 4a), and is above
random (0.5) even for the least successful models. For each case,
we report the average accuracy for all seven tomography models,
which ranges between 0.56 and 0.76. Average accuracy decreases
with decreasing initial model age, is r0.61 when the initial model
age is 100 Myr ago or younger and when Ra is large (7.8 108),
and is between 0.61 and 0.71 when the basal layer is purely
thermal or o2.54% chemically denser than ambient mantle
(Fig. 4a; Table 1). The average accuracy is Z0.71 for all other
cases. The average accuracy for GyPSuM-S (ref. 30) (0.65) is
lower than for other tomography models (between 0.69 and 0.75),
which might reflect that GyPSuM-S (ref. 30) is an inversion for
geodynamic and mineral physics constraints in addition to the
seismic constraints used in other tomography models.
Origin of the Perm Anomaly. Having established that the
present-day lower mantle structure is well reproduced by mantle
flow, we investigate the dynamics leading to the formation of the
Perm-like anomaly. The model high-temperature clusters
(Fig. 2e,f) correspond to temperatures that are B10% higher than
ambient at 2,677 km depth (Fig. 1b). Following the evolution of
temperature at 2,677 km depth in the reference case (Fig. 5a,c,e)
reveals that despite its present-day proximity with the
African LLSVP, the incipient Perm-like anomaly formed
B190 Myr ago centred on 100°W/60°N (Fig. 5a), between three
long-lived subduction systems: Mongol-Okhotsk along Eurasia
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
3
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
Table 1 | Input parameters and output metrics of model cases.
Case
1
1f
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
Input
Ra
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 05
1.3E þ 07
7.8E þ 08
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
7.8E þ 07
g0 ðrÞ
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,10-100
1,1,1,100
1,1,10,100
1,1,1,100
1,0.1,1,100
1,1,1,100
1,1,1,100
1,1,1,100
1,1,1,10-100
1,1,1,10-100
1,0.1/1,1,100
1,0.1,1,100
1,1,1,10-100
1,0.1,1,10-100
1,1,1,100
dqch
4.24%
4.24%
3.39%
2.54%
1.70%
0.85%
N/A
4.24%
4.24%
4.24%
4.24%
4.24%
4.24%
1.70%
4.24%
4.24%
4.24%
4.24%
4.24%
4.24%
4.24%
4.24%
4.24%
4.24%
4.24%
1.70%
1.70%
N/A
a0 (Myr ago)
230
230
230
230
230
230
230
230
200
150
100
50
230
230
230
300*
230
230
230
230
230
230
230
230
230
230
230
230
Output
zislab (km)
1,400
1,400
1,400
1,400
1,400
1,400
1,400
1,400
1,400
1,400
1,400
1,400
425
425
425
500
800
1,100
1,200
1,750
1,200
1,200
1,200
1,400
1,400
1,200
1,200
1,200
zi90 (km)
425
425
425
425
425
425
425
425
425
425
425
425
425
425
425
425
425
425
660
660
425
660
660
425
425
425
660
425
C
N
N
N
N
N
N
N
N
N
N
N
N
N
N
N
Y
N
N
N
N
N
N
N
N
N
N
Y
N
R
C
C
C
C
C
C
C
F
F
F
F
F
F
F
C
D
C
C
B
B
C
C
C
C
E
A
B
C
PS
2
2
2
2
2
2
2
2
1
1
1
1
1
1
1
1
1
2
2
2
0
2
0
2
1
2
2
2
AccG
0.76
0.75
0.76
0.75
0.70
0.69
0.69
0.76
0.76
0.68
0.61
0.56
0.77
0.72
0.74
0.75
0.75
0.75
0.74
0.73
0.75
0.74
0.60
0.75
0.75
0.69
0.71
0.70
AccP
0.67
0.60
0.75
0.81
0.79
0.78
0.72
0.57
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
N/A
0.54
0.61
0.59
N/A
0.70
N/A
0.59
N/A
0.78
0.78
0.70
aP (Myr ago)
191
191
190
190
190
190
190
190
150
90
19
N/A
129
139
120
220
141
190
170
190
N/A
79
N/A
190
190
200
190
189
Ra is the Rayleigh number, Z0(r) is a pre-factor defined with respect to the reference viscosity Z0 for four layers: above 160 km, between 160 and 310 km depth, between 310 and 660 km depth and below
660 km depth (where 10-100 indicates that the reference viscosity linearly increases with depth across the lower mantle, and 0.1/1 indicates that the reference viscosity of the asthenosphere is 0.1
under oceanic plates and 1 under continental plates), drch is the chemical density of the basal layer, derived from the buoyancy ratio B, a0 is the age at which the model starts (* indicates that the same
boundary conditions are repeated between 300 and 230 Myr ago), zislab is the initial slab depth, zi90 is the depth from which the initial dip angle of slabs changes from 45° to 90°, G indicates whether a
phase change is considered at 660 km depth, R is the reconstruction, PS is a score indicating with value 0 if a Perm-like anomaly is not predicted, 1, if it is, and 2 if it is predicted and separate from the
African LLSVP, AccG is the global accuracy, AccP is the accuracy in the Perm region, and ap is the age from which the Perm-like anomaly exists. Parameters in bold are different from the reference case. See
Methods for more details.
(geodynamically preferable than along Central Asia31) to the
west, northern Tethys to the south, and east Asia to the east
(Fig. 5a,b). In this tectonic scenario32,33, subduction to the west of
the Perm-like anomaly ceases when the Mongol-Okhotsk Ocean
closes 150 Myr ago. The Perm-like anomaly then migrates
B1,500 km westward as pushed by descending slabs12,15,23
subducting under east Asia (Fig. 5c–f). The coherent translation
of the discrete, Perm-like anomaly allows us to estimate an
average motion rate of 1 cm year 1 over the last 150 Myr. This
ongoing westward flow is compatible with SKS–SKKS splitting
measurements revealing anisotropy with a fast east–west
direction in the lowermost mantle under eastern Europe and
western Russia18, potentially due to lattice preferred-orientation
of post-perovskite34. Moreover, the prediction is consistent with
deformation on the eastern boundary of the Perm Anomaly and
the presence of high seismic velocity structures to the east of the
Perm Anomaly that also reveal anisotropy in SKS–SKKS splitting
measurements18. Together with seismic observations18 and
previous models13–16, our results challenge the long-term fixity
and rigidity of deep-mantle thermochemical structures.
Discussion
The genesis of the Perm-like anomaly within a long-lived, closed
subduction network with a perimetre B22,000 km (Fig. 5a) could
explain why a single small LSVP18 is observed seismically, and
why only one isolated thermochemical pile forms in our models.
Although the geometry and timing of past plate boundaries is
4
increasingly uncertain back in geological times due to decreasing
amounts of preserved ocean floor35, this tectonic setting is unique
for the last 230 Myr. We find that the Perm-like anomaly is only
separate at present if slabs are inserted to depths 4800 km
(PS ¼ 2 in Table 1) in the initial condition at 230 Myr ago.
This suggests that the subduction network would have been
established at the latest between B330 and B280 Myr ago,
depending on slab sinking rates7. For the reference case, the
present-day Perm-like anomaly is chemically distinct (B550 km
in thickness based on 50% dense material), high temperature
(B850 km in thickness based on mantle 20% hotter than
ambient) and may actively contribute to mantle upwelling
(Fig. 5f). The predicted chemical anomaly is consistent with the
B500 km thickness of the Perm Anomaly inferred from seismic
images3, but the thermal anomaly may be overestimated in the
model. In contrast to models in which the Perm-like anomaly is
similar to LLSVPs (Fig. 2d), global seismic models have not
reported shear-velocity anomalies in the Perm Anomaly that are
as low as in the LLSVPs3 (Fig. 2b), although caution must be
exercised as the amplitudes of seismic anomalies are often poorly
constrained tomographically36. One possibility to explain the
apparent smaller shear-velocity reduction in the Perm Anomaly3
(Fig. 2b) is that it could be compositionally different from the
LLSVPs: decreasing the density of the basal layer results in a
better match of the Perm Anomaly, but a poorer global match
(Figs 3 and 4; Table 1).
Our reconstructions of past mantle flow link the formation of
the Perm Anomaly to the history of subduction around east and
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
a
b
60°
1,000
Low-VS
High-VS
Perm
1,200
60°
1,400
30°
0°
1,600
0°
Depth (km)
30°
1,800
2,000
2,200
−30°
2,400
−30°
2,600
−60°
−60°
2,800
–2.0 –1.5 –1.0 –0.5
S40RTS
0.0
0.5
1.0
1.5
0.5
1.0
1.5
dln VS (%)
c
d
60°
60°
1,000
Low-VS
High-VS
Perm
1,200
1,400
30°
1,600
0°
0°
Depth (km)
30°
1,800
2,000
2,200
−30°
−30°
2,400
2,600
−60°
−60°
2,800
–2.0 –1.5 –1.0 –0.5
Case 1f
0.0
dln VS (%)
e
f
60°
60°
1,000
High-T
Low-T
Perm
1,200
1,400
30°
1,600
0°
0°
Depth (km)
30°
1,800
2,000
2,200
−30°
−30°
2,400
2,600
−60°
−60°
Case 1
2,800
0
0.25
0.5
0.75
1.0
T (non-dimensional)
Figure 2 | Clustering of lower mantle seismic tomography and predicted mantle temperature. (a) High-velocity (blue) and low-velocity (red) regions
between 1,000 and 2,800 km depth for seismic tomography model S40RTS (ref. 27). (b) Seismic velocity profiles in high-velocity and low-velocity regions
for S40RTS (ref. 27). The solid curves are the mean, and the transparent envelopes are the associated standard deviation, of the global low-velocity cluster
(red), global high-velocity cluster (blue) and low-velocity cluster for the separate Perm-like anomaly (black). (c,e). Same as a but for seismically filtered28
case 1f (c) and case 1 (e). (d,f). Same as b but for seismically filtered28 case 1f (d) and case 1 (f). In a,c,e, the solid gray contour indicates a value of five and
the dashed gray contour a value of one in a vote map for tomography models3. Present-day coastlines are shown in black.
central Asia between 230 and 150 Myr ago. The formation of the
Perm Anomaly would have occurred several tens of million years
after this subduction network was established, as slabs slowly
sank. The Perm Anomaly appears between 220 and 80 Myr ago
depending on the Rayleigh number, the initial slab depth and the
initial model age (Methods, Table 1). Although tectonic uncertainties increase back in geological time, some reconstructions
suggest that a closed network of subduction zones B20,000 km in
perimeter might have been established around the MongolOkhotsk Ocean 410 Myr ago37, in which case the Perm Anomaly
might have existed for much of the Phanerozoic. Structures
similar to the Perm Anomaly are likely to have existed earlier in
Earth’s history, controlled by past subduction zone configurations.
Conceptual6 and geodynamic22,23,24 models suggest that
plumes mostly rise from deep thermochemical structures to
form LIPs at Earth’s surface. The reconstructed location of
theB258-Myr ago Emeishan LIP falls within the network formed
by the Mongol-Okhotsk, northern Tethys, and East Asia
subduction zones between 230 and 150 Myr ago33 (Fig. 5a). In
contrast, the reconstructed location of the B251 Myr ago
Siberian Traps does not reconstruct within this subduction
network between 230 and 150 Myr ago (the reconstructed
location of the Siberian Traps is outside the region shown in
Fig. 5a). Because the models show that the Perm anomaly
originated within this subduction network, we propose
that the Emeishan LIP is a possible product of the Perm
Anomaly, in contrast to the Siberian Traps3,8. These competing
hypotheses could be tested in future convection models
including mantle plumes23, contrary to the models presented
here, and based on tectonic reconstructions extending into the
Paleozoic, but that do not assume that the Emeishan LIP
originated from the Pacific LLSVP, contrary to existing
reconstructions8,37.
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
5
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
a
40°
b
80°
75°
60°
60°
0°
45°
–60°
Case 1/S40RTS
d
c
75°
60°
60°
0°
45°
–60°
Case 2/S362ANI
e
75°
f
60°
60°
0°
45°
–60°
Case 5/S362ANI
20°
40°
60°
80°
Figure 3 | Spatial match between predicted lower mantle structure and that inferred from seismic tomography. Orange (true positive) indicates hightemperature cluster for the model and low-velocity cluster for the tomography, gray (true negative) indicates low-temperature cluster for the model and
high-velocity cluster for the tomography, green (false positive) indicates high-temperature cluster for the model and high-velocity cluster for the
tomography and blue (false negative) indicates low-temperature cluster for the model and low-velocity cluster for the tomography. Present-day coastlines
are shown in black. Results are shown for case 1 and S40RTS (ref. 27) (a,b), case 2 and S362ANI (ref. 38) (c,d), e,f, and case 5 and S362ANI (ref. 38) (e,f).
b,d,f show results in the region between 10°–80°E and 40°–75°N that includes the Perm Anomaly.
The tectonic configuration of a B22,000 km network of longlived (480 Myr) subduction zones around east and central Asia
before 150 Myr ago, unique in the last 200 Myr ago, led to the
formation of a single, well-defined and isolated thermo-chemical
anomaly. Numerical models reproduce this past natural
experiment, and the coherent westward motion of the
discrete Perm-like anomaly allows us to quantify an average
motion of 1 cm year 1 since the Mongol-Okhotsk Ocean closed
150 Myr ago.
Methods
Paleogeographically constrained dynamic Earth models. We solve the
equations for incompressible convection in a spherical domain with finite-elements
using the code CitcomS39, modified as described in ref. 25 to progressively
assimilate the velocity of tectonic plates, the age of the ocean floor and the location
and polarity of subduction zones determined in one million year intervals from
global plate tectonic reconstructions32,33 with continuously closing plates40. This
semi-empirical approach ensures our computations represent Earth’s imposed
tectonic history, allowing us to reconstruct the history of deep mantle flow over the
last 230 Myr. This approach is guided by the current intractability of computing
time-dependent models of Earth’s plate–mantle system with the resolution
required to dynamically achieve tectonic-like features, including one-sided
subduction41 and conserve the energy of the system simultaneously. Here we give a
summary of the governing parameters and model setup that are further described
in refs 25,42.
The vigour of convection
is defined by the Rayleigh number
Ra ¼ a0 r0 g0 DTh3M ðk0 Z0 Þ, where a0 ¼ 310 5 K 1 is the coefficient of
6
thermal expansion, r0 ¼ 4; 000kg m 3 is the density, g0 ¼ 9:81m s 2 is the
gravity acceleration, DT ¼ 2; 825K is the temperature change across the mantle,
hM ¼ 2; 867km is the thickness of the mantle, k0 ¼ 110 6 m2 s 1 is the thermal
diffusivity, Z0 ¼ 11021 Pa s is the viscosity, and the subscript ‘0’ indicates
reference values. With the values listed above, Ra ¼ 7:8107 . These values are
varied between model cases such that Ra varies between 7:8105 and 7:8109
(Table 1).
We approximate the Earth’s mantle as a Newtonian fluid in
which viscosity
varies
with depth
and
temperature
following
Z ¼ Z0 ðr Þ exp EZ R T þ TZ EZ R Tb þ TZ , where Z0 ðr Þ is a pre-factor
defined with respect to the reference viscosity Z0 for four layers: above 160 km,
between 160 and 310 km depth, between 310 and 660 km depth and below 660 km
depth, in the lower mantle. Values of Z0 ðr Þ for each layer are given as commaseparated lists in Table 1, where ‘10-100’ indicates that the reference viscosity
linearly increases with depth from 10 to 100 throughout the lower mantle, and
0.1/1 indicates that the reference viscosity of the asthenosphere is 0.1 under oceanic
plates and 1 under continental plates. EZ is the activation energy taken as
100 kJ mol 1 in the upper mantle and 30 kJ mol 1 in the lower mantle,
R ¼ 8.31 J mol 1 K 1 is the universal gas constant, T is the dimensional
temperature, TZ ¼ 452K is a temperature offset and Tb ¼ 1; 685K is the ambient
mantle temperature. The activation energy and temperature offset are chosen to
limit variations in viscosity to three orders of magnitude across the range of
temperatures without imposing a yield stress. Such lateral viscosity contrasts are
lower than expected to occur within the solid Earth41, but they can be computed
with a resolution that allows us to compute time-dependent mantle flow models. A
phase change G at 660 km depth, as described in Flament et al.42 is considered in
some model cases (Table 1).
For the initial condition and progressive data assimilation, the thickness and
temperature of the lithosphere are derived using a half-space cooling model and the
synthetic age of the ocean floor25, and simplified tectonothermal ages for the
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
a
Global
0.81
1
1f
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
SEMUCB-WM1 SAVANI
0.78
0.75
0.69
0.66
Accuracy
0.72
0.63
0.60
0.57
i
(either 425 or 660 km; Table 1). Slabs are initially twice as thick in the lower
z90
mantle compared with their thickness in the upper mantle, to account for advective
thickening in the more viscous lower mantle.
The model consists of 129 129 65 12 E13 106 nodes, which with a
radial mesh refinement that gives average resolutions ofB50 50 15 km at the
surface,B28 28 27 km at the CMB, andB40 40 100 km in the mid-mantle.
We investigate the influence of relative and absolute plate motions across five
global tectonic reconstructions (R in Table 1). Reconstruction A, which uses the
absolute plate motions of ONeill et al.4 between 0 and 100 Myr ago and that of
Steinberger and Torsvik5 before 100 Myr ago, is described in Seton et al.32 and
extended from the last 200 Myr ago to the last 230 Myr ago33. Continuouslyclosing plate polygons40 and ages of the ocean floor32,33, necessary to assimilate
plate reconstructions in mantle flow models with the method described in
Bower et al.25 are available to us back to 230 Myr ago. Reconstruction B includes
changes to relative plate motions in the Arctic region33 and uses the absolute
plate motions of Torsvik et al.29 between 0 and 70 Myr ago and that of
Steinberger and Torsvik5 before 105 Myr ago, with interpolation between the
two absolute plate motion models between 70 and 105 Myr ago. Reconstruction
C includes changes to relative plate motions in Southeast Asia33 compared with
reconstruction B. Reconstruction D includes changes to relative plate motions in
the western Tethys33 compared with reconstruction C. Reconstruction E is based
on the absolute plate motions of van der Meer et al.7 and on the same relative
plate motions as reconstruction D. Reconstruction F uses the same absolute plate
motion model as reconstruction B, and relative plate motions as described in
Muller et al.33
0.54
S40RTS
b
GYPSUMS S362ANI
HMSL-S SAW24B16
Perm
1
0.84
1f
2
0.78
3
4
6
0.66
7
17
0.60
18
19
Accuracy
0.72
5
0.54
21
23
0.48
25
0.42
26
Cluster analysis of lower mantle structure. We use cluster analysis to objectively
classify a set of points on the surface into groups of points with similar variations in
a given property with depth. For each flow model case (or tomography model),
temperature (or seismic velocity) profiles are treated as 196,596 independent
vectors of 31 coordinates specifying the temperatures (or seismic velocities)
sampled at 31 depths between 1,000 and 2,800 km3, with an average resolution
of 58 km. Each vector corresponds to equally-spaced locations on Earth’s surface
(average distance B0.45°). The vectors are grouped into two clusters using
k-means clustering43, a procedure that keeps the variance in squared Euclidean
distance between vectors small within each cluster. We use the scientific Python
implementation of the k-means algorithm (http://docs.scipy.org/doc/scipy/
reference/generated/scipy.cluster.vq.kmeans2.html). The average and standard
deviation of temperature profiles for each cluster are shown in Fig. 2 for the
reference case, along with the average and standard deviation of temperature
profiles for the high-temperature cluster in the Perm region.
27
SEMUCB-WM1 SAVANI
S40RTS
GYPSUMS S362ANI
HMSL-S SAW24B16
Figure 4 | Quantitative match between predicted lower mantle structure
and seismic tomography models. Accuracy (true positive area plus true
negative area over total area), with which the geographical distribution of
clusters of mantle temperature between 1,000 and 2,800 km depth
predicted by 27 mantle flow model cases (case 1f is seismically-filtered case
1) reproduce the geographical distribution of clusters of seismic velocity
anomalies between 1,000 and 2,800 km depth for seven S-wave
tomography models. (a) global accuracy, (b) accuracy in the region
between 10°–80°E and 40°–75°N that includes the Perm Anomaly, for
cases predicting a Perm-like anomaly separate from the model African
LLSVP (PS ¼ 2 in Table 1). See Table 1 for values of the global and regional
accuracies averaged over the seven considered tomography models.
continental lithosphere42. The global thermal structure of slabs is constructed from
the location of subduction zones and from the age of the ocean floor25. The global
thermal structure of the lithosphere and of subducting slabs is assimilated in the
dynamic models in 1 Myr increments, to 350 km depth at subduction zones25. In
the initial condition, subduction zones are inserted in the mantle assuming a
descent rate of 3 cm year 1 in the upper mantle and 1.2 cm year 1 in the lower
mantle7. Subduction zones that appear during the model are progressively inserted
to 350 km depth based on the age of subduction initiation and on the plate
convergence rate.
The initial condition, derived from the tectonic reconstruction at 230 Myr ago,
includes a basal layer just above CMB, which is either purely thermal or
thermochemical (Table 1). The layer is 113 km thick, which represents 2% of the
volume of the mantle, consistent with the seismically inferred value2. The
composition of that layer is modeled using tracers42 and its chemical density drch is
varied between þ 0.85% and þ 4.24% (Table 1) by changing the buoyancy ratio
B ¼ drch =ðr0 a0 DT Þ between 0.1 and 0.5, with increment 0.1. Slabs are initially
i
inserted down to a depth of either zslab
(varied between 425 km and 1,750 km;
Table 1), or to the depth derived from their initiation age and sinking rates if that
i
i
depth is shallower than zslab
, with a dip of 45° down to z90
and a dip of 90° below
Seismic filtering. We seismically filter our reference case 1 following
Ritsema et al.28 to verify if lower mantle features apparent in the model temperature
field would be resolved by global tomography27. We consider that both temperature
and composition variations cause variations in shear velocity. For the thermal
contribution, we determine wave speed variations (dVS) scaling departures from
average temperature at each depth (dT) using dVS/dT ¼ 7.0 10 5 km s 1 K 1
(ref. 44). For the chemical contribution, we determine wave speed variations (dVS)
from variations in the composition field at each depth (dC) using dVS/dC ¼ 4 drch
(ref. 45), where drch is the density of the chemically distinct basal layer. Thermal and
compositional contributions to wave speed variations are jointly considered following
dVs ¼ dVs/dT dT þ fc dVs/dC dC, where fc ¼ 0.05 is the contribution of
composition to total wave speed anomalies.
Quantification of model success. We quantify how well the clusters obtained for
each model case reproduce the global geographic distribution of clusters obtained
for global S-wave tomography models27,30,38,46–49 by computing the accuracy
Acc ¼ ðTP þ TNÞ=A, where TP is the area of true positives, TN the area of true
negatives and A the total area. Three examples of the spatial distribution of true
positives, true negatives, false positives and false negatives with respect to
tomography models are shown in Fig. 3. The accuracy is computed for each model
case, both globally and in the region between 10°–80°E and 40°–75°N that includes
the Perm Anomaly, against seven global S-wave tomography models (Fig. 4):
SAW24B16 (ref. 46), HMSL-S (ref. 47), S362ANI (ref. 38), GyPSuM-S (ref. 30),
S40RTS (ref. 27), Savani (ref. 48), SEMUCB-WM1 (ref. 49). Values of the global
and regional accuracies averaged over the seven tomography models are reported
in Table 1.
Age of the Perm-like anomaly. We report the age ap from which the separate
Perm-like anomaly is 4190 km thick (based on a mantle temperature 10% higher
than ambient, and with model output everyB10 Myr) in Table 1.
Data availability. Maps of the geographic distribution of tomography and
flow model clusters reported are available at https://www.earthbyte.org/originevolution-perm-anomaly/. The computer code that supports the findings of this
study is available from the corresponding author upon reasonable request.
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
7
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
a
b
60°
191 Ma
°
120° 1
60
°
00
2,6
1,3
80
0
30°
00
c
d
60°
111 Ma
°
e
00
00
2,6
1,3
40
80° 120
°
0
30°
f
60°
0 Ma
180°
120°
0.0
0
0.1
350
0.2
0.3
700
0.4
1,050
0.5
1,400
0.6
1,750
0.7
2,100
0.8
0.9
2,450
00
2 cm yr–1
2 cm yr–1
00
60°
2,6
0°
40° 8
0°
1,3
0°
0
30°
1.0 T
2,800 T (°C)
Figure 5 | Formation of the Perm Anomaly. Predicted mantle temperature, flow and composition for case 1 at 2,677 km depth (a,c,e) and along crosssections (green lines) at 60° latitude (b,d,f). Results are shown at 191 Myr ago (a,b), 111 Myr ago (c,d) and present-day (e,f). Reconstructed subduction
locations are shown as red lines with triangles on the overriding plate, reconstructed mid-oceanic ridges and transform faults are shown as yellow lines,
and reconstructed coastlines are shown in gray. The brown contours indicate 50% concentration of dense material. In a,c,e, the blue contours indicate
subducting plates with temperature 10% lower than ambient mantle at 489 km depth. The green circle in a is the location of the B258 Myr ago
Emeishan LIP, reconstructed at 191 Myr ago.
References
1. Garnero, E. J. & McNamara, A. K. Structure and dynamics of Earth’s lower
mantle. Science 320, 626–628 (2008).
2. Hernlund, J. W. & Houser, C. On the statistical distribution of seismic velocities
in Earth’s deep mantle. Earth Planet. Sci. Lett. 265, 423–437 (2008).
3. Lekic, V., Cottaar, S., Dziewonski, A. & Romanowicz, B. A. Cluster analysis of
global lower mantle tomography: a new class of structure and implications for
chemical heterogeneity. Earth Planet. Sci. Lett. 357, 68–77 (2012).
4. O’Neill, C., Müller, R. D. & Steinberger, B. On the uncertainties in hot spot
reconstructions and the significance of moving hot spot reference frames.
Geochem. Geophys. Geosyst. 6, Q04003 (2005).
5. Steinberger, B. & Torsvik, T. H. Absolute plate motions and true polar wander
in the absence of hotspot tracks. Nature 452, 620–623 (2008).
6. Torsvik, T. H., Burke, K., Steinberger, B., Webb, S. J. & Ashwal, L. D. Diamonds
sampled by plumes from the core-mantle boundary. Nature 466, 352–355
(2010).
7. van der Meer, D. G., Spakman, W., van Hinsbergen, D. J. J., Amaru, M. L. &
Torsvik, T. H. Towards absolute plate motions constrained by lower-mantle
slab remnants. Nat. Geosci. 3, 36–40 (2010).
8. Torsvik, T. H. et al. Deep mantle structure as a reference frame for movements
in and on the Earth. Proc. Natl Acad. Sci USA 111, 8735–8740 (2014).
9. Dziewonski, A. M., Lekic, V. & Romanowicz, B. A. Mantle anchor structure: an
argument for bottom up tectonics. Earth Planet. Sci. Lett. 299, 69–79 (2010).
10. Burke, K. & Torsvik, T. H. Derivation of large igneous provinces of the past 200
million years from long-term heterogeneities in the deep mantle. Earth Planet.
Sci. Lett. 227, 531–538 (2004).
11. Bunge, H.-P. et al. Time scales and heterogeneous structure in geodynamic
Earth models. Science 280, 91–95 (1998).
12. McNamara, A. K. & Zhong, S. Thermochemical structures beneath Africa and
the Pacific Ocean. Nature 437, 1136–1139 (2005).
13. Zhang, N., Zhong, S., Leng, W. & Li, Z. X. A model for the evolution of the
Earth’s mantle structure since the Early Paleozoic. J. Geophys. Res. Solid Earth
115, B06401 (2010).
14. Tan, E., Leng, W., Zhong, S. & Gurnis, M. On the location of plumes and lateral
movement of thermochemical structures with high bulk modulus in the 3-D
compressible mantle. Geochem. Geophys. Geosyst. 12, Q07005 (2011).
8
15. Bower, D. J., Gurnis, M. & Seton, M. Lower mantle structure from
paleogeographically constrained dynamic Earth models. Geochem. Geophys.
Geosyst. 14, 44–63 (2013).
16. Zhong, S. & Rudolph, M. L. On the temporal evolution of long-wavelength
mantle structure of the Earth since the early Paleozoic. Geochem. Geophys.
Geosyst. 16, 1599–1615 (2015).
17. Lynner, C. & Long, M. D. Lowermost mantle anisotropy and deformation along
the boundary of the African LLSVP. Geophys. Res. Lett. 41, 3447–3454 (2014).
18. Long, M. D. & Lynner, C. Seismic anisotropy in the lowermost mantle near the
Perm Anomaly. Geophys. Res. Lett. 42, 7073–7080 (2015).
19. Ni, S., Tan, E., Gurnis, M. & Helmberger, D. V. Sharp sides to the African
superplume. Science 296, 1850–1852 (2002).
20. Wang, Y. & Wen, L. Geometry and P and S velocity structure of the ‘African
Anomaly’. J. Geophys. Res. Solid Earth 112, B05313 (2007).
21. Davies, D. R., Goes, S. & Lau, H. C. P. in The Earth’s Heterogeneous Mantle
(eds Khan, A. & Deschamps, F.) 441–477 (Springer, 2015).
22. Davies, D. R. et al. Reconciling dynamic and seismic models of Earth’s lower
mantle: the dominant role of thermal heterogeneity. Earth Planet. Sci. Lett. 353,
253–269 (2012).
23. Hassan, R., Flament, N., Gurnis, M., Bower, D. J. & Müller, R. D. Provenance
of plumes in global convection models. Geochem. Geophys. Geosyst. 16,
1465–1489 (2015).
24. Steinberger, B. & Torsvik, T. H. A geodynamic model of plumes from the
margins of Large Low Shear Velocity Provinces. Geochem. Geophys. Geosyst. 13,
Q01W09 (2012).
25. Bower, D. J., Gurnis, M. & Flament, N. Assimilating lithosphere and slab
history in 4-D Earth models. Phys. Earth Planet. Inter. 238, 8–22 (2015).
26. Jaupart, C., Labrosse, S. & Mareschal, J. in Treatise on Geophysics: Mantle
Dynamics (ed. Bercovici, D.) 253–303 (Elsevier, 2007).
27. Ritsema, J., Deuss, A., van Heijst, H. J. & Woodhouse, J. H. S40RTS: a degree-40
shear-velocity model for the mantle from new Rayleigh wave dispersion,
teleseismic traveltime and normal-mode splitting function measurements.
Geophys. J. Int. 184, 1223–1236 (2011).
28. Ritsema, J., McNamara, A. K. & Bull, A. L. Tomographic filtering of
geodynamic models: implications for model interpretation and large-scale
mantle structure. J. Geophys. Res. Solid Earth 112, B01303 (2007).
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
ARTICLE
NATURE COMMUNICATIONS | DOI: 10.1038/ncomms14164
29. Torsvik, T. H., Müller, R. D., van der Voo, R., Steinberger, B. & Gaina, C.
Global plate motion frames: toward a unified model. Rev. Geophys. 46, RG3004
(2008).
30. Simmons, N. A., Forte, A. M., Boschi, L. & Grand, S. P. GyPSuM: a joint
tomographic model of mantle density and seismic wave speeds. J. Geophys. Res.
Solid Earth 115, B12310 (2010).
31. Fritzell, E. H., Bull, A. L. & Shephard, G. E. Closure of the Mongol–Okhotsk
Ocean: insights from seismic tomography and numerical modelling. Earth
Planet. Sci. Lett. 445, 1–12 (2016).
32. Seton, M. et al. Global continental and ocean basin reconstructions since
200Ma. Earth-Sci. Rev. 113, 212–270 (2012).
33. Müller, R. D. et al. Ocean basin evolution and global-scale plate reorganization
events since pangea breakup. Annu. Rev. Earth Planet. Sci. 44, 107–138 (2016).
34. Ford, H. A., Long, M. D., He, X. & Lynner, C. Lowermost mantle flow at the
eastern edge of the African Large Low Shear Velocity Province. Earth Planet.
Sci. Lett. 420, 12–22 (2015).
35. Torsvik, T. H., Steinberger, B., Gurnis, M. & Gaina, C. Plate tectonics and net
lithosphere rotation over the past 150My. Earth Planet. Sci. Lett. 291, 106–112
(2010).
36. Helmberger, D. V. & Ni, S. in Earth’s Deep Mantle: Structure, Composition, and
Evolution Vol. 160 Geophysical Monograph Series 63–81 (American
Geophysical Union, 2005).
37. Domeier, M. & Torsvik, T. H. Plate tectonics in the late Paleozoic. Geosci.
Front. 5, 303–350 (2014).
38. Kustowski, B., Ekström, G. & Dziewoński, A. M. Anisotropic shear-wave
velocity structure of the Earth’s mantle: a global model. J. Geophys. Res. Solid
Earth 113, B06306 (2008).
39. Zhong, S., McNamara, A., Tan, E., Moresi, L. & Gurnis, M. A benchmark study
on mantle convection in a 3-D spherical shell using CitcomS. Geochem.
Geophys. Geosyst. 9, Q10017 (2008).
40. Gurnis, M. et al. Plate tectonic reconstructions with continuously closing plates.
Comput. Geosci. 38, 35–42 (2012).
41. Stadler, G. et al. The dynamics of plate tectonics and mantle flow: from local to
global scales. Science 329, 1033–1038 (2010).
42. Flament, N. et al. Topographic asymmetry of the South Atlantic from global
models of mantle flow and lithospheric stretching. Earth Planet. Sci. Lett. 387,
107–119 (2014).
43. MacQueen, J. in Proceedings of the Fifth Berkeley Symposium on Mathematical
Statistics and Probability 281–297 (Oakland, CA, USA, 1967).
44. Forte, A. M. & Mitrovica, J. X. Deep-mantle high-viscosity flow and
thermochemical structure inferred from seismic and geodynamic data. Nature
410, 1049–1056 (2001).
45. Steinberger, B. Effects of latent heat release at phase boundaries on flow in the
Earth’s mantle, phase boundary topography and dynamic topography at the
Earth’s surface. Phys. Earth Planet. Inter. 164, 2–20 (2007).
46. Mégnin, C. & Romanowicz, B. A. The three-dimensional shear velocity
structure of the mantle from the inversion of body, surface and higher-mode
waveforms. Geophys. J. Int. 143, 709–728 (2000).
47. Houser, C., Masters, G., Shearer, P. & Laske, G. Shear and compressional
velocity models of the mantle from cluster analysis of long-period waveforms.
Geophys. J. Int. 174, 195–212 (2008).
48. Auer, L., Boschi, L., Becker, T. W., Nissen-Meyer, T. & Giardini, D. Savani:
a variable resolution whole-mantle model of anisotropic shear velocity
variations based on multiple data sets. J. Geophys. Res. Solid Earth 119,
3006–3034 (2014).
49. French, S. W. & Romanowicz, B. A. Whole-mantle radially anisotropic shear
velocity structure from spectral-element waveform tomography. Geophys. J. Int.
199, 1303–1327 (2014).
Acknowledgements
This research was undertaken with the assistance of resources from the National
Computational Infrastructure (NCI), which is supported by the Australian Government.
N.F. and R.D.M. were supported by ARC IH130200012. S.W. was supported by SIEF RP
04-174. M.G. and D.J.B. were supported by Statoil ASA and by the NSF under grants
CMMI-1028978, EAR-1161046 and EAR-1247022. Figures were constructed using the
Generic Mapping Tools and matplotlib. We thank T.C.W. Landgrebe for advice in
calculating accuracy and sensitivity, D. Steinberg for advice on cluster analysis, J. Ritsema
for sharing the seismic filter for S40RTS, and three anonymous reviewers for comments
that improved the quality of the manuscript.
Author contributions
N.F. ran and analysed the models, S.W. instigated the study and implemented the cluster
analysis, R.D.M. and M.G. developed the concepts of the study and D.J.B., M.G. and N.F.
developed the framework to assimilate tectonic reconstructions in CitcomS. All authors
contributed to the writing of the manuscript, led by N.F.
Additional information
Competing financial interests: The authors declare no competing financial
interests.
Reprints and permission information is available online at http://npg.nature.com/
reprintsandpermissions/
How to cite this article: Flament, N. et al. Origin and evolution of the deep thermochemical structure beneath Eurasia. Nat. Commun. 8, 14164 doi: 10.1038/ncomms14164
(2017).
Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in
published maps and institutional affiliations.
This work is licensed under a Creative Commons Attribution 4.0
International License. The images or other third party material in this
article are included in the article’s Creative Commons license, unless indicated otherwise
in the credit line; if the material is not included under the Creative Commons license,
users will need to obtain permission from the license holder to reproduce the material.
To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/
r The Author(s) 2017
NATURE COMMUNICATIONS | 8:14164 | DOI: 10.1038/ncomms14164 | www.nature.com/naturecommunications
9