# Geodynamics of a global plate reorganization from topological data analysis

Alexandre Janin, Nicolas Coltice, Nicolas Chamot-Rooke, Julien Tierny

► To cite this version:

Alexandre Janin, Nicolas Coltice, Nicolas Chamot-Rooke, Julien Tierny. Geodynamics of a global plate reorganization from topological data analysis. Nature Geoscience, 2025, 18 (10), pp.1041-1047. 10.1038/s41561-025-01772-7. hal-05338073

HAL Id: hal-05338073

https://hal.science/hal-05338073v1

Submitted on 30 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Geodynamics of a global plate reorganization from topological data analysis

Received: 16 August 2024

Accepted: 18 July 2025

Published online: 20 August 2025

<figure><img src="image_3.png" /></figure>

Check for updates

Alexandre Janin<fnref n="1" />,²✉, Nicolas Coltice<fnref n="1" />,³, Nicolas Chamot-Rooke<fnref n="1" /> & Julien Tierny⁴

Earth's lithosphere is fragmented into tectonic plates that move relative to each other and the mantle, following a dynamic equilibrium that evolves over time. Sporadically, abrupt kinematic changes occur that typically lead to regional-scale plate reconfigurations, and hypothetically to global-scale events, affecting the motion of all major plates. Identifying and understanding such global plate reorganizations has proven challenging because of the geometric nature of plate tectonic theory and the scarcity of preserved kinematic changes. Here we combine Earth-like global mantle convection with topological data analysis to show that global plate reorganization events can emerge from mantle-lithosphere interactions. Unlike regional changes, during these rare global events, the entire plate network shifts to a new dynamic equilibrium by the formation and destruction of plate boundaries at global scale, affecting the dynamics of the entire mantle. Such events represent transitions between two stable states in the coupled mantle-plates system. Our findings suggest that Earth may have experienced a similar global tectonic reorganization during the mid-Eocene around 47 million years ago. Our new approach provides a quantitative framework to unveil mantle-lithosphere interactions from global geodynamic models.

Plate tectonics successfully describes the instantaneous motion of lithospheric plates, as demonstrated by global kinematic models<fnref n="1" />⁻³. However, plate motion is not constant, and sporadically, tectonic plates experience abrupt changes in the norm and/or direction of their motion⁴,⁵. These abrupt kinematic changes, or plate reorganization events, are an important feature of Earth's plate motion⁶,⁷, well recorded in the oceanic lithosphere⁸⁻¹¹. Although generally limited to one plate and its immediate surrounding, the observation of strong temporal concordances of regional kinematic changes involving together a large fraction of the plate network, one in the mid-Cretaceous¹²⁻¹⁴ and the other in the mid-Eocene¹⁵⁻¹⁹, has led to the hypothesis that plate reorganizations occasionally occur on a global scale, affecting the motion of major lithospheric plates.

How can global reorganizations be identified unequivocally? What is the chronology of events through a reorganization? How long

do they take? What forces drive them? Answering such key questions has so far been limited by our approach and the observations available. Indeed, the study of plate reorganizations is impeded by the non-dynamical but geometrical nature of plate tectonic theory, the rarity of kinematic changes preserved in the geological record and the difficulty of identifying them. To overcome these hurdles, we combine a new methodology coupling topological data analysis²⁰⁻²² and geodynamics with a fully dynamical state-of-the-art mantle convection model computed by ref. 23 that self-consistently generates Earth major tectonic features. This approach allows us to detect automatically and time-track plates emerging at the model surface. It offers an alternative to plate reconstruction models in the study of plate kinematic changes, avoiding the limitations linked to the sparsity of observations. It further allows us to explore the dynamic links between the lithosphere and the mantle during reorganizations.

<footnote n="1">¹Laboratoire de Géologie-CNRS UMR 8538, École normale supérieure-PSL University, Paris, France. ²Department of Earth and Environmental Sciences, Boston College, Chestnut Hill, MA, USA. ³Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Géoazur, Valbonne, France. ⁴CNRS, Sorbonne Université, Laboratoire d'Informatique de Paris 6, LIP6, Paris, France. ✉e-mail: alexandre.janin@bc.edu</footnote>
a

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

**Fig. 1 | Regional and global kinematic changes.** Surface maps of the model and synthetic hotspot tracks. a, Model surface at 290 Myr (38 Myr after discussed 252-Myr RPRE). b, Model surface at 190 Myr (40 Myr after the end of the GPRE). The synthetic hotspot tracks (labelled from A to I) represent the tracks of physical particles advected with the lithosphere and emitted every 2 Myr at the same location (black open circles), fixed relative to the mantle. The colour bars indicate the time at which the particles were generated. Stars highlight the particles produced at *t*<sub>ref</sub>; *t*<sub>ref</sub> = 252 Myr, is the time of the discussed RPRE

b

<figure><img src="image_6.png" /></figure>

(a); *t*<sub>ref</sub> = 136 Myr, is the time of the main event within the GPRE (b) (text). The time series surrounding the maps show the temporal evolution of the azimuth (in blue) and the L2 norm (in red) of the velocity (in NNR: no net rotation of the lithosphere) of the oldest particle for each track. The vertical orange dashed lines indicate *t*<sub>ref</sub>. The temporal extent of the GPRE is shown as a light orange area on the times series of b. Green arrows indicate other local kinematic changes. Divergence and absolute vorticity fields are displayed below the detected boundaries. Cont., continental.

The model analysed here is the only 3D spherical mantle convection calculation published yet that self-consistently generates the major tectonic features of Earth at Earth-like convective vigour²³. It represents a unique high-resolution numerical solution providing 1,500 Myr of coupled evolution between tectonics and mantle dynamics with an average surface horizontal resolution of 23 km and less than 4,000 years per time step. It successfully reproduces major features and complexity of the present-day Earth tectonics such as supercontinent cycles (including extro and introversion), one-sided subductions, slab roll-back and advance, slab stagnation around the 660-km boundary, back-arc extension, concurrent rifting and collision of continents, interactions of plumes with mid-ocean ridges and slabs, diffuse and sharp plate boundaries²³ (Fig. 1, Extended Data Figs. 1 and 2 and Methods). Performing a plate tectonic analysis automatically on such a model is a topological and computing challenge. We develop here a new methodology (Methods) to automatically generate a plate tessellation every 2 Myr using the surface velocity field and a topological and kinematic analysis. Each plate is then a rigid polygon surrounded by a high-velocity gradient and is tracked over time assuming independent time steps. This new methodology enables us to document the kinematic history of all detected tectonic plates over 474 Myr and follow their fragmentation and fusion using a global evolution graph comprising 30,000 vertices (approximately 126 plates detected per time step). By analysing surface tectonics and mantle dynamics at this scale, we can address

the key questions above through an integrated view of the coupled
mantle–lithosphere system.

**Kinematic change on a global scale**

