# The convergence history of India-Eurasia records multiple subduction dynamics processes

## Abstract

During the Cretaceous, the Indian plate moved towards Eurasia at the fastest rates ever recorded. The details of this journey are preserved in the Indian Ocean seafloor, which document two distinct pulses of fast motion, separated by a noticeable slowdown. The nature of this rapid acceleration, followed by a rapid slowdown and then succeeded by a second speedup, is puzzling to explain. Using an extensive observation dataset and numerical models of subduction, we show that the arrival of the Reunion mantle plume started a sequence of events that can explain this history of plate motion. The forces applied by the plume initiate an intra-oceanic subduction zone, which eventually adds enough additional force to drive the plates at the anomalously fast speeds. The two-stage closure of a double subduction system, including accretion of an island arc at 50 million years ago, may help reconcile geological evidence for a protracted India-Eurasia collision.

## INTRODUCTION

The rapid plate motion of India toward Eurasia remains a major tectonic puzzle. Present-day plates move at rates less than 100 mm/year, but for a period of 20 million years (Ma) in the Cretaceous, the Indian plate moved at rates higher than 180 mm/year (14). Previous studies have primarily focused on explaining how subduction systems can exceed typical plate speeds (∼70 mm/year) (5, 6) and have overlooked the detailed shape of the India-Eurasia convergence history. Because the main driving forces, ridge push and slab pull, limit plate speeds (7), other mechanisms mechanisms such as plume push (4, 8) and double subduction (9, 5) have been explored. The concurrence of superfast spreading in the Indian Ocean at ∼65 Ma, with the maximum outpouring of Deccan flood basalts (10) and a slowdown of Africa, has led to the suggestion that India’s speedup was caused by the arrival of the Reunion plume at the base of the Indian lithosphere (4, 8). Numerical models indicate that the Reunion plume head could have contributed a few centimeters per year to the velocity increase of the Indian plate, but the effect is ephemeral and cannot explain the 20-Ma duration observed for fast convergence of India-Eurasia (8). The second hypothesis, double subduction, is based on the fact that the southern Eurasian margin is riddled with fragmented oceanic remains with different geological and geochemical signatures (5, 9, 1114). This indicates that multiple subduction systems operated within the Neo-Tethys since 130 Ma, when the Indian plate separated from Gondwana. Thus, fast Indian plate motion may result from the combined slab pull of two northward dipping subduction zones (15), and numerical models confirmed that same-polarity double subduction systems converge faster than single subduction systems (5, 6, 16), although none reported speeds in excess of 120 mm/year. Furthermore, previous studies do not address how double subduction systems initiate but rather begin with already mature systems.

We propose that these processes (plume push and double subduction) were coupled (illustrated in Fig. 1A), and detailed features of India-Eurasia’s convergence history (Fig. 1B) reflect a convolution of their effects. Specifically, we hypothesize that if the plume push force associated with the Reunion plume created sufficient compressional stress to initiate a subduction zone within the Neo-Tethys ocean (stage 1), then this incipient double subduction system could allow a brief episode of fast convergence with rates reaching ∼180 mm/year (stage 2). As the plume activity wanes, the system would transition (stage 3) to a mature double subduction system capable of sustaining high convergence rates for the observed duration of 20 Ma (stage 4) until the collisions of the intraoceanic arc (stage 5) and Indian continent (stage 6) occur.

## RESULTS