Plates detected at the surface of the model show first-order similarities with Earth tectonics in size, shape and motion (Fig. 1 and Extended Data Fig. 2). Figure 1 tracks the motion of larger plates during two selected time windows using synthetic hotspot tracks, that is, tracks of Lagrangian particles generated every 2 Myr at the same location in the mantle and advected with plate motion. Similar to Earth, the motion of model's plates is quasi-stationary over periods that can last several tens of millions of years (Fig. 1), consistent with our understanding of long-lived plate forces²⁴,²⁵, separated by abrupt transitions that purely convective mechanisms cannot explain alone and during which the norm and/or the direction of the plate motion changes⁸⁻¹¹. A typical example of an abrupt event is the sudden fragmentation at time 252 Myr of the plate bearing synthetic hotspot tracks B and C responsible for a local kinematic change that affects both the azimuth and the norm of the regional velocity field (Fig. 1a and Extended Data Fig. 3). The reorganization is brief, less than 2 Myr, and the motion of the other plates remains unaffected. This type of regional plate reorganization event (RPRE) is analogous to documented Earth's kinematic changes⁸⁻¹¹ and appears to occur commonly at the surface of the model and the surface of Earth (Fig. 1).

In contrast to these regional events, one major global event occurs
at 136 Myr (Fig. 1b) involving changes in the motion of all the major
<figure><img src="image_7.png" /></figure>

**Fig. 2 |** *Time series of global model parameters.*** *a,* Temporal evolution of the number of plates *N*<sub>plate</sub> and the surface of deformed areas *S*<sub>deformed</sub> (*S*<sub>tot</sub>, total surface). *b–d,* Time series of boundary lengths (plates and blocks), number of mantle plumes detected in the upper mantle (after ref. 49, Model 5) and net rotation rate of the lithosphere<sub>MM</sub>*ω*<sub>litho</sub> relative to a mean mantle reference frame,

<figure><img src="image_8.png" /></figure>

respectively. *e,* Time series of the mean L2 norm of the velocity in a mean mantle reference frame $\langle\|\mathbf{v}_{MM}\|\rangle$ at six different depths. Vertical orange dashed lines mark $t_{ref}$ for the GPRE and the RPRE at 252 Myr. The light orange area indicates the temporal extent of the GPRE.

plates between 124 and 150 Myr. The detection of such a global plate reorganization event (GPRE) shows that global kinematic changes can arise self-consistently from the dynamics of the coupled mantle-lithosphere system, suggesting that GPREs are likely to occur on Earth, as hypothesized from plate reconstruction models<sup>15,16,18</sup>. On the basis of our plate motion analysis, we identified only one GPRE over the 474 Myr period, suggesting that GPREs, unlike RPREs, are rare. This is consistent with the existence of only two potential candidates for GPREs since the breakup of Pangaea<sup>12,15,18</sup>. The time series of surface motions show that unlike short-lived RPREs, the total duration of the GPRE, defined as the instability period in the velocity of plates, lasts ~20–25 Myr (Fig. 1b). Nevertheless, the motion of the major plates responds synchronously (within 2 Myr) to a peak event at 136 Myr (Fig. 1b), suggesting that the GPRE represents a transition phase between two tectonic stable states culminating in a sharp reorganization event. Thus, the apparent temporal extent proposed for the mid-Cretaceous (~111–100 million years ago (Ma) (refs. 12,13)) and mid-Eocene (~53–40 Ma (refs. 17,18), with a potential main event at 47 Ma (ref. 26)) GPREs might relate to the processes leading to a GPRE rather than being uncertainties in the timing of its surface expression.

## Reorganizations of the plate network

In the convection model, the GPRE identified from plate kinematics corresponds to a local maximum in the number of plates (Fig. 2a). Its onset is marked by a rapid increase in the length of boundaries, especially converging ones (Fig. 2b), whereas its end coincides with a decrease in the length of the boundaries (Fig. 2b and Extended Data Fig. 4a). The abrupt decrease in the surface of deformed areas (Fig. 2a) highlights a regularization of the localization of the deformation into narrower and more stable plate boundaries. No such change in the number of plates nor the length of the boundaries was observed for the 252 Myr RPRE (Fig. 2a,b). This observation suggests that GPREs are linked to a reorganization of the plate network structure by the formation (plate fragmentation) and the destruction (plate coalescence) of plate boundaries.

However, global parameters such as plate number or total length of the boundaries may be biased by intense local readjustments, questioning their interpretation at a global scale. To investigate more deeply the evolution of the plate network through the formation and destruction of plate boundaries in space and time, we computed a new type of graph representing the lineage of the tectonic plates in the convection model (Fig. 3, Methods and Extended Data Figs. 5–9). This representation shows that the plate network can be described between times *t* and *t* + *t* as a set of plate groups (Methods), each group being surrounded by stable continuous plate boundaries. On the temporal evolution graph of the plate network (Fig. 3), these groups take the form of spindles in which plates exchange fragments only with other plates within the group through fragmentation and coalescence. For example, between 220–278 Myr, the two plates that originally carried the synthetic hotspot tracks F and G (Fig. 3a) are part of the same group and ultimately merge to form a larger plate. This group is disjointed from the group containing the large plate carrying the synthetic hotspot tracks E and H, although the two groups share several common boundaries (Fig. 1a).

On the temporal plate graph, the RPRE at 252 Myr shows relative stability of the global plate network, the main event being the destabilization of the plate group containing the fragmenting plate (Fig. 3a). On the opposite, the GPRE is characterized by interactions between most plate groups, forming a knot in the graph (Fig. 3b) that coincides with the reorganization of the plate velocities (Fig. 1b). The network of plates becomes highly interconnected during this event with plates exchanging fragments through the creation and the destruction of plate boundaries, confirming here that the entire plate network reorganizes globally. Furthermore, our results show that the organization of the plate network during a GPRE is a transient state that leads to a new stable state in response to a change in the force balance at a global scale.

## Mantle flow reorganizations

Quantifying the dynamic interactions between Earth's mantle and lithosphere is challenging due to our poor knowledge of the history
<figure><img src="image_9.png" /><figcaption>**Fig. 3 | Graph of the temporal evolution of the plate network in the model.**
a, Graph centred on the 252-Myr RPRE. The red circle indicates the plate which is fragmenting, at the origin of the RPRE. b, Graph centred on the GPRE. Time moves forward from top to bottom. Plates are represented by either a pentagon if their continental surface is greater than 1/3 of their total surface or by a circle otherwise. The size of the node scales with the area of the plate. Seven present-day Earth plates are given as reference (PA, Pacific; NU, Nubia; AU, Australia; NZ, Nazca; AR, Arabia; CO, Cocos; JF, Juan de Fuca). Letters from A to I indicate and localize on the graphs the origins of synthetic hotspot tracks of Fig. 1.</figcaption></figure>

Several lines starting from a plate represent a fragmentation. Several lines converging on a plate represent a fusion. In a perfectly static configuration, a plate at the time *t* is connected vertically only to two plates: itself at the time *t*−1 and *t* + 1. Plates at *t* = *t*<sub>ref</sub> are emphasized by a thicker edge line. Plates at time *t* are coloured according to their velocity *v*<sub>t</sub> relative to their velocity *v*<sub>tref</sub> at *t* = *t*<sub>ref</sub>. A plate with no connection at *t* = *t*<sub>ref</sub> is shown in yellow (only small plates). Velocities are in the mean mantle reference frame (no net rotation of the mantle). For clarity, isolated small plates have been removed from the graphs. Supplementary Fig. 1 provides the full graph.

of Earth's mantle²⁷,²⁸. Our convection model, with its self-organization of mantle flow and plate tectonics, provides a unique window on mantle–lithosphere interactions.

During the detected GPRE, the mean velocity norm at various depths mirrors the first-order evolution of plate velocities (Fig. 2e). Globally, the onset of the GPRE is marked by a velocity increase, peaking at 136 Myr, followed by a sharp decline (~20% at the surface; Fig. 2e). The average second invariant of the strain rate evolves similarly (Extended Data Fig. 4b), indicating mantle deformation and flow reconfiguration. This mantle dynamic reconfiguration during the GPRE coincides with a 40% reduction in net lithosphere rotation after its peak during the main GPRE event (Fig. 2d). Thus, the global reorganization of the plate network during the GPRE coincides with a major reorganization of mantle dynamics from the surface to its base. Conversely, the RPRE at 252 Myr does not correlate with global

velocity change, highlighting that the surface extent of a GPRE reflects its depth extent.

The GPRE is followed by a 50% increase in the number of plumes at 150 Myr (Fig. 2c). We interpret this plume pulse as the result of an increase in vigour of the return flow triggered by a slab superflux during the GPRE, visible in the increase of the length of convergent plate boundaries (Fig. 2b), as proposed by ref. 29 for the mid-Cretaceous plume pulse 130 Ma. The delay of ~14 Myr between the main GPRE event and the onset of the plume pulse is consistent with the ascent of hot buoyant material from the base of the mantle at an average velocity of ~21 cm yr⁻¹, similar to Earth estimates between 23 and 54 cm yr⁻¹ (refs. 30,31).

The synchronous response over the whole mantle through a GPRE provides new arguments to re-evaluate the nature of the mid-Eocene event around 47 Ma. The change in the motion of the Pacific plumes
<figure><img src="image_10.png" /></figure>

**Fig. 4** | **Three-dimensional views of the model over the GPRE.** **a–c,** Snapshots of the model 10 Myr before the onset of the GPRE (**a**, 114 Myr), during the GPRE (**b**, *t*<sub>ref</sub> 136 Myr) and 10 Myr after the end of the GPRE (**c**, 160 Myr). The temperature (on the sphere) and the horizontal velocities (closed arrows) are represented at 580-km depth. Mantle plumes and cold slabs are visualized with two isotherms. Surface plate boundaries are represented by light grey lines and the extent of continents by transparent brown surfaces. Open arrows indicate

the following highlighted features: (1) region of subduction cessation: coalescence of several subduction fronts, (2) subduction propagation, (3) formation of a long semi-continuous convergent boundary, (4) subduction inversion. The orientation of the 3D views is the same for all three snapshots and is fixed relative to a NNR reference frame. The inset surface map in **b** (same projection as Fig. 1) locates with a blue line the central meridian of the 3D view.

relative to each other³²⁻³⁵ and to the Indo-Atlantic plumes³⁶⁻³⁸ has been interpreted as a main argument against a global plate motion change during the mid-Eocene. Indeed, the observed kinks in the Pacific hotspot tracks do not reflect only the motion of the Pacific plate relative to a fixed mantle but also the motion of plumes themselves relative to the mantle³²⁻³⁵. Furthermore, absolute plate motion (APM) models based solely on the Indo-Atlantic plumes do not predict any notable change in plate motion around 47 Ma when reprojected on the Pacific³⁷,³⁹. However, these Indo-Atlantic APM models accurately predict the Pacific hotspot tracks within the uncertainty ellipses from the present to approximately 47 Ma, but begin to deviate substantially from the actual Pacific hotspot tracks around 47 Ma. Besides, our analysis of the GPRE in the convection model shows that the mantle responds synchronously and globally to a GPRE. Thus, in contrast to previous studies³⁷, during a GPRE we expect a radical change in the predictions of hotspot tracks from APM models derived from distinct hotspot populations (the Indo-Atlantic and Pacific hotspots being above two different mantle domains⁴⁰,⁴¹). We propose that the reorganization of plate motion during the mid-Eocene was accompanied by a reorganization of mantle flow, explaining the loss of coherence between Pacific and Indo-Atlantic APM models around chron 21 and the synchronous motion change for individual plumes such as Hawaii-Emperor.

## Origins of global kinematic changes

The timescales for both regional and global events are too short to correspond to a convective overturn⁴². However, lithospheric-scale events capable of triggering strain localization processes such as slab detachment, subduction initiation or cessation or plate fragmentation, are considered as potential mechanisms for Earth RPRES⁶,¹⁰,⁴³⁻⁴⁶, similar, for instance, to the model's RPRE at 252 Myr (Fig. 1), following a plate fragmentation (Fig. 3 and Extended Data Fig. 3). On the other hand, the mechanisms capable of triggering a global plate motion change remain unclear. A drastic change in plate boundary forces has been a favoured scenario¹⁶,⁴⁷ but remains speculative in the absence of observations or modelling results on a global scale¹⁸,⁴²,⁴⁷,⁴⁸.

The time series (Fig. 2) and the temporal plate graph (Fig. 3) show that the GPRE in the model is linked to the creation and destruction of plate boundaries, especially convergent plate boundaries. Snapshots of the model over the GPRE (Fig. 4) show that it coincides with the development of a hemispheric-scale convergent plate boundary separating the main continental fragments (between 50 and 120 °E; Fig. 4 and Extended Data Fig. 3). This plate boundary evolves in just -25 Myr from a series of individual and isolated local subduction zones to a long and semi-continuous convergence zone. This reconfiguration leads to a change in the forces acting on adjacent plates, as illustrated by the rapid rotation of the velocity vectors, accompanied by local subduction initiations/propagations and cessations (Fig. 4). Similarly, an increase in subduction zone initiation is suspected on Earth around the mid-Eocene event⁴⁵,⁴⁶, particularly in the Pacific Ocean, indicating changes in the balance of plate forces, with subduction initiations possibly both driving and resulting from the plate reorganization. We do not observe a spatial propagation of the reorganization during the GPRE identified in the model. This suggests a strong mechanical coupling between the components of the mantle-lithosphere system, leading to a global and instantaneous (<2 Myr) interconnection between them, contrary to the chain reaction mechanism proposed by ref. 14.

We propose here that rapid changes in plate boundary forces at this long convergent boundary are then transmitted to the entire plate network through both the rigid lithosphere and the underlying mantle. The observation of synchronous reorganization of mantle flow and surface motions indicates that the change in mantle dynamics is not purely passive but has accompanied the change in plate motion and plays a role in up-scaling the kinematic change from local to global. This role of the mantle in the development of the GPRE coupled with forces transmitted through the rigid lithosphere, provides a possible mechanical explanation for the link between the India-Eurasia collision and the Tonga-Kermadec and Izu-Bonin-Mariana subduction initiations in the Western Pacific Ocean just before 47 Ma (ref. 26), suggesting with the arguments we presented on the APM models that the mid-Eocene event is likely to be a GPRE.
Our results show that the temporal evolution of plate motion follows a dynamic equilibrium punctuated by phases of rapid transitions. Plates respond to a set of forces acting on them. In most cases, far-field forces are not able to destabilize plates or their boundaries, and plate motion remains constant for long periods. Sporadically, plates reorganize regionally due to local changes in plate boundary forces. Only during rare events, the GPREs, the entire Earth's plate network reorganizes globally together with mantle flow. These events represent a transient phase of the plate motion and plate network evolution between two more stable states. A drastic change in boundary forces at hemispheric- or global-scale convergent plate boundary is likely to trigger a GPRE with forces instantaneously (<2 Myr) transmitted through both the lithosphere and the mantle to the entire mantle-lithosphere system.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-025-01772-7.