We compiled India-Eurasia convergence data from more than 30 published studies in the last 20 years that either produced new estimates or reinterpreted older data due to method improvements. We collapse the data in 58 plots, with a total of 70 unique convergence profiles (fig. S6). The most accurate recordings of India’s plate motion, shown in Fig. 2 (A and B), are from shiptrack data from the Central Indian Ridge (CIR; India-Africa) and Southeast Indian Ridge (SEIR; India-Antarctica) in the Indian Ocean (2, 4, 17). The detailed features of this high-resolution data occur consistently across methods and studies (fig. S7). We use this extended dataset to quantify the magnitude and durations of each of the proposed stages with our results shown as colored age intervals in Fig. 2 (A and B). Convergence data show that India moved steadily northeast relative to fixed Eurasia at 80 to 90 mm/year from ∼72 to 68 Ma, when it started accelerating to 180 mm/year at ∼66 Ma (1, 3, 8). This period includes increased spreading rates in the Indian Ocean (18), a corresponding slowdown of Africa (4), and culminates in a burst of superfast spreading [up to 220 mm/year (17)] from ∼66 to 63 Ma coincident with the time of maximum outpouring of Deccan flood basalts (10). The sudden drop to 120 mm/year at 63 to 62 Ma has been largely overlooked, despite being noted by some studies (17, 19, 20). From ∼62 to 50 Ma, a second period of faster spreading (up to 160 mm/year) precedes a prolonged slowdown from ∼50 to 45 Ma. The spreading azimuths at CIR and SEIR remain constant during this slowdown, until they abruptly change at ∼43 Ma (1, 2), indicating further changes in plate driving forces.

We test the hypothesized scenario with two-dimensional (2D) physically consistent numerical models of single and double subduction based on a previous model (16). Additional details of the numerical setup, parameters, and boundary conditions are given in Materials and Methods and Supplementary Materials (figs. S1 and S2). We consider a scenario in which an existing single subduction system (at the Eurasian margin initiated by at least 103 Ma) (21) is placed into compression by an increased convergence rate of the (left) Indian plate (mimicking the effect of plume push), which initiates an intraoceanic double subduction zone along a preexisting weak zone. The transition from stage 0 to 1 occurs by applying a specified velocity to the Indian plate, based upon previous models of subduction initiation (22, 23), which demonstrate that compressive stresses are necessary at initiating subduction zones. The preexisting weakness in our model can be interpreted as E-W oriented extensional faults created during seafloor spreading. This is based on recent observational and numerical work on intraoceanic subduction initiation mechanisms, which suggests inversion of spreading ridges to trenches (2426). In particular, extensional detachment faults can effectively localize deformation due to their weakness and control the nucleation of a subduction zone parallel to the former spreading center. However, other preexisting weaknesses could be crustal heterogeneities (i.e., microcontinents), transform faults, or lateral differences in the age and strength of the lithosphere (23, 27).

The results of our reference model are plotted with blue lines in Fig. 2 (A and B) and closely match both the magnitude and temporal variation of the observed convergence rates. The model produces the key features hypothesized in Fig. 1 including two distinct peaks of convergence rate with an intervening slowdown. Figure 3 shows the evolution of the reference model, at time steps corresponding to the colored markers in Fig. 2A. Stages 1 to 3 (Fig. 3A) represent a plume push–driven system in which the plume controls the dynamics of both plates. To initiate an intraoceanic subduction, the push must exceed the subduction rate of the right plate (16), as a lower intensity push would lead to a failed subduction scenario (fig. S3). After the intraoceanic subduction is initiated, the subduction of the left (Indian) plate dominates the system, which helps the plates remain attached to each other, and rapidly close the ocean basin (16). As plume activity wanes (stage 3), there is a considerable slowdown in velocities, marking the transition to a double subduction, with gradually increasing rates. Stages 4 to 6 (Fig. 3B) represent a double subduction–driven system with decreasing separation between the slabs causing an increase in dynamic pressure, until the flow is shut down in that area. The system slows down when the middle plate is consumed and arc-continent collision occurs. Stage 6 (Fig. 3B) shows a complex slab distribution in the mantle.

The reference model was chosen to best fit the available data in terms of amplitude and duration. However, changing the model parameters modifies the convergence profiles and aspects of subduction dynamics. For this reason, we analyze a suite of 40 simulations, with variations of mantle viscosity structure and plume push profiles (figs. S3 and S5). The abrupt change in the force balance (from dominantly plume push to sustained double subduction), marked by velocity drop, is a robust occurrence in all simulations despite various plume push magnitudes. However, the shapes of the two velocity peaks are strongly affected by the nature of the forcing and the viscosity structure, with only 50% of the double subduction systems exhibiting both peaks in excess of 100 mm/year.