## References

1. DeMets, C., Gordon, R. G., Argus, D. & Stein, S. Current plate motions. *Geophys. J. Int.* **101**, 425–478 (1990).

2. Gordon, R. G. & Stein, S. Global tectonics and space geodesy. *Science* **256**, 333–342 (1992).

3. DeMets, C., Gordon, R. G. & Argus, D. F. Geologically current plate motions. *Geophys. J. Int.* **181**, 1–80 (2010).

4. Heirtzler, J. R., Dickson, G., Herron, E., pitman III, W. & Le Pichon, X. Marine magnetic anomalies, geomagnetic field reversals, and motions of the ocean floor and continents. *J. Geophys. Res.* **73**, 2119–2136 (1968).

5. Menard, H. & Atwater, T. Changes in direction of sea floor spreading. *Nature* **219**, 463–467 (1968).

6. Iaffaldano, G. & Bunge, H.-P. Rapid plate motion variations through geological time: observations serving geodynamic interpretation. *Annu. Rev. Earth Planet. Sci.* **43**, 571–592 (2015).

7. Bercovici, D., Tackley, P. & Ricard, Y. in *Treatise on Geophysics* 2nd edn (ed Schubert, G.) 271–318 (Elsevier, 2015).

8. DeMets, C., Iaffaldano, G. & Merkouriev, S. High-resolution neogene and quaternary estimates of nubia-eurasia-north america plate motion. *Geophys. J. Int.* **203**, 416–427 (2015).

9. DeMets, C. & Merkouriev, S. High-resolution reconstructions of South America plate motion relative to Africa, Antarctica and North America: 34 ma to present. *Geophys. J. Int.* **217**, 1821–1853 (2019).

10. DeMets, C. & Merkouriev, S. Detailed reconstructions of india-somalia plate motion, 60 ma to present: implications for somalia plate absolute motion and india-eurasia plate motion. *Geophys. J. Int.* **227**, 1730–1767 (2021).

11. DeMets, C., Merkouriev, S. & Sauter, D. High resolution reconstructions of the southwest Indian Ridge, 52 ma to present: implications for the breakup and absolute motion of the Africa plate. *Geophys. J. Int.* **226**, 1461–1497 (2021).

12. Matthews, K. J., Seton, M. & Müller, R. D. A global-scale plate reorganization event at 105–100 ma. *Earth Planet. Sci. Lett.* **355**, 283–298 (2012).

13. Olierook, H. K. et al. Timing and causes of the mid-cretaceous global plate reorganization event. *Earth Planet. Sci. Lett.* **534**, 116071 (2020).

14. Gürer, D., Granot, R. & van Hinsbergen, D. J. Plate tectonic chain reaction revealed by noise in the cretaceous quiet zone. *Nat. Geosci.* **15**, 233–239 (2022).

15. Rona, P. A. & Richardson, E. S. Early cenozoic global plate reorganization. *Earth Planet. Sci. Lett.* **40**, 1–11 (1978).

16. Patriat, P. & Achache, J. India-Eurasia collision chronology has implications for crustal shortening and driving mechanism of plates. *Nature* **311**, 615–621 (1984).

17. Whittaker, J. et al. Major Australian-Antarctic plate reorganization at Hawaiian-Emperor Bend time. *Science* **318**, 83–86 (2007).

18. 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).

19. Gaina, C. & Jakob, J. Global Eocene tectonic unrest: possible causes and effects around the North American plate. *Tectonophysics* **760**, 136–151 (2019).

20. Edelsbrunner, H. & Harer, J. *Computational Topology: An Introduction* (AMS, 2009).

21. Zomorodian, A. J. Topology for computing. in *Algorithms and Theory of Computation Handbook* 2nd edn (eds Atallah, M. J. & Blanton, M.) 82–112 (CRC Press, 2010).

22. Carlsson, G. & Vejdemo-Johansson, M. *Topological Data Analysis with Applications* (Cambridge Univ. Press, 2022).

23. Coltice, N., Husson, L., Faccenna, C. & Arnould, M. What drives tectonic plates? *Sci. Adv.* **5**, eaax4295 (2019).

24. Forsyth, D. & Uyeda, S. On the relative importance of the driving forces of plate motion. *Geophys. J. Int.* **43**, 163–200 (1975).

25. Lithgow-Bertelloni, C. & Richards, M. A. Cenozoic plate driving forces. *Geophys. Res. Lett.* **22**, 1317–1320 (1995).

26. Gaastra, K. M., Gordon, R. G. & Woodworth, D. T. Quantification of pacific plate hotspot tracks since 80 ma. *Tectonics* **41**, e2021TC006772 (2022).

27. Conrad, C. P., Steinberger, B. & Torsvik, T. H. Stability of active mantle upwelling revealed by net characteristics of plate tectonics. *Nature* **498**, 479–482 (2013).

28. Coltice, N. & Shephard, G. E. Tectonic predictions with mantle convection models. *Geophys. J. Int.* **213**, 16–29 (2018).

29. East, M., Müller, R. D., Williams, S., Zahirovic, S. & Heine, C. Subduction history reveals cretaceous slab superflux as a possible cause for the mid-Cretaceous plume pulse and superswell events. *Gondwana Res.* **79**, 125–139 (2020).

30. Poore, H., White, N. & Jones, S. A neogene chronology of iceland plume activity from V-shaped ridges. *Earth Planet. Sci. Lett.* **283**, 1–13 (2009).

31. Turcotte, D. & Schubert, G. *Geodynamics* (Cambridge Univ. Press, 2014).

32. Tarduno, J. A. et al. The emperor seamounts: southward motion of the Hawaiian hotspot plume in Earth's mantle. *Science* **301**, 1064–1069 (2003).

33. Tarduno, J., Bunge, H.-P., Sleep, N. & Hansen, U. The bent Hawaiian-Emperor hotspot track: inheriting the mantle wind. *Science* **324**, 50–53 (2009).

34. Bono, R. K., Tarduno, J. A. & Bunge, H.-P. Hotspot motion caused the hawaiian-emperor bend and llsvps are not fixed. *Nat. Commun.* **10**, 3370 (2019).

35. Hu, J., Gurnis, M., Rudi, J., Stadler, G. & Müller, R. D. Dynamics of the abrupt change in Pacific plate motion around 50 million years ago. *Nat. Geosci.* **15**, 74–78 (2022).

36. Doubrovine, P. V., Steinberger, B. & Torsvik, T. H. Absolute plate motions in a reference frame defined by moving hot spots in the Pacific, Atlantic, and Indian Oceans. *J. Geophys. Res. Solid Earth* **117**, B09101 (2012).

37. Maher, S., Wessel, P., Müller, R., Williams, S. & Harada, Y. Absolute plate motion of Africa around Hawaii-Emperor Bend time. *Geophys. J. Int.* **201**, 1743–1764 (2015).