Analysis of the entire suite of numerical simulations, shown in Fig. 2 (C and D), provides important distinctions between the shapes of these two peaks. First, we quantify the numerical results using a recently developed metric, the double subduction factor (DSF), which is the ratio of the volume of subducting slab on the left and the total volume subducted (Materials and Methods and figs. S4 and S5) (16). The DSF indicates the state of the system at a given time, with DSF = 0 and DSF = 1 for single subduction systems, while small DSF is representative of incipient double subduction systems, and DSF ∼ 0.55 indicates a mature double subduction system in steady state. The DSF is a useful metric, because the evolving force balance in a double subduction system results from the nonlinear interaction of subduction dynamics of multiple slabs and plates and is difficult to quantify. Second, we also normalize the convergence results by calculating the speedup in double subduction systems relative to single subduction rates, thus removing the dependency on model choices and extract the physics. This way, primary forces have the same magnitude in single and double subduction systems, and differences arise due to the double subduction character.

Time evolution data for all simulations are compressed and collapsed in Fig. 2 (C and D) using the two metrics (DSF and speedup). Results show that both forced and unforced (free) double subduction systems can achieve “fast velocities” (>100 mm/year) but not in the same stage of evolution. Forced subduction systems can achieve high convergence rates at all times in evolution (provided enough forcing), while unforced double subduction systems achieve high rates only at maturity (DSF ∼ 0.4 to 0.5). In addition, the speedup increase is approximately three times higher for forced convergence compared to unforced double subduction for the same stage of evolution. This suggests that the convergence signal given by the two processes are different: A forced system (plume push) gives a sharp (fast speedup) signal, while a free double subduction produces a gradual (slower speedup) convergence signal.

## DISCUSSION AND CONCLUSION

Analysis in Fig. 2D suggests that the sharp rise in velocities from 80 to 180 mm/year at ∼72 to 66 Ma is best explained by the effect of the Reunion plume, and not by an evolving double subduction as proposed before (5). Moreover, these anomalous rates are likely to induce subduction initiation, because models of successful subduction initiation require only an excess of 10 to 20 mm/year rates in the system (22). A key prediction of our hypothesis, however, is India’s deceleration at 63 to 62 Ma, which has been interpreted as a reorganization of plate boundary forces (19, 17). This predicted slowdown is a robust feature in the Indian Ocean spreading record, despite the signal being smoothed in the magnetic anomalies as oceanic crust forms and the inherent difficulty for such an event to be recorded by methods that use more indirect proxies (see Materials and Methods and the Supplementary Materials). Although we do not model the dynamical effect of plume-lithosphere interaction, we observe that the slowdown due to force balance redistribution is independent of the magnitude and duration of the forcing applied to the system (i.e., observed for lower intensity plume push).

We propose that the Reunion plume led to the formation of a double subduction system that kept the convergence rates high for another 10 to 15 Ma, until arc-continent collision and continental collision slowed down the entire system to present-day rates. The long-term India-Eurasia convergence consumed the equatorial Meso- and Neo-Tethyan ocean basins, leaving slab remnants in the upper- and mid-mantle imaged with seismic tomography (11, 2830). Subduction of two slabs could explain the seismic anomalies observed beneath the Indian Ocean distributed over a large area.

Our models also support a two-stage India-Eurasia collision model (5, 9, 15). This two-stage sequence starts with a “soft” collision at ∼50 Ma, perhaps the intraoceanic arc above the middle plate accreting onto the Eurasian margin, followed by a “hard” collision between continental India and Eurasia as early as 43 Ma (1, 2) or as late as ∼25 to 20 Ma (31, 32). The Kohistan-Ladakh Arc terranes in northwest Himalayas are a candidate for the intraoceanic arc (5), because some studies suggest a collisional age with the Lhasa terrane of 50 Ma (12). However, some evidence suggests that the Kohistan portion had accreted as early as ∼100 to 75 Ma (33), while others put the Kohistan arc being attached to Eurasia by 68 Ma (32). The western portion of Indus-Yarlung suture zone, which hosts a belt of Jurassic to Early Cretaceous ophiolites, may include other potential candidates for this island arc (34, 35). One such terrane could be the Xigaze ophiolite (36), derived from the Xigaze forearc located between the Gangdese arc (Eurasia) and the trench to the south, estimated to have subducted around 55 to 50 Ma ago.