38. Koivisto, E. A., Andrews, D. L. & Gordon, R. G. Tests of fixity of the Indo-Atlantic hot spots relative to Pacific hot spots. *J. Geophys. Res. Solid Earth* **119**, 661–675 (2014).
39. O'Neill, C., Müller, 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).

40. Burke, K., Steinberger, B., Torsvik, T. H. & Smethurst, M. A. Plume generation zones at the margins of large low shear velocity provinces on the core–mantle boundary. *Earth Planet. Sci. Lett.* **265**, 49–60 (2008).

41. 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).

42. Richards, M. A. & Lithgow-Bertelloni, C. Plate motion changes, the Hawaiian-Emperor Bend, and the apparent success and failure of geodynamic models. *Earth Planet. Sci. Lett.* **137**, 19–27 (1996).

43. Faccenna, C., Becker, T. W., Lallemand, S. & Steinberger, B. On the role of slab pull in the cenozoic motion of the pacific plate. *Geophys. Res. Lett.* **39**, L03305 (2012).

44. Domeier, M. et al. Intraoceanic subduction spanned the pacific in the late Cretaceous–Paleocene. *Sci. Adv.* **3**, eaao2303 (2017).

45. Crameri, F. et al. A transdisciplinary and community-driven database to unravel subduction zone initiation. *Nat. Commun.* **11**, 3750 (2020).

46. Lallemand, S. & Arcay, D. Subduction initiation from the earliest stages to self-sustained subduction: Insights from the analysis of 70 Cenozoic sites. *Earth Sci. Rev.* **221**, 103779 (2021).

47. Lithgow-Bertelloni, C. & Richards, M. A. The dynamics of Cenozoic and Mesozoic plate motions. *Rev. Geophys.* **36**, 27–78 (1998).

48. Patriat, P., Sloan, H. & Sauter, D. From slow to ultraslow: a previously undetected event at the southwest Indian Ridge at ca. 24 ma. *Geology* **36**, 207–210 (2008).

49. Arnould, M., Coltice, N., Flament, N. & Mallard, C. Plate tectonics and mantle controls on plume dynamics. *Earth Planet. Sci. Lett.* **547**, 116439 (2020).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. 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 the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025
# Methods

## Mantle convection modelling

In this study, we analyse the three-dimensional, spherical mantle convection model developed by ref. 23 to investigate the dynamic interactions between Earth's tectonic plates and the underlying mantle during plate reorganizations. This model is uniquely designed to push the limits of self-consistent plate generation by incorporating Earth-like convective vigour, a pseudo-plastic rheology and deformable continental rafts. It further accounts for extreme viscosity gradients—necessitating high-resolution simulations (average surface horizontal resolution of 23 km and less than 4,000 years per time step) and a weak, lubricating upper oceanic crust, both features being central in generating one-sided subduction (Extended Data Fig. 1). The deep structure of the modelled mantle includes a viscosity jump at 660 km by a factor of 30, as suggested by geoid⁵⁰ and post-glacial rebound studies⁵¹. The combined effect of this feature and the strong temperature dependence generate self-consistently a low viscosity asthenosphere. In the model, both deep dense material at the core–mantle boundary and a decrease of the thermal expansivity structures the plume network and dynamics⁴⁹.

## Plate rigidity criterion

Plate tessellation requires an objective and quantitative diagnostic to assess the rigidity of a piece of the sphere and determine whether it can be interpreted as a tectonic plate.

Plate tectonics provides a framework for the study of plate dynamics and theorizes that the motion of each plate consists of a single rigid-body rotation. We solve a linear inverse problem to find the optimal Euler pole of all the plates of a given tessellation and the associated misfit. We define two different measures of the stiffness of the plate *i*, the plateness $P_{1,i}$ and $P_{2,i}$, following the formalism of ref. 52:

$$P_{1,i}(\boldsymbol{\omega}) = \frac{1}{N} \sum_{k=1}^{N} \frac{\mathbf{v}_{p,k}(\boldsymbol{\omega}) \cdot \mathbf{v}_{d,k}}{\| \mathbf{v}_{p,k}(\boldsymbol{\omega}) \| \| \mathbf{v}_{d,k} \|} \quad (1)$$

$$P_{2,i}(\boldsymbol{\omega}) = 1 - \frac{1}{N} \sum_{k=1}^{N} \frac{\| \mathbf{v}_{p,k}(\boldsymbol{\omega}) - \mathbf{v}_{d,k} \|}{\| \mathbf{v}_{d,k} \|} \quad (2)$$

where $\mathbf{v}_{d,k}$ are the $N$ data velocity vectors at the surface of the plate *i* and $\mathbf{v}_{p,k}(\boldsymbol{\omega})$ are their predictions computed from the angular velocity vector $\boldsymbol{\omega}$ (the operator $\| \cdot \|$ corresponds here to the L2 norm). Plateness $P_1$ is a measure of the good alignment between the velocities predicted by a single rigid-body rotation and the data. Plateness $P_2$ quantifies the magnitude of the difference between the predicted velocities and the data. Taken together, $P_1$ and $P_2$ give an objective measurement of the stiffness of a plate, that is, the ability of its surface velocity field to be described by a single rigid-body rotation.

Following the above formulation, we compute for each plate *i* the plateness $P_{1,i}(\boldsymbol{\omega}_i)$ and $P_{2,i}(\boldsymbol{\omega}_i)$ where $\boldsymbol{\omega}_i$ is the optimal Euler pole describing the velocity field of the plate *i*. However, as plateness defined by ref. 52 are average properties at the scale of a plate, they are not designed to diagnose a possible high local misfit. To test for spatial variations in plateness, we decompose each plate *i* into *M* fragments *j* with a characteristic size *S* and we evaluate the local plateness $P_{1,j}(\boldsymbol{\omega}_i)$ and $P_{2,j}(\boldsymbol{\omega}_i)$ for each fragment *j*. In this study, we define a plate *i* as rigid if:

$$\left\{ \begin{array}{l} P_{1,i}(\boldsymbol{\omega}_i) \ge P_{1c} \\ P_{2,i}(\boldsymbol{\omega}_i) \ge P_{2c} \\ P_{1,i,j}(\boldsymbol{\omega}_i) \ge P_{1c} \quad \forall j \in [1, M] \\ P_{2,i,j}(\boldsymbol{\omega}_i) \ge P_{2c} \quad \forall j \in [1, M] \end{array} \right. \quad (3)$$

where $P_{1c}$ and $P_{2c}$ are critical plateness. The size *S* of the plate fragments represents the spatial resolution of our rigidity criterion, here

$S = 4.15 \times 10^6 \text{ km}^2$ (S - 0.81% of Earth's surface. We choose the values of critical plateness according to a comparison with Earth. We applied our tessellation workflow (described below) on the present-day Earth surface velocity field model from GSRM v.2.1⁵³ and defined $P_{1c} = 0.90$ and $P_{2c} = 0.80$ according to our tessellation, which best fits the plate models commonly used for Earth.⁵⁴ (Extended Data Figs. 2b,d and 5). Extended Data Fig. 5 shows that the tessellation of the Earth velocity model with a stronger rigidity criterion ($P_{1c} = 0.90, P_{2c} = 0.90$) leads to a higher sensitivity to small-scale variations in the deformed areas and an overfit of the imposed rigid polygons in GSRM v.2.1.

The plateness we defined and use (equations (1) and (2)) are designed to be (1) used at the scale of a plate and not globally and (2) be used on present-day Earth plate to generalize its usage. Thus, our plateness differs from those previously introduced in the literature⁵⁵,⁵⁶, designed as a global diagnostic to compare numerical models with each other and not applicable on Earth.

## Tessellation

Whether for the Earth or a 3D mantle convection model, the plate tessellation method we propose in this study follows the same procedure in two main steps: (1) the automatic detection of polygon (potentially plate) boundaries based on the segmentation of the norm of the surface velocity gradient $\|\nabla\mathbf{v}_s\|$ using topological data analysis methods and (2) the adaptive selection of polygon boundaries to maximize the rigidity of the detected polygons and the uniqueness of their solid body rotation.

The automatic detection of polygon boundaries (step 1) is based on the computation of the Morse complex⁵⁷,⁵⁸ of $\|\nabla\mathbf{v}_s\|$ after its topological simplification⁵⁹,⁶⁰, according to a minimum persistence threshold $p_{\min}$ on its pairs of critical points. In this framework, the polygon boundaries are the integral lines linking local maxima with saddle points of $\|\nabla\mathbf{v}_s\|$ (ascending Morse complexes), that is, the crest lines of $\|\nabla\mathbf{v}_s\|$ (Extended Data Fig. 6). We use the visualization software Paraview⁶¹ (v.5.10.0) and the open-source software platform TTK (Topology ToolKit⁶²,⁶³, v.1.1) to process this first step according to the following workflow (Paraview filters are given in quotes):

a. Load and extract the surface of the velocity field $\mathbf{v}_s$ (Extended Data Fig. 6a, filter 'ExtractSurface').

b. Extract the outer connected component of the sphere surface (filters 'Connectivity' and 'Threshold') and apply filters to avoid further mesh issues (filters 'CleantoGrid' and 'Tetrahedralize').

c. Compute the L2 norm of $\mathbf{v}_s$ (Extended Data Fig. 6b, filters 'Gradient' and 'Calculator').

d. Find all the critical points of $\|\nabla\mathbf{v}_s\|$ and compute the topological persistence of each pair (Extended Data Fig. 6c, filter 'TTKPersistenceDiagram').

e. Apply a threshold to the topological persistence to keep only the pairs of critical points whose persistence *p* is larger than a given minimum $p \ge p_{\min}$, (Extended Data Fig. 6d, filter 'Threshold').

f. Compute the topological simplification of $\|\nabla\mathbf{v}_s\|$ according to the reduced persistence diagram obtained after step e. (Extended Data Fig. 6e, filter 'TTKTopologicalSimplification'). The comparison between zoom panels of the Extended Data Fig. 6b,e emphasizes an area where these fine changes can be appreciated.

g. Compute the ascending Morse complex of the simplified $\|\nabla\mathbf{v}_s\|$ (Extended Data Fig. 6f, filter 'TTKMorseSmaleComplex').

The parameter $p_{\min}$ conditions the detection limit of the diffuse plate boundaries and acts as a noise-reduction parameter. The higher $p_{\min}$, the higher the persistence of critical point pairs must be to be considered as belonging to a polygon boundary. As a result, the higher $p_{\min}$ is, the fewer polygons there are (Extended Data Fig. 7).
The aim of the step 1 of our plate tessellation process is to compute the segmentation of a given surface velocity field for *N* different values of *p*<sub>min</sub> (from *p*<sub>min,i=0</sub> to *p*<sub>min,i=N</sub>, the minimum and the maximum values considered, respectively) whose effect on the surface segmentation ranges from over- to under-tessellation (Extended Data Fig. 7). In this study, we choose *N* = 24 to fine tesselate the intermediate states (Extended Data Figs. 7 and 8b).

Step 2 consists in combining segmentations with different values of topological persistence to obtain a tessellation where *p*<sub>min</sub> varies laterally. Extended Data Fig. 7 illustrates the following workflow:

a. Select a random point **A** on the sphere (red cross on Extended Data Fig. 7).

b. Select the tessellation computed at the highest minimum persistence value $p_{min,i=N}$.

c. Inverse the velocity field of the polygon $\Gamma_i$ carrying the point **A** for the tessellation computed at $p_{min,i}$.

d. Test if $\Gamma_i$ is rigid according to the rigidity criteria defined above (Methods, 'Plate rigidity criterion').

d1. If yes, then save the geometry of $\Gamma_i$ and return to a. (excluding areas already mapped) until the entire surface of the sphere has been mapped.

d2. Else, return to step c. decreasing the minimum persistence value to $p_{min,i=i-1}$.

The resulting tessellation map (Extended Data Fig. 8a) corresponds to a set of polygons whose geometry is defined at the highest value of $p_{min}$ for which they can be defined as rigid. This step ensures that the smallest number of rigid polygons bounded by velocity gradients is obtained, thus constraining the uniqueness of their solid body rotation. The values of $p_{min,i=0}$ and $p_{min,i=N}$ are chosen to minimize the area of polygons defined at these values on the final tessellation maps (located in the distribution queue; Extended Data Fig. 8b).

## Distinction between plates and deformable blocks

The result of the tessellation is a set of polygons mapping the surface of the sphere. We make a clear distinction here between plates (rigid entities) and blocks (deformable) by analogy with Earth, based on a threshold of strain rate within polygons. A polygon is considered deformable (orange in Fig. 1 and Extended Data Figs. 2 and 5) if more than 90% of its surface has a strain rate greater than $5 \times 10^{-9} \text{ yr}^{-1}$, a threshold corresponding to a low strain rate limit for Earth diffuse plate boundaries in GSRM v.2.1<sup>53</sup>. We use the same threshold to compute the deformed surface area (Fig. 2a) from the second invariant of the strain rate at the surface of the model.

## Time tracking

The time tracking of polygons is a separate step from tessellation. The time tracking is computed once all the tessellation maps have been calculated for all the time steps of interest. We track polygons over time by tracking their centroids defined as their centres of figure embedded in their surface. The Extended Data Fig. 9 illustrates our workflow. At each time step, we derive a persistence diagram pairing for each polygon its centroid (maxima) with the point belonging to its edge that is furthest from the polygon's centroid. The persistence of the pair is defined as the haversine distance between the two points of the pair. Centroids are tracked using a Lifted Wasserstein Matcher<sup>64</sup> (filter 'TTKTrackingFromPersistenceDiagrams' in TTK).

## Plate graphs

Plate graphs (for example, Fig. 3) are directed acyclic graphs that illustrate the evolution and the lineage of the plate network in space and time (excluding deformable blocks). A plate graph $G(V, E)$ connects a set of vertices *V* representing here all the plates at the surface of the model over a given period *T* by a set of edges *E*. We represent two plates

$\Gamma_A$ and $\Gamma_B$ defined at *t* and *t* + *δt* ((*t*, *t* + *δt*) ∈ *T*<sup>2</sup>), respectively, as linked if their surfaces overlap at *t* + *δt* over at least a fraction *κ* of their respective areas. This condition can be written as follows:

$$ \min \left( \frac{S_{\Gamma_A}}{S_{\Gamma_B}}, \frac{S_{\Gamma_B}}{S_{\Gamma_A}} \right) \geq \kappa \quad (4) $$

where $S_{\Gamma_A}$ is the area of overlapping (that is, $S_{\Gamma_A} = S(\Gamma_A \cap \Gamma_B)$) and $S_{\Gamma_A}$ and $S_{\Gamma_B}$ are the areas of plates $\Gamma_A$ and $\Gamma_B$, respectively. $S_{\Gamma_A}$ is corrected for the rotation of $\Gamma_A$ between *t* and *t* + *δt* considering that the kinematics of $\Gamma_A$ at *t* is constant over *δt*. $\kappa$ is a graphical parameter that controls the density of the displayed connection allowing to keep only the most probable plate connections. The larger $\kappa$ is, the fewer connections will be displayed on the graph. Here we choose to display connections between plates sharing at least 40% of their surfaces between two adjacent time steps ($\kappa$ = 40%) to graphically optimize the number of connections displayed, between highlighting the main structure of the plate graph and displaying the secondary connections.

We used the open-source graph visualization software Graphviz<sup>65</sup> and the graph description language DOT<sup>66</sup> for the graphical representation. Plates are organized along the *x* axis to minimize the intersections of the graph's edges. For the clarity of the graph, the bending of each edge is weighted by a coefficient proportional to the squared root of the area of the smallest plates they linked, to graphically favour vertical connections between two large plates. The full plate graph<sup>67</sup> of the analysis of the convection model over 474 Myr is given in Supplementary Fig. 1.

## Plate groups

We define here the set of plate groups between the times *t* and *t* + *τ* as the set of induced subgraphs formed from a larger plate graph $G(V, E)$ keeping only the vertices (that is, plates) between the time *t* and *t* + *τ*. A plate group is a set of vertices of $G(V, E)$ that are connected by edges between the time *t* and *t* + *τ*.

## Data availability

The full plate graph is available via Zenodo at https://doi.org/10.5281/zenodo.15863930 (ref. 67).

## Code availability

The topological data analysis has been performed using the visualization software Paraview (https://www.paraview.org/) and the open-source software platform TTK (Topology ToolKit, v.1.1, https://topology-tool-kit.github.io/). The numerical toolkit MAPT<sup>3</sup> (Multi-disciplinary and Automatic Plate Tessellation and Time tracking Toolkit) developed to produce the results of this paper is available via Github at https://github.com/AlexandrePFJanin/MAPT3. The mantle convection code StagYY is the property of P. Tackley and Eidgenössische Technische Hochschule (ETH) Zürich. It is available on request from P. Tackley (paul.tackley@erdw.ethz.ch). Minor edits to the code can be provided by N.C. (nicolas.coltice@geoazur.unice.fr) on request.

## References

50. Ricard, Y., Richards, M., Lithgow-Bertelloni, C. & Le Stunff, Y. A geodynamic model of mantle density heterogeneity. *J. Geophys. Res. Solid Earth* **98**, 21895–21909 (1993).

51. Kaufmann, G. & Lambeck, K. Mantle dynamics, postglacial rebound and the radial viscosity profile. *Phys. Earth Planet. Int.* **121**, 301–324 (2000).

52. Alisic, L., Gurnis, M., Stadler, G., Burstedde, C. & Ghattas, O. Multi-scale dynamics and rheology of mantle flow with plates. *J. Geophys. Res. Solid Earth* **117**, B10402 (2012).

53. Kreemer, C., Blewitt, G. & Klein, E. C. A geodetic plate motion and global strain rate model. *Geochem. Geophys. Geosyst.* **15**, 3849–3889 (2014).
54. Bird, P. An updated digital model of plate boundaries. *Geochem. Geophys. Geosyst.* **4**, 1027 (2003).

55. Weinstein, S. A. & Olson, P. L. Thermal convection with non-newtonian plates. *Geophys. J. Int.* **111**, 515–530 (1992).

56. Tackley, P. J. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations 1. pseudoplastic yielding. *Geochem. Geophys. Geosyst.* **1**, 1021 (2000).

57. Milnor, J. *Morse Theory* (Princeton University Press, 1963).

58. Robins, V., Wood, P. J. & Sheppard, A. P. Theory and algorithms for constructing discrete morse complexes from grayscale digital images. *IEEE Trans. Pattern Anal. Mach. Intell.* **33**, 1646–1658 (2011).

59. Edelsbrunner, H., Letscher, D. & Zomorodian, A. Topological persistence and simplification. In *Proc. 41st Annual Symposium on Foundations of Computer Science*, 454–463 (IEEE, 2000).

60. Tierny, J. & Pascucci, V. Generalized topological simplification of scalar fields on surfaces. *IEEE Trans. Visualization Comput. Graphics* **18**, 2005–2013 (2012).

61. Ahrens, J., Geveci, B. & Law, C. *Paraview: an end-user tool for large data visualization* Technical Report LA-UR-03-1560 (Los Alamos National Laboratory, 2005).

62. Tierny, J., Favelier, G., Levine, J. A., Gueunet, C. & Michaux, M. The topology toolkit. *IEEE Trans. Visualization Comput. Graphics* **24**, 832–842 (2017).

63. Masood, T. B. et al. An overview of the Topology ToolKit. in *Topological Methods in Data Analysis and Visualization VI*. Mathematics and Visualization (eds Hotz, I. et al). https://doi.org/10.1007/978-3-030-83500-2_16 (Springer, 2021).

64. Soler, M., Plainchault, M., Conche, B. & Tierny, J. Lifted Wasserstein matcher for fast and robust topology tracking. In *2018 IEEE 8th Symposium on Large Data Analysis and Visualization (LDAV)* 23–33 (IEEE, 2018).

65. Ellson, J., Gansner, E. R., Koutsofios, E., North, S. C. & Woodhull, G. Graphviz and dynagraph—static and dynamic graph drawing tools. in *Graph Drawing Software*. Mathematics and Visualization (eds Jünger, M. & Mutzel, P.) 127–148 (Springer, 2004); https://doi.org/10.1007/978-3-642-18638-7_6

66. Gansner, E. R., Koutsofios, E., North, S. C. & Vo, K.-P. A technique for drawing directed graphs. *IEEE Trans. Software Eng.* **19**, 214–230 (1993).

67. Janin, A., Coltice, N., Chamot-Rooke, N. & Tierny, J. Geodynamics of a global plate reorganisation from topological data analysis. *Zenodo* https://doi.org/10.5281/zenodo.15863930 (2025).

68. Crameri, F. Scientific colour maps: perceptually uniform and colour-vision deficiency friendly. *Zenodo* https://doi.org/10.5281/zenodo.1243862 (2018).

69. Seton, M. et al. Global continental and ocean basin reconstructions since 200 ma. *Earth Sci. Rev.* **113**, 212–270 (2012).

70. Hasterok, D. et al. New maps of global geological provinces and tectonic plates. *Earth Sci. Rev.* **231**, 104069 (2022).

## Acknowledgements

N.C. received support from a European Research Council (ERC) advanced grant, award number 101142672. J.T. received support from a European Research Council (ERC) consolidator grant, award number 863464 (https://erc-tori.github.io/). We are grateful to J. Escartin, F. Garel, C. Kreemer, L. Le Pourhiet, C. Prigent and D. Sauter for fruitful discussions. We thank C. Doucerain and B. Pauget for their helpful contributions. We used the Scientific Colour Maps⁶⁸ to prevent misleading representation of the data and exclusion of readers with colour-vision deficiencies.

## Author contributions

A.J., N.C. and N.C.-R. designed the study. A.J. and J.T. developed the methodology. A.J. conducted the analysis, produced the figures and took the lead in writing the manuscript. A.J., N.C. and N.C.-R. interpreted the results. All authors contributed to the editing and revision process.

## Competing interests

The authors declare no competing interests.

## Additional information

*Extended data* is available for this paper at
https://doi.org/10.1038/s41561-025-01772-7.

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s41561-025-01772-7.

**Correspondence and requests for materials** should be addressed to
Alexandre Janin.

**Peer review information** *Nature Geoscience* thanks Fabio Crameri and Lukas Fuchs and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Stefan Lachowycz, in collaboration with the *Nature Geoscience* team.

**Reprints and permissions information** is available at
www.nature.com/reprints.
<figure><img src="image_11.png" /></figure>

**Extended Data Fig. 1 | Asymmetrical subductions.** a, Surface map of the model at 190 Myr (same as in Fig. 1b). b–e, Slices in the viscosity field of the model at 190 Myr, along great circles and across four different subduction zones. Slices go from the surface to the base of the mantle. The surface extent of the slices is indicated on a by the four coloured thick lines A-A', B-B', C-C' and D-D'. White

arrows represent the velocity field of the model displayed in the upper plate reference frame and reprojected on the slices. The 'pins' illustrate graphically the plate used as reference. The thick black arrows indicate the location of the trenches detected by the tessellation method discussed in this paper. Stars highlight the 660-km depth discontinuity.
<figure><img src="image_12.png" /></figure>

a

log₁₀(εᵢ) (in log₁₀(yr⁻¹))

<figure><img src="image_13.png" /></figure>

c

Lithosphere:

Cont. Interior

Cont. Belt

Oceanic

Deformed
polygons

<figure><img src="image_14.png" /></figure>

b

log₁₀(εᵢ) (in log₁₀(yr⁻¹))

<figure><img src="image_15.png" /></figure>

d

Major plate/block boundaries:

Convergent

Divergent

Strike-slip

**Extended Data Fig. 2** | Plate tessellation of the three-dimensional spherical mantle convection model and the present-day Earth. **a–b**, NNR surface velocity fields and second invariant of the strain rate tensor (ϵ̇_ll) of the mantle convection model at 96 Myr and the model for the present-day Earth from GSRM v.2.1³³,

respectively. **c–d**, Plate tessellation of the mantle convection model and the present-day Earth, respectively. See Methods for deformed polygons. Continental belts drawn from ref. 69 and interiors from ref. 70.
<figure><img src="image_16.png" /><figcaption>b (252 Myr) ★</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>c (272 Myr)</figcaption></figure>

<figure><img src="image_18.png" /></figure>

Lithosphere:

Cont. Interior

Cont. Belt

Oceanic

Deformed
polygons

<figure><img src="image_19.png" /><figcaption>e (136 Myr) ★</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>f (156 Myr)</figcaption></figure>

<figure><img src="image_21.png" /></figure>

Plate/block
boundaries:

<figure><img src="image_22.png" /></figure>

-1.5 0 1.5

Divergence
(x10⁻¹⁵ s⁻¹)

<figure><img src="image_23.png" /></figure>

0 7

[Vorticity]
(x10⁻¹⁵ s⁻¹)

**Extended Data Fig. 3** | Surface maps of model snapshots. a-c, Model snapshots at three different times: 20 Myr before the discussed RPRE, during the RPRE and 20 Myr after the end of the RPRE, respectively. d-f, Model snapshots at three different times: 20 Myr before the main event of the GPRE (8 Myr before the onset

of the GPRE), during the main GPRE event, and 20 Myr after the main GPRE event (6 Myr after the end of the GPRE), respectively. Black stars indicate $t_{ref}$. Divergence and absolute vorticity fields are displayed below the detected boundaries.
<figure><img src="image_24.png" /></figure>

**Extended Data Fig. 4** | Complementary time series of global model parameters.

a, Temporal evolution of the total length of boundaries *that is* plates and deformed areas (running mean with a time window of 6 Myr). b, Time series of the

<figure><img src="image_25.png" /></figure>

average of the second invariant of the strain rate tensor at six different depths. Vertical orange dashed lines mark *t*<sub>ref</sub> for the GPRE and the RPRE at 252 Myr. The light orange area indicates the temporal extent of the GPRE.
<figure><img src="image_26.png" /><figcaption>**Extended Data Fig. 5** | Present-day Earth plate tessellations. **a**, Favoured tessellation with $P_{1c}$ = 0.90 and $P_{2c}$ = 0.80 (Extended Data Fig. 2d). **b**, Tessellation with a stricter rigidity criterion $P_{1c}$ = 0.90 and $P_{2c}$ = 0.90.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Extended Data Fig. 6** | Sequential tessellation workflow for the surface of a 3-dimensional spherical mantle convection model. See Methods, Tessellation, step 1 for more details. a, Input of the surface velocity field. b, Computation of the norm of the velocity gradient. c, Extraction of all the critical points and computation of the associated topological persistence for each critical points pair. d, Application of a threshold on the persistence. e, Computation of the</figcaption></figure>

topologically simplified norm of the velocity gradient. On the zoom panel, the white arrows emphasized areas where changes can be seen compared to the zoom panel in b. f, Ascending Morse complex showing the resulting tessellation. Polygon boundaries (that is separatrices of the Morse complex) are displayed in white). The operator $|| \cdot ||$ represents the L2 norm, *indicates a normalisation (used for display only) that is $|| v_{surf} ||^* = || v_{surf} || / \max( || v_{surf} ||)$.
<figure><img src="image_28.png" /><figcaption>**Extended Data Fig. 7** | Combination of segmentations with different values of topological persistence. See Methods, Tessellation, step 2 for more details. Illustration of the iterative workflow considering only four different values of $p_{min}$</figcaption></figure>

for clarity. Thick black arrows indicate the sequence of the different steps. Filled polygons are the polygons $\Gamma_i$ carrying the point represented by the red cross at $p_{min,i}$, coloured according to the value of $p_{min,i}$.
<figure><img src="image_29.png" /></figure>

**Extended Data Fig. 8** | Tessellation with lateral variation of $p_{min}$·a, Result of the tessellation after the iterative process illustrated in Extended Data Fig. 7 for the mantle convection model at the same time step as in Extended Data Fig. 2 (t=96 Myr). The polygons are coloured according to the value of $p_{min}$ at which

<figure><img src="image_30.png" /></figure>

they are considered rigid (or the $p_{min,0}$ if the polygon is not rigid at $p_{min,0}$).

b, Histogram showing the distribution of $p_{min}$ values on the surface of the sphere, averaged over 400 Myr. The excess of surface area at $p_{min,0}$ is due to a large continental plate bounded by persistent plate boundaries.
<figure><img src="image_31.png" /><figcaption>**Extended Data Fig. 9** | **Time tracking workflow.** See Methods, Time tracking for more details. **a–b**, Computation of polygon centroids at each time step. **c**, Resulting polygon centroid tracks, shown here only over 20 Myr for clarity. Polygons are then tracked over time by tracking their centroids.</figcaption></figure>