Furthermore, our hypothesis is consistent with recent interpretations that Greater India was not an extension of the Indian continent but instead was a microcontinent separated from northern cratonic India by an ocean basin (31, 32). Numerical models show that subduction of large continental areas leads to slowdown or arrest of convergence/subduction (37), and this is inconsistent with subduction of a continental Greater India. Moreover, the locations and polarities of both subduction zones in our models are more consistent with resumed subduction south of the margin of Greater India.

In summary, we propose a logical sequence of events in the Neo-Tethys ocean basin that can explain the convergence history of India-Eurasia between 70 and 50 Ma. We hypothesize that the Reunion plume exerted a significant force at the base of the Indian lithosphere for a couple of million years, which likely induced initiation of a secondary subduction, and that the sudden drop in convergence history at 63 to 62 Ma is due to the waning of the plume. Remnants of closing a double subduction system may persist today, with India indenting into Eurasia at a high postcollisional rate of 50 mm/year and supporting the rise of the largest orogen on Earth (37).

## MATERIALS AND METHODS

### Numerical modeling of formation and evolution of double subduction

The equations of momentum and mass conservation (assuming incompressibility and neglecting thermal diffusion) are solved using the parallel 3D finite difference code LaMEM (Lithosphere and Mantle Evolution Model) (37, 38), capable of simulating lithospheric deformation while simultaneously taking mantle flow and an internal free surface into account. We use a pseudo-2D Cartesian domain in an approach similar to (16). The lithosphere and mantle material is assumed to behave as a continuous medium deforming by steady-state creep over long time intervals (39). This is expressed by a variable viscosity constitutive relationship

$τij=2ηε̇ij$

, where η is the Newtonian viscosity, constant for each material phase;

$ε̇ij=12(∂vi∂xj+∂vj∂xi)$

is the deviatoric strain rate tensor; and i and j represent spatial directions following the Einstein summation convention. A Lagrangian marker-in-cell method (40, 41) is used for accurately tracking distinct material domains (42), as they undergo extensive deformation due to creeping flow. We also use an internal free surface, using the “sticky-air” approach, with a free surface stabilization algorithm (43) that allows the development of topography.

### Model setup

We performed 40 2D numerical simulations of single subduction and double subduction system with influx boundary conditions (table S1 and fig. S1). The model domain is 6000 km across and 2048-km deep. The computational domain has a variable grid spacing, with higher resolution in the upper mantle asthenosphere and close to the subduction trenches (minimum and maximum grid spacings: [xmin, xmax] = 3 km, 12.5 km, [zmin, zmax] = 5 km, 22.6 km). Free-slip boundary conditions are imposed on all boundaries and sticky-air on top of the plates.

The initial model setup and material parameters are similar to the ones used in (16) (fig. S1A). The model consists of a large oceanic plate (4000 km long) subducting beneath an upper plate (2000 km long). The lithosphere has a thickness of 80 km with a 20-km-thick core and 15-km weak crust. While the average ocean crust is only about half as thick, we consider this a parameterization of strength weakening with depth due to the sediment layer and hydrated oceanic lithosphere. The sticky-air layer has a height of 60 km. A 20-km-thick weak zone dipping 30° in the lithospheric mantle that mimics a weak subduction fault and extends to 110-km depth is introduced in double subduction simulations in the middle of the oceanic plate at x = 2000 km. Single subduction simulations lack this weak zone. Material parameters are the ones used in (16), where the mantle asthenosphere has a reference density (ρ0) and viscosity (η0), and the plates are 85 kg/m3 denser and have a variable viscosity structure (slab, 500 × η0; strong core, 5000 × η0; weak crust, η0; and weak zone, η0).

The model setup is built in such a way to resemble the closure of the Neo-Tethys during the India-Asia convergence. Thus, the left subduction system is analogous to the Indo-Australian plate, the middle plate (also referred to as the right plate) is an oceanic plate, and the upper plate represents the Eurasian plate. With this configuration, we assume a long-term subduction at the Eurasian plate margin, consistent with tectonic reconstructions (11). We start from a single subduction setup, where subduction is already initiated (mature), and we stress the system by controlling the convergence rate of the left subduction (i.e., imposing influx/outflux boundary conditions).

### Mantle density and viscosity profiles

The mantle viscosity structure represents one of the largest uncertainty in our models. For this reason, we test four different mantle viscosity profiles and two different density profiles (fig. S2) that affect the sinking velocity of slabs. The viscosity profiles are as follows: a mantle transition zone at 660 km with a viscosity jump of 50 (TZ 660 km) (44) and three profiles derived from (45) (RLLB15). The RLLB15 viscosity profiles are discretized in 10 layers of 200 km thickness. The mantle density profiles are used to calibrate the rheological profiles (i.e., control sinking velocities) in our models: (i) constant density with depth and (ii) density increase of 25 kg/m3 per 200 km. A combination of eight different mantle and density profiles is then obtained (i.e., fig. S3). Model results do not depend on radial viscosity structure.

### Influx/outflux boundary conditions

The effect of plume push is implemented in LaMEM as influx/outflux boundary conditions on the left boundary (i.e., injecting oceanic lithosphere material with a given velocity, Vinflux). To conserve mass, an outflux velocity is calculated and imposed on the mantle part of the left boundary to balance the injected material. Imposing velocity boundary conditions, instead of flux boundary condition, is motivated by previous approaches of subduction initiation studies (23). The nature and the effect of the plume head force are poorly constrained, so we test various magnitudes of the forcing for 9 Ma (controlled convergence stage; fig. S1B), after which we let the system evolve dynamically (free subduction stage with free-slip boundary conditions on the left boundary). The behavior of plume push is supposed to be short-lived high-intensity forcing, so we impose forced convergence for 9 Ma, which also fits the time interval for the first peak in the India-Eurasia convergence data (72 to 63 Ma). Single subduction simulations, used primarily to normalize the results in double subduction simulations, do not have influx/outflux boundary conditions (i.e., free subduction). When calculating speedup, maximum rates in single subduction models are used to normalize double subduction results, thus putting a lower bound on the speedups obtained.

### Double subduction factor

The DSF was first introduced in (16) to better understand the dynamics of double subduction (fig. S4). The DSF is intended to be a measure for slab interaction and is the ratio of the volume of subducting slab on the left and the total volume subducted:

$DSF=VslableftVslableft+Vslabright$

, where

$Vslableft$

is the volume of subducting slab on the left, and

$Vslabright$

is the volume of the slab subducting on the right. The DSF is primarily an indicator of the state of the system at any point in time. If DSF = 0, the system is entirely driven by the right subduction, while a DSF = 1 would mean that the system is driven entirely by the pull of the left slab. These two end-member values actually represent single subduction systems. Values between 0 and 1 are indicative of a double subduction system, and when DSF = 0.5, the pull of both slabs is equal.

The DSF also reflects the evolution of a double subduction from incipient to mature system. For example, when a double subduction system fails to initiate, the majority of the work is done by the right slab and DSF ∼0. An increasing DSF value reflects the participation of the left subduction to the dynamics of the system. Thus, states where DSF ∼0.5 are considered mature double subduction configurations, while states with DSF values closer to 0 and 1 represent incipient double subduction configurations. A result of this is that successful initiation of double subduction can be quantitatively confirmed by the DSF (i.e., 0 < < DSF < <1). Pusok and Stegman (16) found that the mature steady state of a double subduction is when DSF ∼ 0.55 (i.e., slightly more pull from the left slab for the configuration shown in fig. S4). This imbalance is due to the asymmetry of the same-dip double subduction system, where the right plate also acts as the overriding plate for the left slab. This coupling provides an extra back resistive force for the right subduction and a lower overriding plate load for the left subduction. Thus, the right subduction is gradually slowed down by the motion of the left slab, compared to if it was subducting alone (i.e., single subduction). In the simulations presented here (i.e., fig. S5), the DSF either stays close to 0 (single or failed double subduction) or generally increases during forced convergence, as the system transitions from single subduction to stable double subduction. We collapsed the DSF and convergence evolution shown in fig. S5 as average points and variability intervals in Fig. 2 (C and D) for all simulations.

### Compilation of published data on India-Eurasia convergence

We compiled India-Eurasia convergence data published in the last 20 years that either produced new estimates or reinterpreted older data due to method improvements. We obtain a dataset of unique 58 profiles of convergence histories (Ind-Ant, Ind-Afr, Ind-Abs, and Ind-Eur) originating from more than 30 studies (fig. S6). The data were digitized using WebPlotDigitizer (https://automeris.io/WebPlotDigitizer), and the full list of references used is provided in the Supplementary Materials. We primarily focused on recent published data, assuming that the methods have improved over the years (i.e., accuracy is higher and design errors are smaller compared to earlier studies). We do, however, acknowledge the valuable pre-2000 work and include a couple of influential studies [i.e., (1)] that constitute the basis of later studies, but do not modify our observations here.

This is the largest effort to compile such a dataset. Previous reconstructions primarily differ due to (i) focus (motion of India relative to surrounding plates or absolute motion of India relative to mantle) and (ii) method (i.e., Indian Ocean spreading, plate circuit, paleolatitude, and reinterpretations) and variations in method design (i.e., different study areas or inversion algorithms). Then, we have secondary factors that can produce variations in the convergence data: different time focus and resolutions (some studies focus on last 50 Ma, others last 90 Ma, while others go beyond 200 Ma), different parameters and units (millimeters per year or centimeters per year, half-spreading rates), and improvements in the geomagnetic time scales [main ones: CK95 (46), GTS04 (47), and GTS12 (48)].

We thus collapse the data in 58 plots, with a total of 70 unique convergence profiles (fig. S6, 58 + 12 additional profiles plotted in dotted red lines). Of these, 11 profiles also provide confidence intervals or E/W Himalayan syntaxes estimates (green and gray error bars). Each plot in fig. S6 highlights with shaded colors the proposed time intervals in Figs. 1 and 2. These time intervals (or evolution stages) will constitute the basis of our analysis in figs. S7 and S8.

We divide previous reconstructions in six categories based on the primary method: (A) SEIR spreading rates, which give the relative motion between India/Capricorn-Antarctica (9 plots); (B) CIR spreading rates, which give the relative motion between India/Capricorn-Africa/Somalia (8 plots); (C) estimates of the absolute motion of India relative to the mantle derived from the Indian Ocean spreading data (India-Abs, 3 plots); (D) paleolatitude data, which provide estimates on the India-Eurasia convergence (4 plots); (E) plate circuit (India-Eurasia, 19 plots); and (F) plate reconstruction models, which provide various relative or absolute plate motions derived from a combination of observations and computer models (15 plots). Of these, Indian Ocean spreading methods (A to C) generally calculate plate rotations using geophysical data (magnetic anomalies, satellite gravity, fracture zones, synthetic flowlines, or ship lines) and provide the highest spatial and temporal resolution with high confidence intervals (95%). However, these high-resolution estimates are for the motion of India relative to Africa or Antarctica and not necessarily India-Eurasia. Plate circuit methods (E) resolve this issue by using rotation poles from the Indian Ocean plate configuration together with other plates in a global circuit, which results in the net convergence between India-Eurasia, but at the expense of temporal resolution and accuracy. Last, plate reconstruction models combine various observations (i.e., geophysical, geological, and hot spot motions) with computer models, and a prominent software in the last years to provide these types of estimate is GPlates software (http://www.gplates.org) (49).

More details on each of these methods can be found in the individual publications from which these datasets were collected. What we want to highlight in this compilation is that some features of the convergence data are consistent across methods and studies. Figure S7 shows the maximum, mean, and minimum values extracted for every convergence curve (sample) in different evolution stages between 0 and 72 Ma (1, orange; 2, red; 3, yellow; 4, green; 5, blue; 6, gray) and then arranged for every method (A to F). The total number of samples is N = 91 (70 unique profiles plus 22 confidence intervals, and subtracting 1 from plot 49, which only has estimates from E/W Himalayan syntaxes and no mean). The black lines in fig. S7 represent the means, and overall have similar amplitudes across methods, except for paleolatitude estimates (D), which deviate more from the general pattern. The SEIR and CIR spreading data (A and B) are remarkably consistent across studies (with small variations in amplitude and duration) and stages (i.e., stages 3, 4, and 6). The majority of them predict the slowdown at 62 Ma, with incredibly small variations in both the minimum and maximum values (i.e., stage 3).

Figure S8 shows the average convergence model for India for every method (A to F, panels a to f), an overall average (panel g), and a best average (panel h). They all predict a similar pattern: variable orange stage with a high mean—suggesting rapidly changing rates, anomalously high red stage, sharp yellow slowdown, considerable high green rates, first blue slowdown, and second gray slowdown to present-day values. This analysis suggests a general pattern of the convergence data: double high peaks, with a short and sudden drop in between, followed by two slowdowns. We hypothesize that the peaks correspond to a plume push stage and a double subduction stage, with a reorganization of the force balance in between, followed by a two-stage collision, first arc-continent collision and second continental collision. Our preferred convergence model (H; best average) was chosen in the following way: (i) the orange stage represents an increase in velocity, so the average between all methods is more representative; (ii) the peak plume stage (red) uses the data that have the highest resolution and maximum values; (iii) the slowdown at 62 Ma (yellow stage) also uses the highest resolution data, but minimum records registered in that time interval; (iv) the double subduction stage (green) highlights the maximum values calculated; and the last two stages (blue and gray) show averages for the respective time periods.

We note that, for this analysis, we are not in a position to say that one method is better than the other. For example, we tried to reduce the bias due to a particular method, but if different studies used the same method, the results will be biased toward the results of that particular method. Every method has advantages and disadvantages, and we hope that this will motivate future studies to improve both spatial and temporal resolutions. For example, spreading rates in the Indian Ocean (SEIR and CIR) have a higher temporal and spatial resolution and more precise methods, but they need to be combined with equally high quality data on other plates in plate circuit methods to obtain the desired India-Eurasia convergence history. The following references were used in convergence data compilation (figs. S6 to S8): (14, 8, 1721, 3032, 5065). Note that a newly published study (66) that was not included in the data compilation also shows a two-peak convergence profile and a sharp slowdown at 62 Ma.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

1. D. Turcotte, G. Schubert, Geodynamics (Cambridge Univ. Press, ed. 3, 2014).

2. T. Gerya, Introduction to Numerical Geodynamic Modelling (Cambridge Univ. Press, ed. 1, 2009).

3. F. Gradstein, A Geologic Time Scale (Cambridge Univ. Press, 2004).

4. J. G. Ogg, Geomagnetic polarity time scale, in The Geologic Time Scale 2012, F. M. Gradstein, J. G. Ogg, M. D. Schmitz, G. Ogg, Eds. (Elsevier Science, 2012), pp. 85–112.

5. M. Searle, in Treatise on Geophysics (Second Edition), G. Schubert, Ed. (Elsevier, 2015), pp. 469–511.

Acknowledgments: We thank D. Müller and an anonymous reviewer for their helpful comments on the manuscript. Funding: A.E.P. is grateful for postdoctoral research funding from the Cecil H. and Ida M. Green Foundation, UCSD and support from the NSF (grant EAR-1255040). The numerical simulations were performed on Comet (San Diego Supercomputing Center) within the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF grant ACI-1053575. Author contributions: Both authors discussed and developed the central ideas and contributed to the writing of the manuscript. A.E.P. carried out the numerical modeling and analysis. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The Bitbucket version of the numerical code (LaMEM) used can be found in https://bitbucket.org/bkaus/lamem/branch/cvi_test, and the repository containing the input parameters files to reproduce the data can be found in https://bitbucket.org/adina/rep-mconvindia. The full simulation data presented in this study can be provided on request.